问题描述
我想将栅格数据聚合到自定义 shapefile 中的每个多边形。
在这种情况下,我想获得撒哈拉以南非洲次国家区域内的平均城市化程度。
我的科幻是这样的:
> africa_map
Simple feature collection with 543 features and 4 fields
Geometry type: MULTIpolyGON
Dimension: XY
Bounding Box: xmin: -25.36042 ymin: -46.96575 xmax: 63.49391 ymax: 27.66147
Geodetic CRS: WGS 84
First 10 features:
cname ccode regname continent geometry
1 Angola AO Bengo Africa MULTIpolyGON (((13.371 -8.5...
2 Angola AO Benguela Africa MULTIpolyGON (((12.53336 -1...
3 Angola AO Bie Africa MULTIpolyGON (((16.61158 -1...
4 Angola AO Cabinda Africa MULTIpolyGON (((12.78266 -4...
5 Angola AO cuando Cubango Africa MULTIpolyGON (((21.9838 -16...
6 Angola AO cuanza norte Africa MULTIpolyGON (((15.40788 -7...
7 Angola AO cuanza Sul Africa MULTIpolyGON (((13.7926 -11...
或绘制:
另一方面,栅格数据采用这种形式:
> imported_raster
class : RasterLayer
dimensions : 18000,36082,649476000 (nrow,ncol,ncell)
resolution : 1000,1000 (x,y)
extent : -18041000,18041000,-9e+06,9e+06 (xmin,xmax,ymin,ymax)
crs : +proj=moll +lon_0=0 +x_0=0 +y_0=0 +datum=wgs84 +units=m +no_defs
names : GHS_BUILT_LDS1975_GLOBE_R2018A_54009_1K_V2_0
values : 0,100 (min,max)
对于整个地球而言,这些都比需要的要细得多。为了加速计算,我首先聚合了栅格,然后将其转换为 shapefile,每个剩余的栅格像素都转换为 shapefile 中的点几何。然后,这个 shapefile 可以聚合到我的区域边界。诚然,这不是很优雅(并且最终不起作用)。
library(tidyverse)
library(sf)
library(raster)
library(stars)
library(rgdal)
> # aggregate (to 25x25 km)
> imported_raster_ag <- aggregate(imported_raster,fact=25)
>
> # convert to sp
> urbanized = rasterTopolygons(imported_raster_ag)
> # convert to sf
> urbanized_sf <- st_as_sf(urbanized)
# compare projection
st_crs(africa_map)==st_crs(urbanized_sf)
# align projection
urbanized_sf <- st_transform(urbanized_sf,st_crs(africa_map))
urbanized_sf <- urbanized_sf %>% rename(urbanization = GHS_BUILT_LDS1975_GLOBE_R2018A_54009_1K_V2_0)
> urbanized_sf
Simple feature collection with 398872 features and 1 field (with 500 geometries empty)
Geometry type: polyGON
Dimension: XY
Bounding Box: xmin: -179.9999 ymin: -57.40086 xmax: 179.9963 ymax: 82.33738
Geodetic CRS: WGS 84
First 10 features:
urbanization geometry
1 0 polyGON ((-117.1367 82.3373...
2 0 polyGON ((-116.2261 82.3373...
3 0 polyGON ((-115.3156 82.3373...
4 0 polyGON ((-114.405 82.33738...
5 0 polyGON ((-113.4944 82.3373...
6 0 polyGON ((-112.5838 82.3373...
7 0 polyGON ((-111.6732 82.3373...
8 0 polyGON ((-110.7627 82.3373...
9 0 polyGON ((-109.8521 82.3373...
10 0 polyGON ((-108.9415 82.3373...
当时的想法是,我可以沿着其他科幻小说的区域边界聚合这些点。 但是,我收到一条错误消息,我发现的唯一推荐修复只是重复了它。
> urbanized_africa <- aggregate(urbanized_sf["urbanization"],by = africa_map$geometry,mean)
although coordinates are longitude/latitude,st_intersects assumes that they are planar
Error in CPL_geos_binop(st_geometry(x),st_geometry(y),op,par,pattern,:
Evaluation error: IllegalArgumentException: Invalid number of points in LinearRing found 2 - must be 0 or >= 4.
> urbanized_sf_fixed <- sf::st_make_valid(urbanized_sf)
Error in CPL_geos_make_valid(x) :
Evaluation error: IllegalArgumentException: Invalid number of points in LinearRing found 2 - must be 0 or >= 4.
我想在我的转换过程中的某个地方可能会损坏或其他一些更根本的缺陷。将栅格数据聚合为 shapefile 多边形的更优雅、最重要的是功能更强大的工作流是什么?
解决方法
查看 extract
函数
在你的情况下像
africamap$urbanized <- extract(imported_raster,africamap,fun="mean")
africamap
仍应首先投影到光栅投影。如果您有 nodata
值,那么向调用添加 na.rm=TRUE
参数应该是明智之举。
为了缩短处理时间,您还可以将栅格裁剪为非洲。
extract
包中还有一个更快的 exactextractr
版本。您还可以利用 terra
而不是 raster
进行光栅处理。