R矢量地图栅格化（将shapefile转换成raster）
R矢量地图栅格化（将shapefile转换成raster）
R矢量地图栅格化（将shapefile转换成raster）
背景
在处理地图数据时候，经常会碰到shp与raster两种格式。通常r中应用较多的为raster栅格数据。shp文件太大，读取也不方便。逐渐被GeoJSON替代，用sf去处理与读取。 R在读取shp时候，处理，或者画图都会碰到，反应迟钝问题。所以，我们有时候会根据需要，将shp文件转成raster，不仅可视化快，还可方便数据处理与提取。shp文件转成raster主要解决以下问题：
- 根据点经纬度提取shp数值
- 计算到某一位置距离，如河流
- 多个属性的ratser合并输出
image.png
下面就来介绍，如何根据shp文件，转成raster及在转换过程中碰到的一些问题。
案例
利用raster包自带的数据进行演示。读取的是SpatialPolygonsDataFrame，关于如何读取shp文件，可以用rgdal与sf的命令。
关键是 rasterize，rasterize(shape, r, 1)里面有三个主要参数：
- shape是shp文件
- r是要栅格化的范围及像素大小；需要先定义
- 1表示，栅格化后，所有值大小
library(raster)
shape = shapefile(system.file("external/lux.shp", package="raster"))
r = raster(shape, res=0.05)
shape_r = rasterize(shape, r, 1)
plot(shape_r)
plot(shape,add=T)
> shape
class : SpatialPolygonsDataFrame
features : 12
extent : 5.74414, 6.528252, 49.44781, 50.18162 (xmin, xmax, ymin, ymax)
crs : +proj=longlat +datum=WGS84 +no_defs
variables : 5
names : ID_1, NAME_1, ID_2, NAME_2, AREA
min values : 1, Diekirch, 1, Capellen, 76
max values : 3, Luxembourg, 12, Wiltz, 312
> shape_r
class : RasterLayer
dimensions : 15, 16, 240 (nrow, ncol, ncell)
resolution : 0.05, 0.05 (x, y)
extent : 5.74414, 6.54414, 49.43162, 50.18162 (xmin, xmax, ymin, ymax)
crs : +proj=longlat +datum=WGS84 +no_defs
source : memory
names : layer
values : 1, 1 (min, max)
par(mfrow=c(1,2))
# value=1
shape_r = rasterize(shape, r, 1)
plot(shape_r)
plot(shape,add=T)
title(main="value=1")
shape_r
# value=2
shape_r = rasterize(shape, r, 2)
plot(shape_r)
plot(shape,add=T)
title(main="value=2")
shape_r
image.png
变量替换
那如果我们需要根据shp里面的地区数来生成不同的value呢，意思就是，不用地区value不一样，不应该是统一值。
par(mfrow=c(1,2))
# value= ID_2
shape_r = rasterize(shape, r, "ID_2")
plot(shape_r)
plot(shape,add=T)
title(main="value=ID_2")
shape_r
# value= AREA
shape_r = rasterize(shape, r, "AREA")
plot(shape_r)
plot(shape,add=T)
title(main="value=AREA")
shape_r
image.png
NA处理
有时候生成的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.
shape_r = rasterize(shape, r, "ID_2")
par(mfrow=c(1,2))
shape_rc=reclassify(shape_r,cbind(NA,0),right=F)
plot(shape_r)
title(main="value=ID_2")
plot(shape_rc)
title(main="NA==0")
image.png
数值提取
转换成raster最终目的是实现数据的提取。譬如现在有两个点，如何提取对应点上的value。 如果是shp文件，操作比较麻烦点，又是还会提取出NA。转换Raster以后，就更方便了。
image.png
# ponits
par(mfrow=c(1,1))
df=tibble(x=c(6.1,5.9),
y=c(49.7,49.9))
df_sp=df
coordinates(df_sp) <- ~ x + y
proj4string(df_sp) <- CRS("+proj=longlat +datum=WGS84 +no_defs")
#plot
plot(shape)
plot(df_sp,add=T,col="red")
# extract value
extract(shape_r,df_sp)
over(df_sp,shape)
image.png
提高精度
上面的图太模糊了，那我们设置res就好。
par(mfrow=c(1,1))
r = raster(shape, res=0.0005)
shape_r = rasterize(shape, r, "ID_2")
plot(shape_r)
plot(shape,add=T)
title(main="Res=0.0005")结论
res精度提高，运行速度会下降，尤其是遇到很大的shp数据时候。
一般R里面加载shp超过50M，系统就会迟钝。
rasterize里面还可以设置field=1.可以达到同样效果。