How can I rasterize a SpatialPolygonsDataFrame so that I can fill a shapefile with points?

Question

I have a chunk of R code that used to work but does not work anymore and I cannot find the issue. The purpose of the code is to fill a shapefile with regularly spaced points.

My shapefile can be accessed here: https://drive.google.com/drive/folders/1SAbuyIQHevK4fz-0w3TTqpEhz0wKLEII?usp=sharing

If I begin with loading my shapefile:

GUA = raster::shapefile('Guam3BufferPoly.shp')

Then I set a variable for the coordinate reference system for this SpatialPolygonDataFrame:

projGUA = crs(GUA)

Transform to planar crs

putm <- spTransform(GUA, projGUA)

Create a raster (this is where it doesn't work)

ext = extent(putm)
r <- raster(ext, res=500) 

Rasterize the polygon and transform to points

r2 <- rasterize(putm, r)
pts <- rasterToPoints(r2, spatial=TRUE)

Transform the points to lon/lat and plot the results

pts_lonlat <- spTransform(pts, "+proj=longlat +datum=WGS84")
plot(pts_lonlat,pch='*') 

The raster, r, is empty (breaking all the code downstream).

Please let me know if you can help me. And please be kind (this is my first time posting here and I apologize if I have not formatted my question correctly). Thank you!

Answer 1

Welcome!

Basically, it seems like you're working in a coordinate reference system not appropriate for the steps you execute.

library(raster)

GUA <- shapefile("Guam3BufferPoly.shp")

crs(GUA)
#> Coordinate Reference System:
#> Deprecated Proj.4 representation: +proj=longlat +datum=WGS84 +no_defs 
#> WKT2 2019 representation:
#> GEOGCRS["GCS_unknown",
#>     DATUM["World Geodetic System 1984",
#>         ELLIPSOID["WGS 84",6378137,298.257223563,
#>             LENGTHUNIT["metre",1]],
#>         ID["EPSG",6326]],
#>     PRIMEM["Greenwich",0,
#>         ANGLEUNIT["Degree",0.0174532925199433]],
#>     CS[ellipsoidal,2],
#>         AXIS["longitude",east,
#>             ORDER[1],
#>             ANGLEUNIT["Degree",0.0174532925199433]],
#>         AXIS["latitude",north,
#>             ORDER[2],
#>             ANGLEUNIT["Degree",0.0174532925199433]]]

Seems like the crs provided Guam3BufferPoly.prj is not quite up to standard, but I'm not 100 % sure here. However, you're working with lat/lon coordinates on WGS 84. Hence, spTransform(GUA, projGUA) has absolutely no effect, since you're reprojecting your feature to the same crs again.

This leads to a follow-up issue: You're working in degrees as units and specify a raster of 500 units (here: degrees, not meters, as probably intended). And this is exactly what raster() does, even if it makes no sense:

ext <- extent(GUA)
ext
#> class      : Extent 
#> xmin       : 144.5962 
#> xmax       : 145 
#> ymin       : 13.19925 
#> ymax       : 13.70405

r <- raster(ext, res = 500)
r
#> class      : RasterLayer 
#> dimensions : 1, 1, 1  (nrow, ncol, ncell)
#> resolution : 500, 500  (x, y)
#> extent     : 144.5962, 644.5962, -486.2959, 13.70405  (xmin, xmax, ymin, ymax)
#> crs        : NA

To solve this, you should use a projected coordinate reference system suitable for Guam. Not sure how good of a fit this is, you'll probably know better, but let me try WGS 84 / UTM 55 S (EPSG: 32755):

putm <- spTransform(GUA, CRS("+proj=utm +zone=55 +south +datum=WGS84 +units=m +no_defs"))

ext <- extent(putm)
ext
#> class      : Extent 
#> xmin       : 239608.8 
#> xmax       : 283618.4 
#> ymin       : 11460271 
#> ymax       : 11516045

r <- raster(ext, res = 500)
r
#> class      : RasterLayer 
#> dimensions : 112, 88, 9856  (nrow, ncol, ncell)
#> resolution : 500, 500  (x, y)
#> extent     : 239608.8, 283608.8, 11460045, 11516045  (xmin, xmax, ymin, ymax)
#> crs        : NA

r2 <- rasterize(putm, r)
pts <- rasterToPoints(r2, spatial = TRUE)

pts_lonlat <- spTransform(pts, "+proj=longlat +datum=WGS84")
plot(pts_lonlat, pch='*') 

Edit:

You can simplify the process using sf instead of raster and sp IMHO. Individual steps include data import, reprojection to a metric crs, grid creation covering the whole extent of your dataset and "spatial filtering" by intersecting your buffer with the grid centroids.

library(sf)

sf <- st_read("Guam3BufferPoly.shp", crs = "epsg:4326") |> st_transform("epsg:32755")
#> Reading layer `Guam3BufferPoly' from data source 
#>   `Guam3BufferPoly.shp' using driver `ESRI Shapefile'
#> Simple feature collection with 1 feature and 1 field
#> Geometry type: POLYGON
#> Dimension:     XY
#> Bounding box:  xmin: 144.5962 ymin: 13.19925 xmax: 145 ymax: 13.70405
#> Geodetic CRS:  WGS 84

grid <- st_make_grid(sf, cellsize = 500, what = "centers") |> st_intersection(x = sf, y = _) 

grid
#> Simple feature collection with 3839 features and 1 field
#> Geometry type: POINT
#> Dimension:     XY
#> Bounding box:  xmin: 239858.8 ymin: 11460520 xmax: 283358.8 ymax: 11516020
#> Projected CRS: WGS 84 / UTM zone 55S
#> First 10 features:
#>     ID                  geometry
#> 1    1 POINT (246358.8 11460521)
#> 1.1  1 POINT (246858.8 11460521)
#> 1.2  1 POINT (247358.8 11460521)
#> 1.3  1 POINT (247858.8 11460521)
#> 1.4  1 POINT (248358.8 11460521)
#> 1.5  1 POINT (248858.8 11460521)
#> 1.6  1 POINT (249358.8 11460521)
#> 1.7  1 POINT (249858.8 11460521)
#> 1.8  1 POINT (250358.8 11460521)
#> 1.9  1 POINT (250858.8 11460521)

^{Created on 2022-08-18 by the reprex package (v2.0.1)}