R矢量地图栅格化(将shapefile转换成raster)

R矢量地图栅格化(将shapefile转换成raster)

背景

在处理地图数据时候,经常会碰到shpraster两种格式。通常r中应用较多的为raster栅格数据。shp文件太大,读取也不方便。逐渐被GeoJSON替代,用sf去处理与读取。
R在读取shp时候,处理,或者画图都会碰到,反应迟钝问题。所以,我们有时候会根据需要,将shp文件转成raster,不仅可视化快,还可方便数据处理与提取。shp文件转成raster主要解决以下问题:

  1. 根据点经纬度提取shp数值
  2. 计算到某一位置距离,如河流
  3. 多个属性的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.可以达到同样效果。

参考

  1. 栅格化shp数据
  2. Rasterize polygons with R
  3. 替换raster中NA数据
  4. 根据shp裁剪raster地图
  5. [sf裁剪 https://rpubs.com/cyclemumner/423273)
  6. problem cropping sf object to extent of another sf polygon in R
最后编辑于
©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 151,829评论 1 331
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 64,603评论 1 273
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 101,846评论 0 226
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 42,600评论 0 191
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 50,780评论 3 272
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 39,695评论 1 192
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 31,136评论 2 293
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 29,862评论 0 182
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 33,453评论 0 229
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 29,942评论 2 233
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 31,347评论 1 242
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 27,790评论 2 236
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 32,293评论 3 221
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 25,839评论 0 8
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 26,448评论 0 181
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 34,564评论 2 249
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 34,623评论 2 249