4

I'm having an issue rasterizing a shapefile to produce points on a 0.5*0.5 grid. The shapefile represents classifications of risk level (Low-0, Medium-100, High-1000, Very High-1500) of global coral reefs to integrated threats.

I pulled the code from another example that works fine, but when I try it for my data I get nothing from the plot function. See below for the link to the shapefile and my code:

Reefs At Risk: Global Integreated Threats

# Read shapefile into R
library(rgdal)
library(raster)    

int.threat.2030 <- readOGR(dsn = "Global_Threats/Integrated_Future", 
                           layer = "rf_int_2030_poly")

## Set up a raster "template" for a 0.5 degree grid
ext <- extent(-110, -50, 0, 35)
gridsize <- 0.5
r <- raster(ext, res=gridsize)

## Rasterize the shapefile
rr <- rasterize(int.threat.2030, r)

## Plot raster
plot(rr)

Any ideas where I might be going wrong? Is it an issue with the shapefile itself?

Please and thanks!

Robert Hijmans
  • 40,301
  • 4
  • 55
  • 63
RMFish
  • 43
  • 1
  • 1
  • 4

1 Answers1

8

You assumed that the polygons were in lon/lat coordinates, but they are not:

library(raster)
library(rgdal)
p <- shapefile('Global_Threats/Integrated_Future/rf_int_2030_poly.shp')
p

#class       : SpatialPolygonsDataFrame 
#features    : 63628 
#extent      : -18663508, 14601492, -3365385, 3410115  (xmin, xmax, ymin, ymax)
#coord. ref. : +proj=cea +lon_0=-160 +lat_ts=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 
#variables   : 3
#names       :    ID, THREAT, THREAT_TXT 
#min values  :     1,      0,   Critical 
#max values  : 63628,   2000,  Very High

You can either change the projection

pgeo <- spTransform(p, CRS('+proj=longlat +datum=WGS84'))

and then do something like:

ext <- floor(extent(pgeo))
rr <- raster(ext, res=0.5)
rr <- rasterize(pgeo, rr, field=1)

Or keep the orginal CRS and do something like:

ext <- extent(p)
r <- raster(ext, res=50000)  
r <- rasterize(p, r, field=1)
plot(r)

Note that you are rasterizing very small polygons to large raster cells. A polygon is considered 'inside' if it covers the center of a cell (i.e. assuming a case where polygons cover multiple cells). So for these data you would need to use a much higher resolution (and then perhaps aggregate the results). Alternatively you could rasterize polygon centroids.

But none of the above is relevant really, as you are doing this all backwards. The polygons are clearly derived from a raster (look how blocky they are) and the raster is available in the dataset you point to!

So instead of rasterizing, do:

x <- raster('Global_Threats/Integrated_Future/rf_int_2030')
x
#class       : RasterLayer 
#dimensions  : 25456, 80150, 2040298400  (nrow, ncol, ncell)
#resolution  : 500, 500  (x, y)
#extent      : -20037508, 20037492, -6363885, 6364115  (xmin, xmax, ymin, ymax)
#coord. ref. : NA 
#data source : C:\temp\Global_Threats\Integrated_Future\rf_int_2030 
#names       : rf_int_2030 
#values      : 0, 2000  (min, max)
#attributes  :
#   ID  COUNT THREAT_TXT
#    0  80971        Low
#  100 343535     Medium
# 1000 322231       High
# 1500 168518  Very High
# 2000  83598   Critical

Here plotting a part of Palawan:

e <- extent(c(-8990636, -8929268, 1182946, 1256938))
plot(x, ext=e)
plot(p, add=TRUE)

If you need a lower resolution see raster::aggregate. For a different coordinate reference system, see raster::projectRaster.

Robert Hijmans
  • 40,301
  • 4
  • 55
  • 63
  • RobertH I didn't get to thank you for your assistance here. My original code was based off of the first example I came across here on stackoverflow: http://stackoverflow.com/questions/14992197/shapefile-to-raster-conversion-in-r The resolution of the data in my shapefile made it impossible to rasterize to the resolution I specified. I really liked your simplification at the end of your post. Cheers. – RMFish Mar 17 '16 at 18:00
  • To automatize things a bit, I created https://github.com/brry/misc/blob/master/shp2raster.R – Berry Boessenkool Nov 25 '16 at 10:59