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.

8
  • Do you need it to be rasterised? Would cliping the shapefile by a polygon suffice?
    – mnel
    Feb 21, 2013 at 0: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, 2013 at 0: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, 2013 at 0:08
  • why not just overlay the polygon onto the SpatialGridDataFrame you wish to match it to, (use overlay)
    – mnel
    Feb 21, 2013 at 0: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, 2013 at 0:16

1 Answer 1

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

8
  • 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.). Feb 21, 2013 at 0: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. Dec 24, 2015 at 3: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). Dec 25, 2015 at 6:52
  • Great tip! Wow, so rasterize() might have taken weeks or years to finish my 8000x15000 raster, had it been able to at all... Dec 26, 2015 at 7: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)? Jan 20, 2016 at 18:45

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.