1
2
3
4
5
6
7

library(ggplot2)
library(maps)

world_map <- map_data("world")

ggplot(world_map, aes(x = long, y = lat, group = group)) +
  geom_polygon(fill = "lightblue", color = "white")

1
2
3
4
5
6
7
8

library(sp)
library(maptools)

#sids.shp是一个内置的数据，后面会讲到
us_states <- readShapePoly(system.file("shapes/sids.shp", package="maptools")[1],
                           proj4string=CRS("+proj=longlat +datum=WGS84"))

plot(us_states)

1
2
3
4
5

library(leaflet)

leaflet() %>%
  addTiles() %>%
  addMarkers(lng=-74.5, lat=40, popup="New York City")

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19

library(proj4)

# 定义源坐标系（经度-纬度坐标系）
src_crs <- "+proj=longlat +datum=WGS84"
# 定义目标坐标系（局部等面积投影坐标系）
dst_crs <- "+proj=laea +lat_0=40 +lon_0=104"
# 定义待转换的点（经度、纬度）
lon <- 116.4074
lat <- 39.9042

# 转换坐标系
p <- project(cbind(lon, lat), dst_crs)

# 输出转换后的坐标（x、y）
print(p)

# 或者使用sf转换
source_df=st_as_sf(data.frame(lon, lat),coords = c("lon", "lat"), crs = 4326) 
sf::st_transform(source_df,crs="+proj=laea +lat_0=40 +lon_0=104")

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12

# 注意，无法正常运行
library(maps)
library(ggplot2)

# 读取中国地图数据
china_map <- map_data("china")

# 绘制中国地图
ggplot() +
  geom_polygon(data = china_map, aes(x = long, y = lat, group = group), fill = "white", color = "black") +
  theme_void() +
  coord_fixed()

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11

library(sf)
library(ggplot2)
#1.下载数据
dir=getwd()
china_shp=paste0(dir,"/china.json")
if(!file.exists(china_shp))download.file("https://gitcode.net/mirrors/lyhmyd1211/geomapdata_cn/-/raw/master/china.json?inline=false",china_shp)
china <- sf::read_sf(china_shp)

#2.民政部API link: https://mp.weixin.qq.com/s/qj1SRc6D8sgYJYaZzDux6Q
API_pre = "http://xzqh.mca.gov.cn/data/"
china = sf::st_read(dsn = paste0(API_pre, "quanguo.json"), stringsAsFactors=FALSE) 

1

head(china)

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16

## Simple feature collection with 6 features and 10 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 97.1729 ymin: 34.58341 xmax: 126.0685 ymax: 53.33701
## Geodetic CRS:  WGS 84
## # A tibble: 6 × 11
##   adcode name         child…¹ level parent subFe…² adchar center centr…³ acrou…⁴
##   <chr>  <chr>          <int> <chr> <chr>    <int> <chr>  <list> <list>  <list> 
## 1 110000 北京市            16 prov… "{ \"…       0 <NA>   <dbl>  <dbl>   <int>  
## 2 120000 天津市            16 prov… "{ \"…       1 <NA>   <dbl>  <dbl>   <int>  
## 3 130000 河北省            11 prov… "{ \"…       2 <NA>   <dbl>  <dbl>   <int>  
## 4 140000 山西省            11 prov… "{ \"…       3 <NA>   <dbl>  <dbl>   <int>  
## 5 150000 内蒙古自治区      12 prov… "{ \"…       4 <NA>   <dbl>  <dbl>   <int>  
## 6 210000 辽宁省            14 prov… "{ \"…       5 <NA>   <dbl>  <dbl>   <int>  
## # … with 1 more variable: geometry <MULTIPOLYGON [°]>, and abbreviated variable
## #   names ¹​childrenNum, ²​subFeatureIndex, ³​centroid, ⁴​acroutes

1
2
3
4
5
6
7

library(sf)
library(ggplot2)
library(ggspatial)
# plot
ggplot(china)+
  geom_sf()+
   labs(title="Ministry of Civil of PRC",x="Lon",y="Lat") 

1
2
3
4
5
6

#改变crs
# plot
ggplot(china)+
  geom_sf()+ 
  coord_sf(crs = "+proj=laea +lat_0=40 +lon_0=104")+
    labs(title="Ministry of Civil of PRC",x="Lon",y="Lat") 

1
2
3
4
5
6
7
8
9

main_map=ggplot()+
  geom_sf(data = china,fill=NA,size=1,color="black")+ 
  coord_sf(crs = "+proj=laea +lat_0=40 +lon_0=104")+
   labs(title="Ministry of Civil of PRC",x="Lon",y="Lat") +
    annotation_scale(location = "bl") +
      # spatial-aware automagic north arrow
    annotation_north_arrow(location = "tl", which_north = "false",
                             style = north_arrow_fancy_orienteering)+theme_minimal()
main_map

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22

#生成测试数据
anno=data.frame(lon=runif(20,90,120),
                lat=runif(20,30,40),
                type=sample(letters[1:3], 20, replace = TRUE),
                value=runif(20,1,10))
#转换对象
anno_sf <- st_as_sf(anno,coords = c("lon", "lat"),crs = 4326)

p2=ggplot()+
  geom_sf(data = china,fill=NA,size=1,color="black")+
    #添加散点注释
    geom_sf(data = anno_sf,aes(fill=type,size=value),shape=21,colour='black',stroke=.25)+ 
    scale_size(range = c(1,5))+ 
    #添加text
    geom_sf_text(data = anno_sf,aes(label=round(value,1)),size=2,legend="none")+
  coord_sf(crs = "+proj=laea +lat_0=40 +lon_0=104")+
   labs(title="Ministry of Civil of PRC",x="Lon",y="Lat") +
    annotation_scale(location = "bl") +
      # spatial-aware automagic north arrow
    annotation_north_arrow(location = "tl", which_north = "false",
                             style = north_arrow_fancy_orienteering)+theme_minimal()
p2

1
2
3
4
5
6
7
8
9

#例如，这个方法是无法取出90E到120E的部分的，而是很窄的一段
p2+xlim(90,120)

data.frame(lon=c(90,120),lat=c(30,30))%>%
    st_as_sf(coords = c("lon", "lat"), crs = 4326) %>%
    sf::st_transform(df2_sf,crs="+proj=laea +lat_0=40 +lon_0=104")

#发现结果对应的是(-1348988,1539346),这样截取的才是90E到120E的部分。
p2+xlim(-1348988,1539346)

1
2

install.packages("leaflet")
library(leaflet)

 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

# Prepare the text for the tooltip (HTML style):
gre_text <- paste(
  "Type: ", anno$type, "<br/>",
  "Value: ", round(anno$value,2), "<br/>"
  )%>%
  lapply(htmltools::HTML)

type_col=colorFactor(palette = c("red3","green3","blue3"),domain =anno$type)

inter_p=leaflet(anno) %>%
    #添加图层
  addTiles() %>%
    #确定中心点
  setView(lng = 104, lat = 40, zoom = 4)%>%
    #添加散点注释
    addCircleMarkers(~lon, ~lat,
                   fillColor = ~type_col(type), fillOpacity = 0.7, color="white", radius=8, stroke=FALSE,
                   label = gre_text,
                   labelOptions = labelOptions( style = list("font-weight" = "normal", padding = "3px 8px"),
                                                textsize = "13px", direction = "auto")) %>%
    #添加图例
  addLegend( pal=type_col, values=~type, opacity=0.7, title = "Type", position = "bottomright" )

class(inter_p)

inter_p

1
2

#保存结果
htmlwidgets::saveWidget(inter_p,file = "leaflet.html")

 1
 2
 3
 4
 5
 6
 7
 8
 9
10

library(raster)

# 读入一个栅格数据文件
r <- raster("path/to/raster.tif")

# 计算最小值、最大值、均值、方差
min_val <- minValue(r)
max_val <- maxValue(r)
mean_val <- mean(r)
var_val <- var(r)

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15

library(raster)

# 读入两个栅格数据文件
r1 <- raster("path/to/raster1.tif")
r2 <- raster("path/to/raster2.tif")

# 裁剪两个栅格数据文件
r1_crop <- crop(r1, extent(r2))
r2_crop <- crop(r2, extent(r1))

# 合并两个栅格数据文件
r_merge <- merge(r1_crop, r2_crop)

# 将一个栅格数据文件的分辨率重采样为另一个栅格数据文件的分辨率
r_resample <- resample(r1_crop, r2_crop)

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16

library(raster)

# 读入两个栅格数据文件
r1 <- raster("path/to/raster1.tif")
r2 <- raster("path/to/raster2.tif")

# 对两个栅格数据文件进行加、减、乘、除运算
r_add <- r1 + r2
r_sub <- r1 - r2
r_mul <- r1 * r2
r_div <- r1 / r2

# 对一个栅格数据文件进行逻辑与、或、非运算
r_and <- r1 & r2
r_or <- r1 | r2
r_not <- !r1

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13

library(raster)

# 读入一个栅格数据文件
r <- raster("path/to/raster.tif")

# 绘制等高线图
contour(r, add=FALSE)

# 绘制等值线图
image(r, col=terrain.colors(255))

# 绘制渐变色图
plot(r, col=rainbow(255))

1
2
3
4
5

library(raster)

avgtemp <- raster("~/Downloads/wc2.1_10m_tavg/wc2.1_10m_tavg_01.tif")
#全球平均气温数据
plot(avgtemp)

1
2
3

#切片
china_tmp=crop(avgtemp,china)
plot(china_tmp)

1
2
3

#掩膜
china_tmp=mask(china_tmp,china)
plot(china_tmp)

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22

# 将RasterLayer转换为sf对象
sf_data <- st_as_sf(rasterToPolygons(china_tmp), crs = st_crs(china_tmp))
colnames(sf_data)[1]="Average_temperature"

ggplot()+
  geom_sf(data = china,fill=NA,size=1,color="black")+
    #添加栅格数据
    geom_sf(data = sf_data,aes(fill=Average_temperature),shape=22,size=0.5,stroke=0,color=NA)+
    scale_fill_gradient2(low = "blue",mid = "white",high = "red")+
    #使用ggnewscale产生一个新的fill映射，否则会覆盖报错。
    ggnewscale::new_scale_fill()+
    #添加散点注释
    geom_sf(data = anno_sf,aes(fill=type,size=value),shape=21,colour='black',stroke=.25)+ 
    scale_size(range = c(1,5))+ 
    #添加text
    geom_sf_text(data = anno_sf,aes(label=round(value,1)),size=2,legend="none")+
  coord_sf(crs = "+proj=laea +lat_0=40 +lon_0=104")+
   labs(title="Ministry of Civil of PRC",x="Lon",y="Lat") +
    annotation_scale(location = "bl") +
      # spatial-aware automagic north arrow
    annotation_north_arrow(location = "tl", which_north = "false",
                             style = north_arrow_fancy_orienteering)+theme_minimal()