CODEWITHSUNDEEP
rplotraster
Wraf
Wraf

Reputation: 757

How to crop a raster with a different projection

I would like to make a figure with 2 subplots from two geotiff files with different extent and different projection system. I would like to crop the plot based on the Rapideye extent. How should I do ? Below are the details of the file.

SPOT-VGT

class       : RasterLayer 
dimensions  : 8961, 8961, 80299521  (nrow, ncol, ncell)
resolution  : 0.008928571, 0.008928571  (x, y)
extent      : -20, 60.00893, -40.00893, 40  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
data source : /AFRI_VGT_V1.3.tiff 
names       : g2_BIOPAR_WB.GWWR_201305110000_AFRI_VGT_V1.3 
values      : 0, 255  (min, max)

RAPIDEYE

class       : RasterStack 
dimensions  : 14600, 14600, 213160000, 5  (nrow, ncol, ncell, nlayers)
resolution  : 5, 5  (x, y)
extent      : 355500, 428500, 2879500, 2952500  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=utm +zone=36 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 
names       : /rapideye.tif 
min values  :      0,         0,         0,       0,       0 
max values  :           65535,      65535,   65535,         65535,     65535

Upvotes: 2

Views: 1096

Answers (1)

Marie Auger-Methe
Marie Auger-Methe

Reputation: 838

This might not be the most elegant way but might help. I have created two example rasters loosely based on your examples. They have the same projection and extent.

library(raster)
r1 <- raster(nrows=500, ncols=500, 
         ext=extent(c(-20, 60.00893, -40.00893, 40)),
         crs='+proj=longlat +datum=WGS84')
r1[] <- rnorm(500*500,0,1)

r2 <- raster(nrows=50, ncols=50, 
         ext=extent(c(355500, 428500, 2879500, 2952500)),
         crs='+proj=utm +zone=36 +datum=WGS84 +units=m')
r2[] <- rnorm(50*50,0,1)

To be able to crop raster r1 using the extent of raster r2, I am first creating a spatialPolygon from the extent of raster r2, second assigning it the good projection, and third transforming the polygon to the projection of raster r1.

library(rgdal) 
# Make a SpatialPolygon from the extent of r2
r2extent <- as(extent(r2), 'SpatialPolygons')
# Assign this SpatialPolygon the good projection
proj4string(r2extent) <- proj4string(r2)
# Transform the projection to that of r1
r2extr1proj <- spTransform(r2extent, CRS(proj4string(r1)))

Finally, you can crop raster r1 using the polygon r2extr1proj which represents the extent of r2 in the projection of r1. Then plot the two rasters.

r1crop <- crop(r1, r2extr1proj)
layout(matrix(c(1:2), nrow=1))
plot(r1crop)
plot(r2)

Upvotes: 5

Related Questions