ggplot2---绘制地图

最新推荐文章于 2024-09-08 12:11:39 发布

osummer3148

最新推荐文章于 2024-09-08 12:11:39 发布

阅读量1.4k

点赞数

分类专栏： # tidyverse 文章标签： ggplot2 地图

本文链接：https://blog.csdn.net/nikang3148/article/details/122109567

版权

tidyverse 专栏收录该内容

7 篇文章 2 订阅

订阅专栏

文章目录

sf包介绍
加载包
空间数据读取
数据操作
可视化

sf包介绍

主要使用sf包进行空间数据分析。
sf包主要控制simple feature对象，比之前的sp包从数据结构上看更容易懂，也能对空间数据进行大量的操作。

sf包中的对象主要有：

sf 数据框
sfc 列表列
sfg sfc中的任一元素
crs 坐标参考系统

加载包

library(sf)
library(tidyverse)

空间数据读取

以中国省级行政划分为例，数据来源于aliyun的datav。

china_dat <- read_sf("china_province.json")

attr(china_dat, "sf_column")

## [1] "geometry"

class(china_dat)

## [1] "sf"         "tbl_df"     "tbl"        "data.frame"

数据操作

gfc数据绑定

sf对象中有一列（gfc）,专门存贮空间属性数据，在对sf对象数据进行操作时，默认是不能排除gfc列的。

china_dat %>% select(1) %>% knitr::kable()#绑定

Name_Province	geometry
安徽省	POLYGON ((119.1702 31.2955,…
澳门特别行政区	MULTIPOLYGON (((113.5348 22…
北京市	POLYGON ((115.961 40.2643, …
福建省	MULTIPOLYGON (((119.4823 25…
甘肃省	POLYGON ((104.0702 33.67677…
广东省	MULTIPOLYGON (((113.6739 22…
广西壮族自治区	MULTIPOLYGON (((109.083 21…
贵州省	MULTIPOLYGON (((108.302 25…
海南省	MULTIPOLYGON (((110.2262 20…
河北省	MULTIPOLYGON (((117.667 39…
河南省	MULTIPOLYGON (((115.4948 35…
黑龙江省	POLYGON ((124.0832 50.54255…
湖北省	POLYGON ((109.951 31.47014,…
湖南省	MULTIPOLYGON (((112.8936 29…
吉林省	POLYGON ((129.0903 43.54665…
江苏省	MULTIPOLYGON (((118.8514 31…
江西省	POLYGON ((118.2324 28.92183…
辽宁省	MULTIPOLYGON (((122.7548 39…
南海九段线	MULTILINESTRING ((108.338 7…
南海诸岛	MULTIPOLYGON (((119.6371 23…
内蒙古自治区	POLYGON ((102.8302 42.10144…
宁夏回族自治区	MULTIPOLYGON (((106.0464 35…
青海省	POLYGON ((92.40531 39.03601…
山东省	MULTIPOLYGON (((120.1647 35…
山西省	POLYGON ((110.9079 38.7014,…
陕西省	POLYGON ((109.62 33.27277, …
上海市	MULTIPOLYGON (((121.5172 31…
四川省	POLYGON ((103.7797 28.52829…
台湾省	MULTIPOLYGON (((120.4145 24…
天津市	POLYGON ((116.8757 39.07094…
西藏自治区	POLYGON ((79.15081 32.48816…
香港特别行政区	MULTIPOLYGON (((114.1675 22…
新疆维吾尔自治区	POLYGON ((80.77494 43.29664…
云南省	POLYGON ((98.69413 28.21893…
浙江省	MULTIPOLYGON (((121.9817 29…
重庆市	POLYGON ((107.2248 30.22737…

china_dat %>% as_tibble() %>%   select(1) #绑定解除

## # A tibble: 36 x 1
##    Name_Province 
##    <chr>         
##  1 安徽省        
##  2 澳门特别行政区
##  3 北京市        
##  4 福建省        
##  5 甘肃省        
##  6 广东省        
##  7 广西壮族自治区
##  8 贵州省        
##  9 海南省        
## 10 河北省        
## # ... with 26 more rows

可见直接select操作中，gfc列仍然存在。
如果需要去除这种效应，需要把sf类型转换为tibble、或者data.frame类型。

gfc数据操作

china_dat %>% st_geometry() #提取geom列

## Geometry set for 36 features 
## Geometry type: GEOMETRY
## Dimension:     XY
## Bounding box:  xmin: 73.50078 ymin: 3.553559 xmax: 134.7752 ymax: 53.56082
## Geodetic CRS:  WGS 84
## First 5 geometries:

## POLYGON ((119.1702 31.2955, 119.1528 31.29677, ...

## MULTIPOLYGON (((113.5348 22.19369, 113.5325 22....

## POLYGON ((115.961 40.2643, 115.9576 40.25636, 1...

## MULTIPOLYGON (((119.4823 25.35403, 119.4669 25....

## POLYGON ((104.0702 33.67677, 104.0789 33.67631,...

china_dat %>% st_geometry()  %>% class() #sfc

## [1] "sfc_GEOMETRY" "sfc"

china_dat %>% st_geometry()  %>%  attributes() #sfc的属性

## $n_empty
## [1] 0
## 
## $crs
## Coordinate Reference System:
##   User input: WGS 84 
##   wkt:
## GEOGCRS["WGS 84",
##     DATUM["World Geodetic System 1984",
##         ELLIPSOID["WGS 84",6378137,298.257223563,
##             LENGTHUNIT["metre",1]]],
##     PRIMEM["Greenwich",0,
##         ANGLEUNIT["degree",0.0174532925199433]],
##     CS[ellipsoidal,2],
##         AXIS["geodetic latitude (Lat)",north,
##             ORDER[1],
##             ANGLEUNIT["degree",0.0174532925199433]],
##         AXIS["geodetic longitude (Lon)",east,
##             ORDER[2],
##             ANGLEUNIT["degree",0.0174532925199433]],
##     ID["EPSG",4326]]
## 
## $class
## [1] "sfc_GEOMETRY" "sfc"         
## 
## $precision
## [1] 0
## 
## $bbox
##       xmin       ymin       xmax       ymax 
##  73.500775   3.553559 134.775196  53.560815 
## 
## $classes
##  [1] "POLYGON"         "MULTIPOLYGON"    "POLYGON"         "MULTIPOLYGON"   
##  [5] "POLYGON"         "MULTIPOLYGON"    "MULTIPOLYGON"    "MULTIPOLYGON"   
##  [9] "MULTIPOLYGON"    "MULTIPOLYGON"    "MULTIPOLYGON"    "POLYGON"        
## [13] "POLYGON"         "MULTIPOLYGON"    "POLYGON"         "MULTIPOLYGON"   
## [17] "POLYGON"         "MULTIPOLYGON"    "MULTILINESTRING" "MULTIPOLYGON"   
## [21] "POLYGON"         "MULTIPOLYGON"    "POLYGON"         "MULTIPOLYGON"   
## [25] "POLYGON"         "POLYGON"         "MULTIPOLYGON"    "POLYGON"        
## [29] "MULTIPOLYGON"    "POLYGON"         "POLYGON"         "MULTIPOLYGON"   
## [33] "POLYGON"         "POLYGON"         "MULTIPOLYGON"    "POLYGON"

#china_dat %>% st_geometry() %>% .[[3]]  #提取第3行，sfg对象

china_dat %>% st_geometry() %>% .[2]  #提取第2行，sfc对象

## Geometry set for 1 feature 
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 113.5325 ymin: 22.11169 xmax: 113.5983 ymax: 22.21531
## Geodetic CRS:  WGS 84

## MULTIPOLYGON (((113.5348 22.19369, 113.5325 22....

china_dat %>% st_geometry() %>% sapply(length) # 每一行中的多边形数量

##  [1]  2  4  1 45  2 55  5  2  9  4  3  1  3  2  1  4  1 15 15 96  1  2  1 10  1
## [26]  1 10  1 14  3  1 11  1  1 81  1

坐标CRS数据

一般空间数据读取会自带坐标CRS参数,否则需要制定CRS。CRS操作主要有3类：

st_crs 获取CRS
st_set_crs 设置CRS，gfc数据不变，但st_crs返回的参数发生改变
st_transform 转换CRS，gfc数据发生变化

china_dat %>% st_geometry()%>% attributes()

## $n_empty
## [1] 0
## 
## $crs
## Coordinate Reference System:
##   User input: WGS 84 
##   wkt:
## GEOGCRS["WGS 84",
##     DATUM["World Geodetic System 1984",
##         ELLIPSOID["WGS 84",6378137,298.257223563,
##             LENGTHUNIT["metre",1]]],
##     PRIMEM["Greenwich",0,
##         ANGLEUNIT["degree",0.0174532925199433]],
##     CS[ellipsoidal,2],
##         AXIS["geodetic latitude (Lat)",north,
##             ORDER[1],
##             ANGLEUNIT["degree",0.0174532925199433]],
##         AXIS["geodetic longitude (Lon)",east,
##             ORDER[2],
##             ANGLEUNIT["degree",0.0174532925199433]],
##     ID["EPSG",4326]]
## 
## $class
## [1] "sfc_GEOMETRY" "sfc"         
## 
## $precision
## [1] 0
## 
## $bbox
##       xmin       ymin       xmax       ymax 
##  73.500775   3.553559 134.775196  53.560815 
## 
## $classes
##  [1] "POLYGON"         "MULTIPOLYGON"    "POLYGON"         "MULTIPOLYGON"   
##  [5] "POLYGON"         "MULTIPOLYGON"    "MULTIPOLYGON"    "MULTIPOLYGON"   
##  [9] "MULTIPOLYGON"    "MULTIPOLYGON"    "MULTIPOLYGON"    "POLYGON"        
## [13] "POLYGON"         "MULTIPOLYGON"    "POLYGON"         "MULTIPOLYGON"   
## [17] "POLYGON"         "MULTIPOLYGON"    "MULTILINESTRING" "MULTIPOLYGON"   
## [21] "POLYGON"         "MULTIPOLYGON"    "POLYGON"         "MULTIPOLYGON"   
## [25] "POLYGON"         "POLYGON"         "MULTIPOLYGON"    "POLYGON"        
## [29] "MULTIPOLYGON"    "POLYGON"         "POLYGON"         "MULTIPOLYGON"   
## [33] "POLYGON"         "POLYGON"         "MULTIPOLYGON"    "POLYGON"

china_dat %>% st_crs() #默认EPSG4326

## Coordinate Reference System:
##   User input: WGS 84 
##   wkt:
## GEOGCRS["WGS 84",
##     DATUM["World Geodetic System 1984",
##         ELLIPSOID["WGS 84",6378137,298.257223563,
##             LENGTHUNIT["metre",1]]],
##     PRIMEM["Greenwich",0,
##         ANGLEUNIT["degree",0.0174532925199433]],
##     CS[ellipsoidal,2],
##         AXIS["geodetic latitude (Lat)",north,
##             ORDER[1],
##             ANGLEUNIT["degree",0.0174532925199433]],
##         AXIS["geodetic longitude (Lon)",east,
##             ORDER[2],
##             ANGLEUNIT["degree",0.0174532925199433]],
##     ID["EPSG",4326]]

china_dat %>% st_set_crs(3857) #更改CRS属性，不计算，gfc不变

## Warning: st_crs<- : replacing crs does not reproject data; use st_transform for
## that

## Simple feature collection with 36 features and 2 fields
## Geometry type: GEOMETRY
## Dimension:     XY
## Bounding box:  xmin: 73.50078 ymin: 3.553559 xmax: 134.7752 ymax: 53.56082
## Projected CRS: WGS 84 / Pseudo-Mercator
## # A tibble: 36 x 3
##    Name_Province  Code_Province                                         geometry
##  * <chr>          <chr>                                           <GEOMETRY [m]>
##  1 安徽省         34            POLYGON ((119.1702 31.2955, 119.1528 31.29677, ~
##  2 澳门特别行政区 82            MULTIPOLYGON (((113.5348 22.19369, 113.5325 22.~
##  3 北京市         11            POLYGON ((115.961 40.2643, 115.9576 40.25636, 1~
##  4 福建省         35            MULTIPOLYGON (((119.4823 25.35403, 119.4669 25.~
##  5 甘肃省         62            POLYGON ((104.0702 33.67677, 104.0789 33.67631,~
##  6 广东省         44            MULTIPOLYGON (((113.6739 22.66379, 113.6655 22.~
##  7 广西壮族自治区 45            MULTIPOLYGON (((109.083 21.0469, 109.0866 21.03~
##  8 贵州省         52            MULTIPOLYGON (((108.302 25.51802, 108.3148 25.5~
##  9 海南省         46            MULTIPOLYGON (((110.2262 20.06673, 110.2322 20.~
## 10 河北省         13            MULTIPOLYGON (((117.667 39.38564, 117.663 39.41~
## # ... with 26 more rows

china_dat %>% st_transform(3857) #重新投影，gfc变化

## Simple feature collection with 36 features and 2 fields
## Geometry type: GEOMETRY
## Dimension:     XY
## Bounding box:  xmin: 8182069 ymin: 395834.3 xmax: 15003110 ymax: 7087415
## Projected CRS: WGS 84 / Pseudo-Mercator
## # A tibble: 36 x 3
##    Name_Province  Code_Province                                         geometry
##  * <chr>          <chr>                                           <GEOMETRY [m]>
##  1 安徽省         34            POLYGON ((13265965 3671186, 13264033 3671350, 1~
##  2 澳门特别行政区 82            MULTIPOLYGON (((12638632 2534796, 12638384 2533~
##  3 北京市         11            POLYGON ((12908715 4904424, 12908338 4903267, 1~
##  4 福建省         35            MULTIPOLYGON (((13300705 2919292, 13298992 2918~
##  5 甘肃省         62            POLYGON ((11585041 3985482, 11586011 3985421, 1~
##  6 广东省         44            MULTIPOLYGON (((12654123 2591410, 12653183 2593~
##  7 广西壮族自治区 45            MULTIPOLYGON (((12143065 2397471, 12143470 2395~
##  8 贵州省         52            MULTIPOLYGON (((12056128 2939507, 12057548 2941~
##  9 海南省         46            MULTIPOLYGON (((12270323 2280938, 12270989 2280~
## 10 河北省         13            MULTIPOLYGON (((13098626 4777063, 13098181 4780~
## # ... with 26 more rows

st_transform转换后gfc数据发生变化，则可视化地图也会发生改变。

数据有效性

st_is_valid检查gfc数据是否有效，否则可以用st_make_valid进行修正。

st_is_valid(china_dat)

##  [1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [16] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [31] TRUE TRUE TRUE TRUE TRUE TRUE

china_dat %>% st_is_valid(reason=T)

##  [1] "Valid Geometry" "Valid Geometry" "Valid Geometry" "Valid Geometry"
##  [5] "Valid Geometry" "Valid Geometry" "Valid Geometry" "Valid Geometry"
##  [9] "Valid Geometry" "Valid Geometry" "Valid Geometry" "Valid Geometry"
## [13] "Valid Geometry" "Valid Geometry" "Valid Geometry" "Valid Geometry"
## [17] "Valid Geometry" "Valid Geometry" "Valid Geometry" "Valid Geometry"
## [21] "Valid Geometry" "Valid Geometry" "Valid Geometry" "Valid Geometry"
## [25] "Valid Geometry" "Valid Geometry" "Valid Geometry" "Valid Geometry"
## [29] "Valid Geometry" "Valid Geometry" "Valid Geometry" "Valid Geometry"
## [33] "Valid Geometry" "Valid Geometry" "Valid Geometry" "Valid Geometry"

数据变形

基于gfc数据计算距离、面积等。

st_distance(china_dat[c(1,5,10),],china_dat[c(1,8,14,30),]) #求距离

## Units: [m]
##          [,1]     [,2]     [,3]     [,4]
## [1,]      0.0 667571.5 221686.7 440083.9
## [2,] 627042.3 434023.6 505780.6 760186.8
## [3,] 185372.5 954865.4 701016.8      0.0

st_area(china_dat[1,]) #求面积

## 140266534479 [m^2]

st_centroid(china_dat[15,]) #计算重心

## Warning in st_centroid.sf(china_dat[15, ]): st_centroid assumes attributes are
## constant over geometries of x

## Simple feature collection with 1 feature and 2 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: 126.2348 ymin: 43.66413 xmax: 126.2348 ymax: 43.66413
## Geodetic CRS:  WGS 84
## # A tibble: 1 x 3
##   Name_Province Code_Province             geometry
## * <chr>         <chr>         <POINT [arc_degree]>
## 1 吉林省        22             (126.2348 43.66413)

可视化

ggplot结合geom_sf进行绘图。

china_dat %>% ggplot()+geom_sf()

在这里插入图片描述

CRS变换效果

library(patchwork)

## Warning: package 'patchwork' was built under R version 4.1.1

p1=china_dat %>% ggplot()+geom_sf()+labs(title="EPSG4326")
p2=china_dat %>% st_transform(4508) %>% ggplot()+geom_sf()+labs(title="EPSG4508") #CGCS2000中国地图投影效果，这也是大多数纸质出版物地图的国家标准
p1+p2

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-NjPKWFqq-1640245754579)(R-map_files/figure-html/unnamed-chunk-9-1.png)]

空间数据构造及绘图

st_point(c(105,30)) %>% st_sfc() %>%  st_set_crs(4326) #坐标点

## Geometry set for 1 feature 
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: 105 ymin: 30 xmax: 105 ymax: 30
## Geodetic CRS:  WGS 84

## POINT (105 30)

aa=data.frame(long=seq(80,130,by=10),lat=seq(10,60,by=10)) 

bb=aa %>% st_as_sf(coords=c("long","lat"),crs=4326) #构造经纬度点

st_crs(bb)

## Coordinate Reference System:
##   User input: EPSG:4326 
##   wkt:
## GEOGCRS["WGS 84",
##     DATUM["World Geodetic System 1984",
##         ELLIPSOID["WGS 84",6378137,298.257223563,
##             LENGTHUNIT["metre",1]]],
##     PRIMEM["Greenwich",0,
##         ANGLEUNIT["degree",0.0174532925199433]],
##     CS[ellipsoidal,2],
##         AXIS["geodetic latitude (Lat)",north,
##             ORDER[1],
##             ANGLEUNIT["degree",0.0174532925199433]],
##         AXIS["geodetic longitude (Lon)",east,
##             ORDER[2],
##             ANGLEUNIT["degree",0.0174532925199433]],
##     USAGE[
##         SCOPE["Horizontal component of 3D system."],
##         AREA["World."],
##         BBOX[-90,-180,90,180]],
##     ID["EPSG",4326]]

china_dat %>% ggplot()+geom_sf()+
  geom_sf(data=bb,aes(color="red"))

在这里插入图片描述

图形微调

局部放大

#p1为主图

sub_china<-p1+
  coord_sf(xlim=c(105,125),ylim=c(3,30))+
  labs(x="",y="",title="")+
  theme_void()+
  theme(panel.border = element_rect(fill="NA"),legend.position='none')
  

sub_china

在这里插入图片描述

南海诸岛小图

main_china<-p1+
  coord_sf(ylim=c(18,55),xlim=c(70,140))+
  labs(x="Longitude",y="Latitude",fill="Pops")

main_china

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-NHLW7VPH-1640245754581)(R-map_files/figure-html/unnamed-chunk-12-1.png)]

main_china + inset_element(sub_china,0.7, 0.02, 0.99, 0.45)

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-2v4V3qcn-1640245754582)(R-map_files/figure-html/unnamed-chunk-12-2.png)]

1为主图

sub_china<-p1+
coord_sf(xlim=c(105,125),ylim=c(3,30))+
labs(x="",y="",title="")+
theme_void()+
theme(panel.border = element_rect(fill=“NA”),legend.position=‘none’)

sub_china


[外链图片转存中...(img-hJQkZmhV-1640245754581)]<!-- -->




### 南海诸岛小图


```r
main_china<-p1+
  coord_sf(ylim=c(18,55),xlim=c(70,140))+
  labs(x="Longitude",y="Latitude",fill="Pops")

main_china + inset_element(sub_china,0.7, 0.02, 0.99, 0.45)

在这里插入图片描述

osummer3148

关注

0
点赞
踩
7

收藏

觉得还不错? 一键收藏
0
评论
复制链接

分享到 QQ

分享到新浪微博

扫一扫

专栏目录