R矢量地图栅格化（将shapefile转换成raster）

2020-10-23 本文已影响0人 jamesjin63

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.可以达到同样效果。

R矢量地图栅格化（将shapefile转换成raster）