Deutsches R-Forum

Verfasst: **Do Aug 17, 2017 11:53 am**

Hallo (:

nachdem mir hier das letzte Mal bereits gut geholfen wurde bin ich in meinem Projekt inzwischen etwas weiter. Dabei ist eine neue Frage aufgekommen an der ich ein wenig fest hänge, war mir allerdings nicht sicher in welchem Teil des Forums ich es am besten poste. Zunächst einmal habe ich einen kleinen veränderten Beispieldatensatz. Es handelt sich dabei um Messwerte mit Koordinaten der jeweiligen Messstelle.

Code: Alles auswählen

dput(t1)
structure(list(Lfd = 1:20, MstNr = structure(c(16L, 17L, 15L, 
14L, 8L, 7L, 20L, 6L, 5L, 12L, 19L, 9L, 11L, 10L, 4L, 1L, 2L, 
18L, 13L, 3L), .Label = c("3369B", "3370", "3376", "3377", "3378", 
"3379", "3380", "3381", "3382", "3383", "3384", "3385", "3404", 
"3445", "3446", "3452", "3456", "636", "850", "851"), class = "factor"), 
    x = c(3469489.45369692, 3469083.73148078, 3469041.96280914, 
    3468635.21639453, 3468449.03013263, 3468450.69908635, 3469320.1627802, 
    3469424.6028951, 3469282.52506032, 3468388.84707303, 3469333.39947069, 
    3468243.43041729, 3469290.56665505, 3469588.87341771, 3469530.48687964, 
    3469971.09926508, 3469983.19236769, 3469736.45998308, 3470503.41533231, 
    3469088.22260067), y = c(5413581.63214282, 5413280.51059682, 
    5412268.69015278, 5412768.89058611, 5413137.30120601, 5412192.65287435, 
    5412100.59727255, 5414934.89892535, 5411931.3317301, 5413421.79058764, 
    5412878.37172573, 5413334.97727913, 5412852.64747422, 5413525.19319656, 
    5415130.55043546, 5412607.6942914, 5413845.45738916, 5412053.95384853, 
    5413575.24603885, 5415347.22405993), Messwert = c(0.333433333333333, 
    0.333433333333333, 0.334108333333333, 0.333433333333333, 
    0.356583333333333, 0.334401666666666, 0.351333333333333, 
    0.389209649122807, 1.15468717948718, 0.366882745098039, 0.443333333333333, 
    0.368706060606061, 0.35508696969697, 0.430333333333333, 0.497765151515151, 
    0.334573333333333, 0.366733833333333, 0.336677333333333, 
    0.392514583333333, 0.382595913978494)), .Names = c("Lfd", 
"MstNr", "x", "y", "Messwert"), class = "data.frame", row.names = c("1", 
"2", "3", "4", "5", "6", "7", "8", "9", "10", "11", "12", "13", 
"14", "15", "16", "17", "18", "19", "20"))

Was ich tun möchte ist diese auf ein Raster zu plotten und dieses Raster um die Punkte herum abzuschneiden. Den ersten Teil bekomme ich auch ganz gut hin denke ich:

Code: Alles auswählen

#Pakete 

install.packages("sp")
library(sp)

##Raumbezug herstellen
coordinates(t1) <- c("x", "y") 	# x/y in Koordinaten umwandeln
plot(t1)			# testweise plotten


###Raster dimensionieren

min_x = min(t1$x) - 100			# Min x + innerer Rand
min_y = min(t1$y) - 100			# Min y + innerer Rand
max_x = max(t1$x) 			# Max x			
max_y = max(t1$y) 			# Max y
x_length = max(t1$x - min_x) + 100	#  x-Länge + äußerer Rand
y_length = max(t1$y - min_y) + 100	#  y-Länge + äußerer Rand
cellsize = 100				# Zellgröße
ncol = round(x_length/cellsize,0)	# Zeilen
nrow = round(y_length/cellsize,0)	# Reihen

proj4string(t1) = CRS("+proj=utm +ellps=WGS84 +datum=WGS84")	# Projektion

#### Grid erstellen und Klasse zuweisen
grid = GridTopology(cellcentre.offset=c(min_x,min_y),cellsize=c(cellsize,cellsize),cells.dim=c(ncol,nrow))
grid = SpatialPixelsDataFrame(grid,
                              data=data.frame(id=1:prod(ncol,nrow)),
                              proj4string=CRS(proj4string(t1)))
                              
##### Punkte auf Grid plotten

plot(grid)
points(t1)

Leider hab ich keine wirklich Vorstellung wie ich die überschüssigen Teile des Grids nun abschneide. Ich hätte vermutet das man zunächst ein Polygon um die Punkte erstellt welches immer die äußersten Punkte verbindet. Anschließend könnte man das Polygon dazu benutzen die Fläche aus dem Raster zu clipen. Ich würde mich über jede Hilfe freuen (:

mit freundlichen Grüßen

Johann

Verfasst: **Do Aug 17, 2017 11:53 pm**

Andere Idee:

Das raster direkt auf die richtige Größe (=bounding box + Rand) erzeugen:

Code: Alles auswählen

library(sp)
library(raster)
# make spatial
t1_sp <- t1 
coordinates(t1_sp) <- c("x", "y")
proj4string(t1_sp) = CRS("+proj=utm +ellps=WGS84 +datum=WGS84")	
plot(t1_sp)


# get bounding box of point & ad margin of 100
ext <- bbox(t1_sp)
ext[ , 1] <- ext[ , 1] - 100
ext[ , 2] <- ext[ , 2] + 100

# create raster from bounding box
r <- raster(extent(ext), 
            crs =  proj4string(t1_sp),
            resolution=100)
r
# add some random noise (just for color...)
r[] <- rnorm(length(r))

plot(r)
points(t1_sp)

Verfasst: **Sa Aug 19, 2017 12:35 pm**

Hi EDI,

danke erstmal für den Vorschlag! Allerdings bekomme ich wenn ich das umsetze noch immer ein rechteckiges Grid. Ich stelle es mir aber eher so vor dass das Grid nur die Punkte umgibt, als würde man eine Linie um die Punkte herum zeichnen.

Vlt. kann ich es mit einem Beispiel etwas verdeutlichen, auch wenn hier teilweise ein etwas größerer Rand besteht und ich vermute dass die Kontur etwas anders zustande gekommen ist:

Code: Alles auswählen

library(gstat)

data(meuse.grid)
coordinates(meuse.grid) <- c("x", "y")
meuse.grid <- as(meuse.grid, "SpatialPixelsDataFrame")
meuse.grid
plot(meuse.grid, col="yellow")

data(meuse)
coordinates(meuse) <- c("x", "y")
points(meuse$x, meuse$y)

Verfasst: **So Aug 20, 2017 12:10 am**

Dann eine convex hull nehmen (das Polygon, dass durch die äußeren Punkte aufgespannt wird) und dann zuschneiden.

Code: Alles auswählen

library(raster)
library(tidyverse)

# an empty raster of global extent
r <- raster()
# make up some raster values
values(r) <- runif(length(r))

# make up some random coordinates around north america
coords <- cbind( 
  lon = runif(100, min = -120, max = -60),
  lat = runif(100, min = 30, max = 50))

# let's have a look
plot(r)
points(coords, add = TRUE)

# get the convex hull
hull_points <- coords[chull(coords),]

# convert to polygon
hull_polygon <- hull_points %>%
  Polygon() %>%
  list() %>%
  Polygons(1) %>%
  list() %>%
  SpatialPolygons()

# mask the raster
rr <- mask(r, hull_polygon)

# let's have another look
plot(rr)

Quelle:
https://stackoverflow.com/questions/452 ... onvex-hull

Verfasst: **So Aug 20, 2017 8:47 am**

Genial, genau das was ich gesucht hab

danke EDi, auch für den Link!

Deutsches R-Forum

Grid zuschneiden

Grid zuschneiden

Re: Grid zuschneiden

Re: Grid zuschneiden

Re: Grid zuschneiden

Re: Grid zuschneiden