2

I'm writing R script which will allow users to overlay a grid on their data in a new projection with the option of defining the size of the grid cell. This seems to work fine, however, when I plot two grids with different cell sizes over top of each other, the grids do not align (see photo). Any suggestions on how to standardize the position of the grid so grids with different cell sizes will align with each other?

new projection to use

projection <- "+proj=aea +lat_1=-30 +lat_2=-50 +lat=-40 +lon_0=175 +x_0=0 +y_0=0 +ellps=WGS84 +units=m +no_defs"

select size of cell to use

size <- c(250, 1000, 2000, 4000, 8000, 16000, 32000, 64000)

grid.size <- size[6]

read in data

data.df <- "your data set"

turn data into SpatialPoints

data.df <- SpatialPoints(cbind(data$lon, data$lat), proj4string = CRS("+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0")) 

lat/long projection

print(extent(data.df))

plot data in lat/long coordinates

plot(coordinates(data.df))

change lat/long projection to new projection and plot

data.t <- spTransform(data.df, CRS(projection))

data.t.coords <- coordinates(data.t)

plot(data.t.coords, xlim=c((extent(data.t)[1]*1.1), (extent(data.t)[2]*1.1)), ylim=c((extent(data.t)[3]*1.1), (extent(data.t)[4]*0.9)))

plot(data.t.coords, xlim=c((extent(data.t)[1]), (extent(data.t)[2])), ylim=c((extent(data.t)[3]), (extent(data.t)[4])))

define lon0, lon1, lat0, lat1

minx <- min(data.t.coords[,1]) #lon0
maxx <- max(data.t.coords[,1]) #lon1
miny <- min(data.t.coords[,2]) #lat0
maxy <- max(data.t.coords[,2]) #lat1

print(c(minx, maxx, miny, maxy, grid.size))

calculate number of cells in each dimension

Lx <- (maxx - minx)

Ly <- (maxy - miny)

Nbx <- (round((Lx) / grid.size)) #* 1.1 #add 10% buffer to grid

Nby <- (round((Ly) / grid.size)) #* 1.1

build grid using centroid of bottom left cell, cell length and width (equal), number of cells in each dimension; spatialise grid w/new projection

grid <- GridTopology(c(minx, miny), c(grid.size, grid.size), c(Nbx, Nby))

grid <- SpatialGrid(grid, proj4string = CRS(projection))

grid <- as(grid,'SpatialPixels')

grid <- as(grid,'SpatialPolygons')

minx <- grid@bbox[1,1]

maxx <- grid@bbox[1,2]

miny <- grid@bbox[2,1]

maxy <- grid@bbox[2,2]

print(c(minx, maxx, miny, maxy))

plot grid

for(i in 1:length(grid)){
  for(po in 1:length(grid@polygons[[i]]@Polygons)){
    pol <- grid@polygons[[i]]@Polygons[[po]]
    xp <- fortify(pol)[,1]
    yp <- fortify(pol)[,2]
    polygon(xp, yp, col=rgb(0, 0, 0, 0), border='red')
  }
}

enter image description here

  • there's a function in spex to buffer an extent to a given grain, create a raster extent (much easier than working with GridTopology btw) and create an extent from your data with buffer 64000, all your other size values will align if you use that master extent e.g. ex <- spex::buffer_extent(raster::extent(3383838, 4944040, -10, 3838378), 64000) (I just use random values in there for xmin, xmax, ymin, ymax) – mdsumner Nov 11 '19 at 0:05

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.