将多个多边形的shapefile转换为R中的栅格

Question

我在将多边形转换为 R 中的栅格时遇到了很大的麻烦。我想做的是：我有 574 种形状文件（即多边形）。即在属性 table 中它有 574 行（即 FID 在 0 到 573 之间）。可以在此处找到数据的子集：https://drive.google.com/file/d/1AdTChjerCXopE1-PZIAPp5ZXISdq45i8/view?usp=sharing

我想将其转换为光栅。在输出栅格中，我看到最小值和最大值分别为 1 和 574。我怀疑的是：它正在获取单元格中的字段 ID 作为不应该的像素值。单元格值应来自覆盖多边形。任何帮助将不胜感激。下面是示例代码：

library(raster)
library(rgdal)
library(maptools)


# define porjection
projection1 <- CRS ("+proj=longlat +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +no_defs")

# load species polygons for all 574 species
# already in the home directory

sp <- readShapePoly("AllP574sp.shp", proj4string = projection1)

# load a raster file to use as a mirror for rasterize
raster1 <- raster("/data/projects/MeanTemp2050rcp45_BCC_CSM1_1.tif")


r.sp <- rasterize(sp, raster1) # rasterize our species polygon to the same resoluton of loaded raster

writeRaster(r.sp, "/data/projects/all574spRaster", format = "GTiff", overwrite = TRUE)

输出栅格的属性如下：

class       : RasterLayer 
dimensions  : 18000, 43200, 777600000  (nrow, ncol, ncell)
resolution  : 0.008333333, 0.008333333  (x, y)
extent      : -180, 180, -60, 90  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +no_defs 
data source : E:\all574spRaster.tif 
names       : all574spRaster 
values      : 1, 574  (min, max)

Answer 1

当您使用 rasterize 函数时，指定 field 参数很重要，否则默认情况下它会尝试为您创建一个；在您的情况下，它看起来像是从 FID 列创建了一个。

我做了一些猜测以重新生成一组可能与您的相似的多边形工作集。

library(maptools)
library(rgdal)
library(sp)
library(geosphere)

# set seed for duplicatable results
set.seed(1)

# some data that looks a little like yours
BINOMIAL <- c("Controversial chimneyswift", "Dull dungbeetle", 
"Easternmost eel", "Jumping jaeger", "Qualified queenconch")
FID <- 0:(length(BINOMIAL) - 1)
RANGE <-  runif(length(BINOMIAL), min = 118, max = 3875370)
MyData <- cbind.data.frame(FID, BINOMIAL, RANGE)
row.names(MyData) <- FID

# some semi-random polygons in your extent box
ext <- extent(c(-180, 180, -60, 90))
create_polygon <- function(n = 4, lat, lon, r) {
  lengths <- rnorm(n, r, r/3)
  smoother_lengths <- c(sort(lengths), rev(sort(lengths)))
  lengths <- smoother_lengths[sort(sample(n * 2, n))]
  lengths <- rep(lengths[1], length(lengths))
  directions <- sort(runif(n, 0, 360))
  p <- cbind(lon, lat) 
  vertices <- t(mapply(destPoint, b = directions, 
                       d = lengths, MoreArgs = list(p = p)))
  vertices <- rbind(vertices, vertices[1, ])
  sapply(vertices[,1], min, ext@xmax)
  sapply(vertices[,1], max, ext@xmin)
  sapply(vertices[,2], min, ext@ymax)
  sapply(vertices[,2], max, ext@ymin)
  Polygon(vertices) 
}
rand_lats <- runif(nrow(MyData), min = -50, max = 60)
rand_lons <- runif(nrow(MyData), min = -100, max = 100)
rand_sides <- sample(4:20, nrow(MyData), replace = TRUE)
rand_sizes <- rnorm(nrow(MyData), mean = 5e+06, sd = 1e+06)
make_species_polygon <- function(i) {
  p.i <- list(create_polygon(rand_sides[i], rand_lats[i], 
                             rand_lons[i], rand_sizes[i]))
  P.i <- Polygons(p.i, FID[i])
}

polys <- SpatialPolygons(lapply(1:nrow(MyData), make_species_polygon))
spdf <- SpatialPolygonsDataFrame(Sr = polys, data = MyData)

t.shp <- tempfile(pattern = "MyShapefile", fileext = ".shp")
raster::shapefile(spdf, t.shp)

此时在您的临时目录中写入了一个shapefile，其名称存储在变量t.shp中。我希望该 shapefile 成为您真实的任何大 shapefile 的可行副本。所以现在我们可以查看您的代码，它在做什么，以及您希望它做什么：

## now we get into your code
library(raster)
library(rgdal)
library(maptools)

# define porjection
projection1 <- CRS ("+proj=longlat +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +no_defs")

## 
## I don't know what your shapefile looks like exactly, 
## but substituting `t.shp` the tempfile that I created above
## also since the function readShapePoly is deprecated
## instead I use the recommended new function rgdal::readOGR() 
## 
sp <- rgdal::readOGR(t.shp)

##
## I don't know what your tiff file looks like exactly,
## but I can duplicate its characteristics
## for speed I have decreased resolution by a factor of 10
##
raster1 <- raster(nrow = 1800, ncol = 4320, ext)

# rasterize our species polygon to the same resoluton of loaded raster
r.sp <- rasterize(x = sp, y = raster1, field = MyData$RANGE) 

t.tif <- tempfile(pattern = "MyRastfile", fileext = ".tif")
writeRaster(r.sp, t.tif, format = "GTiff", overwrite = TRUE)

结果如下：

raster(t.tif)
class       : RasterLayer 
dimensions  : 1800, 4320, 7776000  (nrow, ncol, ncell)
resolution  : 0.08333333, 0.08333333  (x, y)
extent      : -180, 180, -60, 90  (xmin, xmax, ymin, ymax)
coord. ref. : NA 
data source : [["a file name in your temp directory"]] 
names       : MyRastfile1034368f3cec 
values      : 781686.3, 3519652  (min, max)

结果现在显示取自 RANGE 列而不是 FID 列的值。

将多个多边形的shapefile转换为R中的栅格

convert shapefile of multiple polygons to raster in R

r

polygon

arcgis

shapefile

rasterize