我想用两个具有不同范围和不同投影系统的geotiff文件中的2个子图制作一个图。我想根据Rapideye范围裁剪图。我应该怎么做 ？以下是文件的详细信息。
点到点
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
这可能不是最优雅的方法，但可能会有所帮助。我根据您的示例松散地创建了两个示例栅格。它们具有相同的投影和范围。
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的范围创建一个spacealPolygon，然后为其分配良好的投影，然后将多边形转换为栅格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)
本文收集自互联网，转载请注明来源。
如有侵权，请联系[email protected] 删除。
我来说两句