R ggplot2 与 shapefile 和 csv 数据合并以填充多边形 | 珊瑚贝

R ggplot2 merge with shapefile and csv data to fill polygons


我们每天制作地图,显示我们地区 30 个不同区域的计算温度水平,每个区域根据水平用不同的颜色填充。这张地图看起来像

enter

现在我想将地图生成切换到 R。我已经下载了省和市边界(您可以找到整个西班牙的边界或我所在地区的子集)并设法按照 Hadley 的示例使用 ggplot2 绘制它们。

我还可以生成一个包含两列的 ascii 文件:标识符 (CODINE) 和每日级别。你可以在这里下载。

这是我第一个尝试使用 R 和 ggplot2 绘制 shapefile 的脚本,因此可能会出现错误,并且可以肯定它可以改进,欢迎提出建议。以下代码(基于 Hadley 前面提到的)对我有用:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
> require(“rgdal”)
> require(“maptools”)
> require(“ggplot2”)
> require(“plyr”)

# Reading municipal boundaries

esp = readOGR(dsn=”.”, layer=”lineas_limite_municipales_etrs89″)

muni=subset(esp, esp$PROV1 ==”46″ | esp$PROV1 ==”12″ | esp$PROV1 ==”3″)
muni@data$id = rownames(muni@data)
muni.points = fortify(muni, region=”id”)
muni.df = join(muni.points, muni@data, by=”id”)

# Reading province boundaries

prov = readOGR(dsn=”.”, layer=”poligonos_provincia_etrs89″)

pr=subset(prov, prov$CODINE ==”46″ | prov$CODINE ==”12″ | prov$CODINE ==”03″ )
pr@data$id = rownames(pr@data)
pr.points = fortify(pr, region=”id”)
pr.df = join(pr.points, pr@data, by=”id”)

ggplot(muni.df) + aes(long,lat,group=group) + geom_path(color=”blue”) +
+ coord_equal()+ geom_path(data=pr.df, +
aes(x=long, y=lat, group=group),color=”red”, size=0.5)

这段代码绘制了一张带有所有边界的漂亮地图enter

对于按级别填充的多边形,我尝试按照 http://tormodboe.wordpress.com/2011/02/22/g?y-med-kart-2/ 中的建议进行读取然后合并

level=read.csv(” levels.dat”,header=”T,sep=””)

munlevel=merge(muni.df,level,by=”CODINE”)

但它给出了一个错误

Error en fix.by(by.x, x) : ‘by’ must specify a uniquely valid column

我不熟悉 shapefile,也许我需要了解更多关于 shp 数据属性的信息,才能找到合并两个数据集的正确选择。如何合并数据以便绘制线条(市政边界)然后用级别填充?

  • 可以在 [gis.stackexchange.com/q/131741/9227] 上找到此问题的更新以及地图上的一些额外功能


[注意:这个问题是在一个多月前提出的,所以 OP 可能已经找到了一种不同的方法来解决他们的问题。我在处理这个相关问题时偶然发现了它。包含此答案是希望对其他人有所帮助。]

这似乎是 OP 要求的…

></p>
<p>... 并使用以下代码生成:</p>
<div class=

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
require(“rgdal”)
require(“maptools”)
require(“ggplot2”)
require(“plyr”)

# read temperature data
setwd(“<location if your data file>”)
temp.data        <- read.csv(file =”levels.dat”, header=TRUE, sep=””, na.string=”NA”, dec=”.”, strip.white=TRUE)
temp.data$CODINE <- str_pad(temp.data$CODINE, width = 5, side = ‘left’, pad = ‘0’)

# read municipality polygons
setwd(“<location of your shapefile”)
esp     <- readOGR(dsn=”.”, layer=”poligonos_municipio_etrs89″)
muni    <- subset(esp, esp$PROVINCIA ==”46″ | esp$PROVINCIA ==”12″ | esp$PROVINCIA ==”3″)
# fortify and merge: muni.df is used in ggplot
muni@data$id <- rownames(muni@data)
muni.df <- fortify(muni)
muni.df <- join(muni.df, muni@data, by=”id”)
muni.df <- merge(muni.df, temp.data, by.x=”CODIGOINE”, by.y=”CODINE”, all.x=T, a..ly=F)
# create the map layers
ggp <- ggplot(data=muni.df, aes(x=long, y=lat, group=group))
ggp <- ggp + geom_polygon(aes(fill=LEVEL))         # draw polygons
ggp <- ggp + geom_path(color=”grey”, linestyle=2)  # draw boundaries
ggp <- ggp + coord_equal()
ggp <- ggp + scale_fill_gradient(low =”#ffffcc”, high =”#ff4444″,
                                 space =”Lab”, na.value =”grey50″,
                                 guide =”colourbar”)
ggp <- ggp + labs(title=”Temperature Levels: Comunitat Valenciana”)
# render the map
print(ggp)

解释:

使用 readOGR(…) 导入 R 的形状文件是 SpacialDataFrame 类型,并有两个主要部分:一个多边形部分,其中包含每个多边形上所有点的坐标,以及一个包含每个多边形信息的数据部分(所以,每个多边形一行)。这些可以被引用,例如,使用 muni@polygons 和 muni@data。实用函数 fortify(…) 将多边形部分转换为用 ggplot 组织的用于绘图的数据框。所以基本的工作流程是:

1
2
3
4
5
6
[1] Import temperature data file (temp.data)
[2] Import polygon shapefile of municipalities (muni)
[3] Convert muni polygons to a data frame for plotting (muni.df <- fortify(…))
[4] Join columns from muni@data to muni.df
[5] Join columns from temp.data to muni.df
[6] Make the plot

连接必须在公共字段上完成,这是大多数问题出现的地方。原始 shapefile 中的每个多边形都有一个唯一的 ID 属性。在 shapefile 上运行 fortify(…) 会创建一个基于此的列 id。但是数据部分没有 ID 列。相反,多边形 ID 存储为行名称。所以首先我们必须向 muni@data 添加一个 id 列,如下所示:

1
muni@data$id <- rownames(muni@data)

现在我们在 muni@data 中有一个 id 字段,在 muni.df 中有一个对应的 id 字段,所以我们可以进行连接:

1
muni.df <- join(muni.df, muni@data, by=”id”)


要创建地图,我们需要根据温度水平设置填充颜色。为此,我们需要将 LEVEL 列从 temp.data 连接到 muni.df。在 temp.data 中有一个字段 CODINE 用于标识自治市。现在,muni.df 中还有一个对应的字段 CODIGOINE。但是有一个问题:CODIGOINE 是 char(5),带有前导零,而 CODINE 是整数,这意味着缺少前导零(可能是从 Excel 导入的?)。因此,仅加入这两个字段不会产生匹配项。我们必须首先将 CODINE 转换为带有前导零的 char(5):

1
temp.data$CODINE <- str_pad(temp.data$CODINE, width = 5, side = ‘left’, pad = ‘0’)

现在我们可以根据相应的字段将 temp.dat 连接到 muni.df。

1
muni.df <- merge(muni.df, temp.data, by.x=”CODIGOINE”, by.y=”CODINE”, all.x=T, a..ly=F)

我们使用 merge(…) 而不是 join(…),因为连接字段具有不同的名称,而 join(…) 要求它们具有相同的名称。 (但请注意,join(…) 更快,应尽可能使用)。因此,最后,我们有一个数据框,其中包含用于绘制多边形的所有信息和温度 LEVEL,可用于为每个多边形建立填充颜色。

OP原代码的一些说明:

  • OP 的第一张地图(顶部的绿色地图)标识了”我们地区的 30 个不同区域……”。我找不到识别这些区域的 shapefile。市政府文件确定了 543 个市镇,我看不出有办法将它们分成 30 个区域。此外,温度级别文件有 542 行,每个市(或多或少)一个。

  • OP 正在为市政当局导入线形文件以绘制边界。您不需要这样做,因为 geom_polygon(…) 将绘制(并填充)多边形,而 geom_path(…) 将绘制边界。

    • 嗨,我没有找到解决方案,但我会试试你的代码。你的地图看起来令人印象深刻,正是我想要的。我不得不暂时离开这个问题,因为我之前还有其他问题要解决。非常感谢您的辛勤工作。
    • 如果这对您有用,请考虑选择作为答案(绿色复选标记)。
    • 由于我的 shapefile,我的脚本有一些问题。我再次下载,现在代码运行完美满足我的需求。干得好@jlhoward
    • 这是太棒了。感谢您提供最有用的解释!
    • 我不得不做一些不同的事情。使用 fortify(muni, region = ‘id’) (否则它不使用 id 变量)和 merge (我猜是 dplyr 的不同版本)。请参阅stackoverflow.com/questions/11052544/…
    • 如果 merge 很慢,您可以使用 dplyr 并使用 left_join(muni.df,temp.data,by=c(“CODIGOINE”=”CODINE”)) which is equivalent to merge(x,y,all.x=TRUE,all.y=FALSE)`


    来源:https://www.codenong.com/19791210/

    微信公众号
    手机浏览(小程序)

    Warning: get_headers(): SSL operation failed with code 1. OpenSSL Error messages: error:14090086:SSL routines:ssl3_get_server_certificate:certificate verify failed in /mydata/web/wwwshanhubei/web/wp-content/themes/shanhuke/single.php on line 57

    Warning: get_headers(): Failed to enable crypto in /mydata/web/wwwshanhubei/web/wp-content/themes/shanhuke/single.php on line 57

    Warning: get_headers(https://static.shanhubei.com/qrcode/qrcode_viewid_8815.jpg): failed to open stream: operation failed in /mydata/web/wwwshanhubei/web/wp-content/themes/shanhuke/single.php on line 57
    0
    分享到:
    没有账号? 忘记密码?