如何使用不同的投影裁剪栅格答案

【问题标题】：How to crop a raster with a different projection如何使用不同的投影裁剪栅格
【发布时间】：2013-10-19 17:40:30
【问题描述】：

我想从两个不同范围和不同投影系统的 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)

快速眼

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

【问题讨论】：

标签： r plot raster

【解决方案1】：

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

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)))

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

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

【讨论】：