7

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.

  • I think the terms you should search on are fishnet and vector grid. – PolyGeo Jan 19 '17 at 19:44
  • Thanks @PolyGeo, however those search terms appear quite closely tied in with ArcGIS terminology, while they are what I want- I am looking for a way to do this in R. – Amroco Jan 19 '17 at 20:05
  • 2
    Vector grid is a generic term favored by QGIS users whereas ArcGIS users tend to use fishnet. Both are stored data whereas "grids and graticules" are more often generated only for display (on-the-fly). – PolyGeo Jan 19 '17 at 20:14
  • The term grid is very ambiguous within GIS so much care should be taken when discussing and searching on it. I am confident that what you are describing, and that other Q&A gives a solution to, is how to create a vector grid using R. Whether that vector grid can be polygons or just lines I cannot be certain but, based on that Q&A, it would seem to make this a duplicate. – PolyGeo Jan 19 '17 at 20:22
  • 1
    From the linked answer, just changing the call for coerce to as(r, "SpatialPolygonsDataFrame"), and saving it by rgdal::writeOGR() gives me a polygon file. Is it what you need? – Kazuhito Jan 19 '17 at 20:33
2

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.

Your Answer

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

Not the answer you're looking for? Browse other questions tagged or ask your own question.