2013-10-11 45 views
2

我想用两个不同范围和不同投影系统的geotiff文件生成2个子图的图。我想根据Rapideye范围来裁剪这个地块。 我该怎么办? 以下是该文件的详细信息。如何用不同的投影来裁剪光栅

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 

回答

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,我首先从光栅R2的程度创建spatialPolygon的程度能够作物光栅R1,第二分配给它的良好的投影和第三变换所述多边形的投影光栅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))) 

最后,可以使用多边形r2extr1proj表示R2的r1中的投影范围裁剪光栅R1。然后绘制两个栅格。

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

+1 - 非常棒的回答!谢谢。 –