如何解决从栅格文件中提取数据以与shp文件匹配
更新: 更新:问题是我生产网格单元的方式导致我的shapefile投影错误,从而导致上述错误。我能够产生另一种形式的网格,这些网格与我的投影效果完美匹配。
我成功创建了世界地图(How to generate 10x10km grid cells of all countries?)的栅格网格
library(sf)
library(stars)
library(rnaturalearth)
# Polygon
world = ne_countries(scale = "small",returnclass = "sf")
world = st_transform(world,"+proj=robin +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs")
pol = world
# Make grid
grid = st_as_stars(st_bbox(pol),dx = 10000,dy = 10000)
grid = st_as_sf(grid)
grid = grid[pol,]
grid = st_transform(grid,"+proj=robin +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs")
st_write(grid,"data/grid/grid.shp",driver = "ESRI Shapefile")
我现在正在尝试使用hyde数据(https://dataportaal.pbl.nl/downloads/HYDE/HYDE3.2/baseline.zip)将其与每个单元的种群信息合并。
我下载了Bayseline Hyde数据,并试图通过如下所示的提取函数将生成的网格应用于总体数据:
## use extract function for nightlights
shape <- rgdal::readOGR("data/grid/grid.shp")
source("scripts/extractfunction_hyde.R")
extract_hyde(directory = "data/hyde/hyde_harmonized",shape)
extract_hyde <- function(directory = ".",shp,years = NULL) {
require(raster)
require(velox)
#sname <- shpname
#shpname <- paste0("data/FDI/",shpname)
#shp <- readRDS(shpname)
if (!class(shp) %in% c("SpatialPolygons","SpatialPolygonsDataFrame","SpatialPointsDataFrame")) {
stop(paste("'shp' must be either a SpatialPolygons","SpatialPolygonsDataFrame or SpatialPointsDataFrame"))
}
crs <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0"
shp <- sp::spTransform(shp,sp::CRS(crs))
orig.dir <- getwd()
setwd(directory)
files <- list.files(pattern = "*.asc$")
# Years in which this population data are available:
all.years <- as.numeric(substr(files,6,9)) # The year is characters 4 to 9
# Need to average the years where there are two satellite readings:
double.years <- all.years[duplicated(all.years)]
# If years aren't provided,take all of them:
if (is.null(years)) {
years <- sort(unique(all.years))
}
# Start the output data.frame:
if (class(shp) == "SpatialPolygons") {
df <- data.frame(id = 1:length(shp@polygons))
} else if (class(shp) == "SpatialPolygonsDataFrame") {
df <- data.frame(shp@data)
} else if (class(shp) == "SpatialPointsDataFrame") {
df <- data.frame(shp@data)
}
for (i in seq_along(years)) {
cat("Extracting data for year ",years[i],"...",sep = "")
# If there are two satellite readings in a year,average them first:
if (years[i] %in% double.years) {
both.files <- grep(years[i],files,value = TRUE)
r <- crop(raster(both.files[1]),snap = "out")
r2 <- crop(raster(both.files[2]),snap = "out")
values(r) <- (values(r) + values(r2)) / 2
rm(r2)
r <- velox(r)
# With only one reading in a year,just read in the file normally:
} else {
r <- crop(raster(grep(years[i],value = TRUE)),snap = "out")
r <- velox(r)
}
extract <- r$extract(sp=shp,fun = function(x) mean(x,na.rm = TRUE))
df[[paste0("pop_mean",years[i])]] <- c(extract)
extract <- r$extract(sp=shp,fun = function(x) sum(x,na.rm = TRUE))
df[[paste0("pop_sum",years[i])]] <- c(extract)
cat("Done\n")
}
saveRDS(df,"data/hyde/hyde_grid.Rds")
}
不幸的是,我收到以下错误代码:
[,1] [,2] [,3] [,4]
[1,] -9070131 -8626996 Inf Inf
[2,] -9080131 -8626996 Inf Inf
Error in .spTransform_Polygon(input[[i]],to_args = to_args,from_args = from_args,:failure in Polygons 1584553 Polygon 1 points 3:4
该功能适用于我使用过的所有其他形状文件。例如,我遵循了另一个stackoverflow帖子中的第二个选项,并且效果很好。
我进行了研究,当有不同的预测时,有时会出现此问题。但是我的函数可以调整hyde文件和shape文件的投影。
我认为我收到此错误是因为未正确生成shp文件?
版权声明:本文内容由互联网用户自发贡献,该文观点与技术仅代表作者本人。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如发现本站有涉嫌侵权/违法违规的内容, 请发送邮件至 dio@foxmail.com 举报,一经查实,本站将立刻删除。