[R-sig-Geo] Crop a raster using a shapefile

[Thiago Veloso]

  Thanks a lot to all who helped me.

  This combination of crop and mask was exactly what I needed: crop a raster using not only the extents of a shapefile, but also its contours, its curves.

  Cheers,
  Thiago.

________________________________
From: Mauricio Zambrano-Bigiarini <mauricio.zambrano at jrc.ec.europa.eu>
To: Thiago Veloso <thi_veloso at yahoo.com.br> 
Cc: R-SIG list <r-sig-geo at r-project.org> 
Sent: Monday, June 4, 2012 5:46 AM
Subject: Re: [R-sig-Geo] Crop a raster using a shapefile

On 01/06/12 20:49, Thiago Veloso wrote:
>    Dear all,
>
>    When cropping a raster using a shapefile, through 'crop' function from "raster" package, only the extent is taken into account:
>
>> r<-raster("lai2011361.Lai_1km.tif")
>> extent(r)
>
> class       : Extent
> xmin        : -104.4326
> xmax        : -29.99736
> ymin        : -40.00064
> ymax        : 10
>
>> br<-readShapePoly("/mnt/disco3/MODIS/shapes/brazil.shp")
>
>> extent(br)
>
> class       : Extent
> xmin        : -73.83943
> xmax        : -34.8581
> ymin        : -33.77086
> ymax        : 5.38289
>
>> c<-crop(r,br)
>
>> extent(c)
> class       : Extent
> xmin        : -73.83772
> xmax        : -34.85554
> ymin        : -33.76852
> ymax        : 5.38428
>
>    However, the cropped raster is rectangular. Is there any way to keep also the "shape" of a shapefile (and not only its extent) when performing this kind of crop operation?
>

You may use a combination of crop and mask:

------------ START ---------------
library(raster)

# Reading the shapefile
myshp <- readOGR("myshapefile.shp", layer=basename(myshapefile) )

# Getting the spatial extent of the shapefile
e <- extent(myshp)

# Reading the raster you want to crop
myraster <- raster(myraster.filename)

# Cropping the raster to the shapefile spatial extent
myraster.crop <- crop(raster, e, snap="out")

# Dummy raster with a spatial extension equal to the cropped raster,
# but full of NA values
crop          <- setValues(myraster.crop, NA)

#  Rasterize the catchment boundaries, with NA outside the catchment 
boundaries
myshp.r <- rasterize(myshp, crop)

# Putting NA values in all the raster cells outside the shapefile boundaries
myraster.masked <- mask(x=myraster.crop, mask=myshp.r)

------------ END ---------------

IHTH,

Mauricio

-- 
====================================================
Water Resources Unit
Institute for Environment and Sustainability (IES)
Joint Research Centre (JRC), European Commission
webinfo    : http://floods.jrc.ec.europa.eu/
====================================================
DISCLAIMER:
"The views expressed are purely those of the writer
and may not in any circumstances be regarded as sta-
ting an official position of the European Commission"
====================================================
Linux user #454569 -- Ubuntu user #17469
====================================================
"The only things in life you regret,
Are the risks that you didn't take"
(Anonymous)
>    Thanks in advance,
>    Thiago.
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>