裁剪R中使用shape文件栅格,但保持shape文件的几何

Question

我使用{光栅}夹(或作物)的基础上的不规则shape文件(亚马逊生物群落),但输出总是具有矩形程度的栅格.不过,我需要在shape文件完全相同的几何形状的输出.有小费吗？干杯.

library(raster)
library(rgdal)

myshp <- readOGR("Amazon.shp", layer="Amazon")
e <- extent(myshp)
myraster <- raster("Temperature.tif")
myraster.crop <- crop(myraster, e, snap="out", filename="myoutput.tif")

Answer 1

一种选择是使用 raster::mask()

library(maptools)  ## For wrld_simpl
library(raster)

## Example SpatialPolygonsDataFrame
data(wrld_simpl)
SPDF <- subset(wrld_simpl, NAME=="Brazil")

## Example RasterLayer
r <- raster(nrow=1e3, ncol=1e3, crs=proj4string(SPDF))
r[] <- 1:length(r)

## crop and mask
r2 <- crop(r, extent(SPDF))
r3 <- mask(r2, SPDF)

## Check that it worked
plot(r3)
plot(SPDF, add=TRUE, lwd=2)

Answer 2

terra::包使其变得更简单，您可以在同一步骤中::裁剪和遮罩。

# Load packages
library(maptools) # For geometry
library(terra) # Perform the crop and mask

###--- Preparing polygon and raster ---###

# Example SpatialPolygonsDataFrame
data(wrld_simpl)
polygon <- subset(wrld_simpl, NAME=="Luxembourg")
plot(polygon) # have a look

# Convert from SpatialPolygonsDataFrame to SpatVector (terra package format)
# And create a smaller polygon with buffer (negative to be "inside")
polygon_bf <- buffer(vect(polygon), width= -100000)
plot(polygon_bf, add= T) # have a look on both

# Create a SpatRaster from a file
f <- system.file("ex/elev.tif", package="terra")
r_lux <- rast(f)
plot(r_lux) # have a look on SpatRaster (terra package format)

# See the steps with plot
plot(polygon, add= T)
plot(polygon_bf, add= T)

点击查看中间步骤（所有文件）

########################################################
### Crop and mask by any polygon                     ###
raster_cp <- crop(r_lux, polygon_bf, mask= T)
# Note: if mask= F, the crop will be by extent (box) ###
########################################################

### Check the results
plot(raster_cp)
plot(polygon_bf, lwd=1, add=T)

点击查看最终输出

此外，（使用光栅包）如果您想使用简单的几何图形（ig box）执行，则可以直接放置范围的坐标：

e <- as(extent(c(xmin= -16, xmax= -7.25, ymin= 4, ymax= 12.75)), 'SpatialPolygons')
crs(e) <- "+proj=longlat +datum=WGS84 +no_defs"
r <- crop(my_raster, e)