2

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
4

1 回答 1

3

这可能不是最优雅的方式，但可能会有所帮助。我根据您的示例松散地创建了两个示例栅格。它们具有相同的投影和范围。

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)

为了能够使用栅格 r2 的范围裁剪栅格 r1，我首先从栅格 r2 的范围创建一个 spatialPolygon，然后为其分配良好的投影，然后将多边形转换为栅格 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)))

最后，您可以使用表示 r2 在 r1 投影中的范围的多边形 r2extr1proj 裁剪栅格 r1。然后绘制两个栅格。

r1crop <- crop(r1, r2extr1proj)
layout(matrix(c(1:2), nrow=1))
plot(r1crop)
plot(r2)
于 2013-10-11T16:22:02.013 回答