Wednesday 12 April 2017

polygon - Generate rectangular fishnet or vector (grid cells) shapefile in R?


I would like to create a polygon of 20 x 10 degree graticules from 180, -180 to 90, -90 degrees.



To do this in QGIS, I followed:
Vector creation tools > Create grid
[Grid extent: -180, 180, -90, 90]
[Horizontal spacing: 20, vertical spacing: 10]
And selected 'rectangular polygon' rather than 'rectangular grid'.


I got the purple grid/vector which was 18 rectangles across (360/20) in the picture below. enter image description here


Purple: grid produced by QGIS, 18 rectangles across; Black: grid by code below, 19 rectangles across.


I then tried to replicate this grid using R. Similar to Generating grid shapefile in R?, (however, covering the whole extent of long and lat, and resulting in a rectangular grid).


I used:


    r <- raster(extent(matrix( c(-180, -90, 180,  90), nrow=2)), nrow=20, ncol=10, 

crs = "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0")
r[] <- 1:ncell(r)
summary(r)
print(r)
plot(r)

I also need the product as a polygon class and not spatial lines or grid.


    r2 = as(r, "SpatialPolygonsDataFrame")
class(r2)
writeOGR(r2, dsn=getwd(), layer="trial", driver="ESRI Shapefile",

overwrite_layer=T)

I shouldn't end up with 19 rectangles across if this code was using degrees, so I am not sure how that code works, or where I have gone wrong.


I also need this as a polygon, not lines, so tried:


grid_20_line = SpatialLines2PolySet(grid_20)
grid_20_poly = PolySet2SpatialPolygons(grid_20_line)

Many square vector's (such as 10 x 10 degree cells) like this can be downloaded, but I would like to know the options for making rectangles in R for specific degrees.


And to do this in QGIS is very simple, so I am curious if there is an R solution.



Answer




Your raster-sp workflow worked fine for me, not sure how you got 19 columns either. If it wasn't a typo, it might have been an rgdal bug on write; that's the only related package that I can see a recent update for. Anyway here's a couple of examples with the newer sf package:


library(sf)

# make an object the size and shape of the output you want
globe_bb <- matrix(c(-180, 90,
180, 90,
180, -90,
-180, -90,
-180, 90), byrow = TRUE, ncol = 2) %>%
list() %>%

st_polygon() %>%
st_sfc(., crs = 4326)

# generate grid of 20 x 10 tiles (200 in total, each 18 x 18 degrees)
globe_grid_18x18 <- st_make_grid(globe_bb, n = c(20, 10),
crs = 4326, what = 'polygons') %>%
st_sf('geometry' = ., data.frame('ID' = 1:length(.)))

# identical to the above:
globe_grid_18x18a <- st_make_grid(globe_bb, cellsize = c(18, 18),

crs = 4326, what = 'polygons') %>%
st_sf('geometry' = ., data.frame('ID' = 1:length(.)))

# you wanted this though, yes? Grid of 20 x 10 tiles (324 tiles, each 20 x 10 degrees)
globe_grid_20x10 <- st_make_grid(globe_bb, n = c(18, 18),
crs = 4326, what = 'polygons') %>%
st_sf('geometry' = ., data.frame('ID' = 1:length(.)))

# e.g. for checking in qgis
st_write(globe_grid_20x10, 'C:/DATA/globe_grid_20x10_sf.gpkg')


The st_sf() step is optional, if you just want to write the grid directly to file.


No comments:

Post a Comment

arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...