该博客介绍了如何使用R语言进行地理信息系统(GIS)的数据处理，包括读取伦敦边界.shp文件和垃圾事件.csv数据，通过sf包进行数据融合，并使用tmap包绘制地图。作者展示了从设置工作路径、安装所需包到数据读取、合并、绘图及文件导出的完整步骤，最后还涉及了数据格式转换和tmap的交互视图。
https://andrewmaclachlan.github.io/CASA0005repo/geographic-information.html0 数据准备
伦敦边界数据 (.shp)
伦敦2008-2020垃圾乱丢事件统计数据 (.csv)
1 基本设置
1.1 设置工作路径（数据所在路径）
setwd("C:/Users/Desktop/GIS/wk1")1.2 采用pacman的p_load批量安装包
install.packages('pacman')
pacman::p_load(sf,tidyverse,RSQLite,tmap,tmaptools)2 读取shp文件
2.1 sf包st_read函数读取shp文件
shape <- st_read("statistical-gis-boundaries-london/ESRI/London_Borough_Excluding_MHW.shp")#利用sf包读取shp文件
Reading layer `London_Borough_Excluding_MHW' from data source
`C:\Users\liumengjie\Desktop\GIS\wk1\statistical-gis-boundaries-london\ESRI\London_Borough_Excluding_MHW.shp'
using driver `ESRI Shapefile'
Simple feature collection with 33 features and 7 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: 503568.2 ymin: 155850.8 xmax: 561957.5 ymax: 200933.9
Projected CRS: OSGB 1936 / British National Grid2.2 查看文件内容
summary(shape)包括7列数据信息和一列地理位置信息
NAME GSS_CODE HECTARES NONLD_AREA
Length:33 Length:33 Min. : 314.9 Min. : 0.00
Class :character Class :character 1st Qu.: 2724.9 1st Qu.: 0.00
Mode :character Mode :character Median : 3857.8 Median : 2.30
Mean : 4832.4 Mean : 64.22
3rd Qu.: 5658.5 3rd Qu.: 95.60
Max. :15013.5 Max. :370.62
ONS_INNER SUB_2009 SUB_2006 geometry
Length:33 Length:33 Length:33 MULTIPOLYGON :33
Class :character Class :character Class :character epsg:NA : 0
Mode :character Mode :character Mode :character +proj=tmer...: 0 2.3 绘图
plot(shape)根据前7列数据得到7个图
如果只想要边界图
plot(shape$geometry)
3 加载csv数据并合并到shp文件中
3.1 读取csv
mycsv <- read_csv("fly-tipping-borough.csv") Rows: 33 Columns: 15
-- Column specification -------------------------------------------------------------------
Delimiter: ","
chr (1): code
dbl (14): 2007-08, 2008-09, 2009-10, 2010-11, 2011-12, 2012-13, 2013-14, 2014-15, 2015-...
i Use `spec()` to retrieve the full column specification for this data.
i Specify the column types or set `show_col_types = FALSE` to quiet this message.3.2 查看文件
mycsv# A tibble: 33 x 15
code `2007-08` `2008-09` `2009-10` `2010-11` `2011-12` `2012-13` `2013-14` `2014-15`
<chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 E09000001 252 6 31 243 563 1492 433 347
2 E09000002 3296 2746 1865 658 2687 2122 1399 3100
3 E09000003 2415 3952 5876 3564 2356 1828 890 219
4 E09000004 406 286 771 998 505 627 732 399
5 E09000005 5726 5161 4990 5715 6713 2232 3189 3926
6 E09000006 224 223 351 200 306 479 653 462
7 E09000007 14607 10232 9476 8662 5541 5962 8281 4837
8 E09000008 378 390 1490 2580 31 108 460 712
9 E09000009 1845 10874 14341 11104 6727 5679 6543 5521
10 E09000010 4491 3942 9680 9350 7262 3595 3269 4937
# ... with 23 more rows, and 6 more variables: 2015-16 <dbl>, 2016-17 <dbl>,
# 2017-18 <dbl>, 2018-19 <dbl>, 2019-20 <dbl>, total <dbl>3.3 按shape中的GSS_CODE列和mycsv中的code列合并数据
shape <- merge(shape,
mycsv,
by.x="GSS_CODE",
by.y="code")3.4 查看合并后的数据
head(shape,n=10)Simple feature collection with 10 features and 7 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: 503568.2 ymin: 155850.8 xmax: 561957.5 ymax: 198355.2
Projected CRS: OSGB 1936 / British National Grid
NAME GSS_CODE HECTARES NONLD_AREA ONS_INNER SUB_2009 SUB_2006
1 Kingston upon Thames E09000021 3726.117 0.000 F <NA> <NA>
2 Croydon E09000008 8649.441 0.000 F <NA> <NA>
3 Bromley E09000006 15013.487 0.000 F <NA> <NA>
4 Hounslow E09000018 5658.541 60.755 F <NA> <NA>
5 Ealing E09000009 5554.428 0.000 F <NA> <NA>
6 Havering E09000016 11445.735 210.763 F <NA> <NA>
7 Hillingdon E09000017 11570.063 0.000 F <NA> <NA>
8 Harrow E09000015 5046.330 0.000 F <NA> <NA>
9 Brent E09000005 4323.270 0.000 F <NA> <NA>
10 Barnet E09000003 8674.837 0.000 F <NA> <NA>
geometry
1 MULTIPOLYGON (((516401.6 16...
2 MULTIPOLYGON (((535009.2 15...
3 MULTIPOLYGON (((540373.6 15...
4 MULTIPOLYGON (((521975.8 17...
5 MULTIPOLYGON (((510253.5 18...
6 MULTIPOLYGON (((549893.9 18...
7 MULTIPOLYGON (((510599.8 19...
8 MULTIPOLYGON (((510599.8 19...
9 MULTIPOLYGON (((525201 1825...
10 MULTIPOLYGON (((524579.9 19...4 绘图
利用tmap包里的qtm函数快速绘图，根据2011-12列数据进行颜色填充
tmap_mode("plot")
qtm(shape,fill = "2012-13")
5 导出文件
5.1 导出GeoPackage(.gpkg)格式
st_write(shape,"london_1.gpkg",
"london_boroughs_fly_tipping",
delete_layer=TRUE)Deleting layer `london_boroughs_fly_tipping' using driver `GPKG'
Writing layer `london_boroughs_fly_tipping' to data source `london_1.gpkg' using driver `GPKG'
Writing 33 features with 21 fields and geometry type Multi Polygon.5.2 将csv文件添加到Geopackage
con <- dbConnect(RSQLite::SQLite(),dbname="london_1.gpkg")##连接数据RSQLite包
dbListTables(con)##查看数据
dbWriteTable(con,'original csv',
mycsv,
overwrite=TRUE)#写入mycsv
dbListTables(con)##查看数据，多出original csv项
dbDisconnect(con)##断开连接6 summary & Appendix
6.1 全部代码
############################
#step1：设置和安装包
setwd("C:/Users/Desktop/GIS/wk1")##设置工作路径
install.packages('pacman')
pacman::p_load(sf,tidyverse,RSQLite,tmap,tmaptools)
############################
#step2：读取shp数据
shape <- st_read("statistical-gis-boundaries-london/ESRI/London_Borough_Excluding_MHW.shp")#利用sf包读取shp文件
#sf包st_read函数读取shp文件
summary(shape)#查看文件内容，包括7列数据信息和一列地理位置信息
plot(shape)#绘图
plot(shape$geometry)#只绘制地理信息
############################
#step3：加载csv数据
mycsv <- read_csv("fly-tipping-borough.csv")
mycsv##查看文件
shape <- merge(shape,
mycsv,
by.x="GSS_CODE",
by.y="code")##按shape中的GSS_CODE列和mycsv中的code列合并数据
head(shape,n=10)##查看新的shape数据
############################
#step4：绘图
tmap_mode("plot")
qtm(shape,fill = "2012-13")##利用tmap包里的qtm函数快速绘图，根据2011-12列数据进行颜色填充
############################
#step5：导出文件
st_write(shape,"london_1.gpkg",
"london_boroughs_fly_tipping",
delete_layer=TRUE)###导出GeoPackage(.gpkg)格式
#将csv添加到Geopackage
con <- dbConnect(RSQLite::SQLite(),dbname="london_1.gpkg")##连接数据RSQLite包
dbListTables(con)##查看数据
dbWriteTable(con,'original csv',
mycsv,
overwrite=TRUE)#写入mycsv
dbListTables(con)##查看数据，多出original csv项
dbDisconnect(con)##断开连接
############################
##Chris622 20210928#########6.2 sp格式和sf格式互换
class(shape)
shape1 <- as_Spatial(shape)#sf到sp
class(shape1)
shape2 <- st_as_sf(shape1)#sp到sf
class(shape2)> class(shape)
[1] "sf" "data.frame"
> shape1 <- as_Spatial(shape)#sf到sp
> class(shape1)
[1] "SpatialPolygonsDataFrame"
attr(,"package")
[1] "sp"
> shape2 <- st_as_sf(shape1)#sp到sf
> class(shape2)
[1] "sf" "data.frame"6.3 tmap_mode的view模式
tmap_mode("view")#动态格式
qtm(shape,fill = "2012-13")view模式下可交互
6.4 st_write也可以导出shp，geojson等格式
st_write(shape,"london_1.geojson",
"london_boroughs_fly_tipping",
delete_layer=TRUE)
扫一扫
930
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
