当您使用该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 列的值。