R矢量地图栅格化（将shapefile转换成raster）
Posted
tags:
篇首语：本文由小常识网(cha138.com)小编为大家整理，主要介绍了R矢量地图栅格化（将shapefile转换成raster）相关的知识，希望对你有一定的参考价值。
参考技术A 在处理地图数据时候，经常会碰到 shp 与 raster 两种格式。通常r中应用较多的为raster栅格数据。shp文件太大，读取也不方便。逐渐被 GeoJSON 替代，用sf去处理与读取。R在读取shp时候，处理，或者画图都会碰到，反应迟钝问题。所以，我们有时候会根据需要，将shp文件转成raster，不仅可视化快，还可方便数据处理与提取。shp文件转成raster主要解决以下问题：
下面就来介绍，如何根据shp文件，转成raster及在转换过程中碰到的一些问题。
利用 raster 包自带的数据进行演示。读取的是 SpatialPolygonsDataFrame ，关于如何读取shp文件，可以用rgdal与sf的命令。
关键是 rasterize， rasterize(shape, r, 1) 里面有三个主要参数：
那如果我们需要根据shp里面的地区数来生成不同的value呢，意思就是，不用地区value不一样，不应该是统一值。
有时候生成的raster里面有NA数据，那么如何替换掉呢，(reclassify)[ http://search.r-project.org/library/raster/html/reclassify.html] 可以实现该过程。主要参数cbind（0,a,b）意思是将0-a的数值全部变成b。
具体参见： ?reclassify
下面我么将NA替换成0，或者，value=12的替换成100.
转换成raster最终目的是实现数据的提取。譬如现在有两个点，如何提取对应点上的value。
如果是shp文件，操作比较麻烦点，又是还会提取出NA。转换Raster以后，就更方便了。
上面的图太模糊了，那我们设置res就好。
res精度提高，运行速度会下降，尤其是遇到很大的shp数据时候。
一般R里面加载shp超过50M，系统就会迟钝。
rasterize 里面还可以设置field=1.可以达到同样效果。
R：将栅格聚合为 shapefile 多边形
【中文标题】R：将栅格聚合为 shapefile 多边形【英文标题】：R: Aggregating raster to shapefile polygons 【发布时间】：2021-08-15 11:05:24 【问题描述】：我想将栅格数据聚合到自定义 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...
我的想法是，到目前为止，我可以将这些点沿其他 SF 的区域边界聚合起来。 但是，我收到一条错误消息，并且我发现的唯一推荐修复只是重复它。
> 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 多边形中有哪些更优雅、更重要的是功能更强大的工作流？
【问题讨论】：
【参考方案1】：看看extract函数
在你的情况下
africamap$urbanized <- extract(imported_raster, africamap, fun="mean")
africamap 仍应首先投影到光栅投影。如果您有nodata 值，那么在调用中添加na.rm=TRUE 参数应该是明智的。
为了缩短处理时间，您还可以将栅格裁剪到非洲。
exactextractr 包中还有更快的extract 版本。您还可以利用terra 而不是raster 进行光栅处理。
【讨论】：
很好的建议，非常感谢！使用推荐的 extractr 包，聚合工作非常迅速，不需要裁剪，我得到了想要的结果。对齐投影：africa_map_mollweide = st_transform(africa_map, crs = "+proj=moll") 并提取栅格数据：africa_map_mollweide$urbanized <- exact_extract(imported_raster_ag, africa_map_mollweide, fun="weighted_mean", weights="area")以上是关于R矢量地图栅格化（将shapefile转换成raster）的主要内容，如果未能解决你的问题，请参考以下文章