中国地图

> library(maps)
> library(mapdata)
> map("china", col = "blue",  panel.first = grid())
> 
> library(mapchina)

> library(tidyverse)
> library(sf)
> head(china)

## Simple feature collection with 6 features and 13 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: 115.4248 ymin: 39.44473 xmax: 116.8805 ymax: 41.05936
## geographic CRS: WGS 84
## # A tibble: 6 x 14
##   Code_County Code_Perfecture Code_Province Name_Province Name_Perfecture
##   <chr>       <chr>           <chr>         <chr>         <chr>          
## 1 110101      1101            11            <U+5317><U+4EAC><U+5E02>        <NA>           
## 2 110102      1101            11            <U+5317><U+4EAC><U+5E02>        <NA>           
## 3 110114      1101            11            <U+5317><U+4EAC><U+5E02>        <NA>           
## 4 110115      1101            11            <U+5317><U+4EAC><U+5E02>        <NA>           
## 5 110111      1101            11            <U+5317><U+4EAC><U+5E02>        <NA>           
## 6 110116      1101            11            <U+5317><U+4EAC><U+5E02>        <NA>           
## # ... with 9 more variables: Name_County <chr>, Pinyin <chr>, Pop_2000 <dbl>,
## #   Pop_2010 <dbl>, Pop_2017 <dbl>, Pop_2018 <dbl>, Area <dbl>, Density <dbl>,
## #   geometry <MULTIPOLYGON [°]>

> js <- china %>%
+ filter(Code_Province == 32)
> ggplot(js) +      geom_sf(aes(fill = rank(Density))) +
+ scale_fill_distiller(palette = "GnBu", direction = 1) +
+         theme_bw() +
+         theme(legend.position = "bottom")

地统计分析

> library(sp)
> library(gstat)
> data(meuse)#Meuse河流域表层土壤重金属的采样数据
> coordinates(meuse) = ~x+y#构建坐标
> 
> data(meuse.grid)#网格
> coordinates(meuse.grid) <- c("x","y")#网格编码
> meuse.grid <- as(meuse.grid,"SpatialPixelsDataFrame")#空间数据转dataframe
> 
> plot(variogram(log(zinc)~1,meuse,alpha=c(0,45,90,135)),type="l") #各向异性判断， 0 北 90 东

> v1 <- variogram(log(zinc)~1, meuse, width=90, cutoff=1300)#变异函数与自变量的关系，平均数不变模型
> v <- fit.variogram(variogram(log(zinc)~1,meuse), vgm("Sph"))#球形模型拟合
> plot(v1,v)

> sk.fit <- krige(log(zinc)~1,meuse,meuse.grid,v,beta=5.9)#简单克里格插值

## [using simple kriging]

> plot(sk.fit)#插值结果地图展示

来个电影放松下https://www.bilibili.com/video/BV1xX4y1V7wd/?spm_id_from=333.788.recommend_more_video.-1