9

I have a shapefile downloaded from the worldwildlife.org for the terrestrial ecoregions of the world. The file can be loaded here: http://worldwildlife.org/publications/terrestrial-ecoregions-of-the-world.

It comes as a standard shape file and I would like to do two things with it. First: take the shapefile from my local directory and clip it to an extent of eastern North America (ext= extent (-95, -50, 24, 63)) 

# Read shapefile using package "maptools"
eco_shp <- readShapeLines("F:/01_2013/Ecoregions/Global/wwf_terr_ecos.shp", 
                          proj4string=CRS("+proj=utm +zone=33 +datum=WGS84")) 


# Set the desired extent for the final raster using package "raster" 
ext <- extent(-95, -50, 24, 63)

I am sure I have to use the rasterize function in the package "raster" but I am still not able to get it work correctly. I would appreciate any suggestions on how to do this.

Andre Silva
  • 4,782
  • 9
  • 52
  • 65
I Del Toro
  • 913
  • 4
  • 15
  • 36
  • Do you need it to be rasterised? Would cliping the shapefile by a polygon suffice? – mnel Feb 21 '13 at 00:03
  • Mnel: I do need it to be rasterized into an .asc file to match up with my other environmental layers. – I Del Toro Feb 21 '13 at 00:06
  • The first step of clipping by a shapfile would be helpful, then that new object could be rasterized and might be less intensive . – I Del Toro Feb 21 '13 at 00:08
  • why not just overlay the polygon onto the SpatialGridDataFrame you wish to match it to, (use `overlay`) – mnel Feb 21 '13 at 00:12
  • Mnel: My model requires that my input dataframes all be in an asc format and have the same spatial extent. If I only wanted to visualize it the "overlay" would work. :) – I Del Toro Feb 21 '13 at 00:16
  • The `overlay` function in R will extract the values at the grid points, you could then save as an ascii file. – mnel Feb 21 '13 at 00:18
  • Im sorry, I am new to working on GIS in R. Is the idea to make a empty raster with the desired extent then overlay the shapefile and extract the values to the points that way? – I Del Toro Feb 21 '13 at 00:26
  • Yes, that's the idea. Have a look at the answer I just posted for an example. – Josh O'Brien Feb 21 '13 at 00:28

1 Answers1

14

You are right to think that you should be using raster (rather than the sp raster Spatial classes) for spatial raster data. You should also use rgdal (rather than maptools) for reading, writing, and otherwise manipulating spatial vector data.

This should get you started:

library(rgdal)
library(raster)

## Read in the ecoregion shapefile (located in R's current working directory)
teow <- readOGR(dsn = "official_teow/official", layer = "wwf_terr_ecos")

## Set up a raster "template" to use in rasterize()
ext <-  extent (-95, -50, 24, 63)
xy <- abs(apply(as.matrix(bbox(ext)), 1, diff))
n <- 5
r <- raster(ext, ncol=xy[1]*n, nrow=xy[2]*n)

## Rasterize the shapefile
rr <-rasterize(teow, r)

## A couple of outputs
writeRaster(rr, "teow.asc")
plot(rr)

enter image description here

Bhargav Rao
  • 50,140
  • 28
  • 121
  • 140
Josh O'Brien
  • 159,210
  • 26
  • 366
  • 455
  • 1
    Glad to hear that. Note that in practice, you'll probably want to use one of your own raster layers as the template raster (`r` in `rasterize(teow, r)`), and may need to do a bit of fiddling to get the proj4strings matched up (although both **raster** and **rgdal** are really good about handling projection metadata.). – Josh O'Brien Feb 21 '13 at 00:48
  • 1
    It should be noted that rasterize() appears to be impossibly unstable and inefficient with large datasets. I recently tried to rasterize a coverage of the city of Seattle at 6ft resolution (the output is a 14 MB GeoTIFF file, 8095 x 14819) and R spent about three hours running 7 threads before "finishing" with no output and no error message either. Using R to generate a blank GeoTIFF raster of the desired extent and resolution, and then running the rasterize operation via gdal (with some help from QGIS) took less than half an hour in a single thread. – Miles Erickson Dec 24 '15 at 03:30
  • 3
    @MilesErickson Indeed. For any but the smallest rasters, I tend to use a wrapper around `gdalUtils::gdal_rasterize`. For 1000-by-1000 rasters, I get a 6-fold speed-up, and for 2000-by-2000 rasters, I get a 12-fold speed-up (from 70 to 6 seconds). – Josh O'Brien Dec 25 '15 at 06:52
  • Great tip! Wow, so rasterize() might have taken weeks or years to finish my 8000x15000 raster, had it been able to at all... – Miles Erickson Dec 26 '15 at 07:25
  • Hey @Josh you created variable `n` but did not use it? Did you mean `r <- raster(ext, ncol=xy[1]*n, nrow=xy[2]*n)`? – Bhargav Rao Jan 20 '16 at 18:45
  • @BhargavRao I did indeed. Thanks! (Fixed it.) – Josh O'Brien Jan 20 '16 at 19:10
  • Nice, @Josh Thanks. I have made a small edit to reflect the code formatting. – Bhargav Rao Jan 20 '16 at 19:22
  • For a raster 79862 by 67313, gdal(python) takes about 100 seconds. The R raster package takes all day and usually crashes. The difference is that gdal is a C application and is optimized...whereas the code for raster->rasterize is written in R and uses way too much memory. – user1269942 Nov 10 '17 at 01:49