中国夜间灯光数据的处理和展示

中国夜间灯光数据的处理和展示

在之前的推文:ggplot2 案例:绘制大规模散点图 我展示了处理夜间灯光数据的一种方法,这种方法将 tiff 文件当作图片读入,然后再为每个像素点添加经纬度坐标。这种方法是可行的但是速度慢且不准确。事实上 tiff 文件可以通过 raster 包的 raster 函数读取为 raster 对象再进行处理,本文使用中国夜间灯光数据演示了这种处理方法的使用。

首先可以从知识下载我提供的附件,里面有一个 lightF182013.tiff 文件,就是中国夜间灯光数据了。使用下面的方式读取该数据:

1
2
3
library(raster)
library(sf)
light <- raster('lightF182013.tiff')

light 对象是个 raster 对象:

1
2
3
4
class(light)
#> [1] "RasterLayer"
#> attr(,"package")
#> [1] "raster"

我们可以通过下面的函数查看所有 raster 对象可用的函数:

1
2
3
4
methods(class = "raster")
#> [1] [ [<- anyNA as.matrix as.raster is.na
#> [7] Ops plot print
#> see '?methods' for accessing help and source code

我们先用 plot() 函数作用一下:

1
plot(light)

将 raster 对象转换成 sf 对象会更方便绘图,前面的教程中我们已经介绍了很多关于 ggplot2 + sf 绘图的教程:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
light_point <- rasterToPoints(light, spatial = TRUE) %>% 
st_as_sf()
light_point
#> Simple feature collection with 604701 features and 1 field
#> geometry type: POINT
#> dimension: XY
#> bbox: xmin: 8187500 ymin: 427500 xmax: 15037500 ymax: 7087500
#> epsg (SRID): NA
#> proj4string: +proj=merc +lon_0=0 +lat_ts=0 +x_0=0 +y_0=0 +a=6378137 +b=6378137 [email protected] +units=m +no_defs
#> First 10 features:
#> lightF182013 geometry
#> 1 4 POINT (13722500 7087500)
#> 2 4 POINT (13717500 7082500)
#> 3 4 POINT (13722500 7082500)
#> 4 4 POINT (13727500 7082500)
#> 5 4 POINT (13732500 7082500)
#> 6 4 POINT (13737500 7082500)
#> 7 4 POINT (13742500 7082500)
#> 8 4 POINT (13747500 7082500)
#> 9 4 POINT (13752500 7082500)
#> 10 4 POINT (13757500 7082500)

画中国地图就一定要带上九段线:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
boundary <- st_read('quanguo_Line.geojson', 
stringsAsFactors = FALSE)
ggplot(light_point) +
geom_sf(aes(color = lightF182013), size = 0.01) +
geom_sf(data = subset(boundary,
QUHUADAIMA == "guojiexian"),
aes(geometry = geometry),
color = "white", size = 0.1) +
scale_color_gradientn(
colours = c("#1E1E1E", "#1E1E1E", "#fcfff7", "#f1ff8f", "#f1e93f", "#fabd08")) +
theme_modern_rc(base_family = cnfont,
plot_title_family = cnfont,
subtitle_family = cnfont,
caption_family = cnfont) +
worldtilegrid::theme_enhance_wtg() +
theme(legend.position = "none") +
labs(title = "中国夜间灯光数据",
subtitle = "TidyFriday Project",
caption = "数据来源:Earth Observation Group - Defense Meteorological Satellite Progam, Boulder | ngdc.noaa.gov\n<https://www.ngdc.noaa.gov/eog/dmsp/downloadV4composites.html>")

安徽夜间灯光数据

首先我们需要一份安徽省的地图数据,然后根据安徽地图从中国夜间灯光地图上截取安徽的部分:

1
2
3
4
ah <- st_read('安徽省.json', stringsAsFactors = FALSE) %>% 
st_transform(projection(light))
ahlight <- mask(light, ah)
plot(ahlight)

把 ahlight 转换成 sf 对象处理:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
ahlight_point <- rasterToPoints(ahlight, spatial = TRUE) %>% 
st_as_sf()
ahlight_point
ggplot(ahlight_point) +
geom_sf(data = ah, size = 0.1) +
geom_sf(aes(color = lightF182013), size = 0.01) +
scale_color_gradientn(
colours = c("#1E1E1E", "#1E1E1E", "#fcfff7", "#f1ff8f", "#f1e93f", "#fabd08")) +
theme_modern_rc(base_family = cnfont,
plot_title_family = cnfont,
subtitle_family = cnfont,
caption_family = cnfont) +
worldtilegrid::theme_enhance_wtg() +
theme(legend.position = "none") +
labs(title = "安徽夜间灯光数据",
subtitle = "TidyFriday Project",
caption = "数据来源:Earth Observation Group - Defense Meteorological Satellite Progam, Boulder | ngdc.noaa.gov\n<https://www.ngdc.noaa.gov/eog/dmsp/downloadV4composites.html>")

分组聚合 —— 计算安徽省每个市的平均夜间灯光亮度

使用 aggregate() 函数进行分组聚合即可:

1
2
ahmean <- aggregate(x = ahlight_point, by = ah, FUN = mean)
plot(ahmean)

各个市的平均夜间灯光亮度值:

1
2
3
4
5
6
7
st_drop_geometry(ahmean) %>% 
cbind(st_drop_geometry(ah)) %>%
as_tibble() %>%
dplyr::select(夜间灯光平均亮度 = lightF182013,
行政区划代码 = adcode,
城市名称 = name) %>%
knitr::kable(align = 'c')
夜间灯光平均亮度 行政区划代码 城市名称
13.605723 340100 合肥市
14.274390 340200 芜湖市
10.032258 340300 蚌埠市
12.111821 340400 淮南市
15.977679 340500 马鞍山市
15.301887 340600 淮北市
12.658537 340700 铜陵市
7.637978 340800 安庆市
6.060311 341000 黄山市
9.328982 341100 滁州市
9.383161 341200 阜阳市
9.791738 341300 宿州市
7.397661 341500 六安市
8.640082 341600 亳州市
6.492240 341700 池州市
7.485030 341800 宣城市

再使用 ggplot2 绘图展示:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
ggplot(ahmean) + 
geom_sf(data = ah, size = 0.1) +
geom_sf(aes(fill = lightF182013), size = 0.01) +
scale_fill_gradientn(
colours = c("#fcfff7", "#f1ff8f", "#f1e93f", "#fabd08")) +
theme_modern_rc(base_family = cnfont,
plot_title_family = cnfont,
subtitle_family = cnfont,
caption_family = cnfont) +
worldtilegrid::theme_enhance_wtg() +
labs(title = "安徽省各市平均夜间灯光亮度",
subtitle = "TidyFriday Project",
caption = "数据来源:Earth Observation Group - Defense Meteorological Satellite Progam, Boulder | ngdc.noaa.gov\n<https://www.ngdc.noaa.gov/eog/dmsp/downloadV4composites.html>",
fill = "夜间平均\n灯光亮度")

可以看出安徽合肥附近的皖江城市群发展的明显更好。

为 leaflet 地图添加夜间灯光数据图层

1
2
3
4
5
6
library(leaflet)
library(leafletCN)
leaflet() %>%
amap() %>%
setView(lng = 103, lat = 36, zoom = 4) %>%
addProviderTiles("NASAGIBS.ViirsEarthAtNight2012")

更多 leaflet 的拓展图层可以从这里获得:http://leaflet-extras.github.io/leaflet-providers/preview/

我们再给中国添加国界线:

1
2
3
4
5
6
subset(boundary, QUHUADAIMA == "guojiexian") %>% 
leaflet() %>%
amap() %>%
setView(lng = 103, lat = 36, zoom = 4) %>%
addProviderTiles("NASAGIBS.ViirsEarthAtNight2012") %>%
addPolylines(color = "white", weight = 1)

我们也可以自己添加夜间灯光数据图层:

1
2
3
4
leaflet() %>% 
amap() %>%
setView(lng = 103, lat = 36, zoom = 4) %>%
addRasterImage(light, colors = c("#1E1E1E", "#1E1E1E", "#fcfff7", "#f1ff8f", "#f1e93f", "#fabd08"))

安徽部分:

opacity = 0.8 用于设定透明度:

1
2
3
4
leaflet() %>% 
amap() %>%
setView(lng = 117.283042, lat = 31.86119, zoom = 6) %>%
addRasterImage(ahlight, colors = c("#1E1E1E", "#1E1E1E", "#fcfff7", "#f1ff8f", "#f1e93f", "#fabd08"), opacity = 0.8)

我们再回顾一下之前提到的给 leaflet 地图添加标题的知识点:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
library(htmltools)
tag.map.title <- tags$style(HTML("
.leaflet-control.map-title {
transform: translate(-50%,20%);
position: fixed !important;
left: 50%;
text-align: center;
padding-left: 10px;
padding-right: 10px;
background: rgba(255,255,255,0.75);
font-weight: bold;
font-size: 22px;
}
"))
title <- tags$div(
tag.map.title, HTML("我的微信好友分布")
)

leaflet() %>%
amap() %>%
setView(lng = 103, lat = 36, zoom = 4) %>%
addRasterImage(light, colors = c("#1E1E1E", "#1E1E1E", "#fcfff7", "#f1ff8f", "#f1e93f", "#fabd08")) %>%
addControl(title, position = "topleft", className = "map-title")

作业

今天有作业了!

  • 作业一:对附件中的 ndvi_201912.tiff(中国年度植被指数(NDVI)空间分布数据集) 数据进行探索性分析;
  • 作业二:将 ndvi_201912.tiff 数据绘制在 leaflet 图层上。
  • 附加题:计算中国每个城市的平均夜间灯光亮度数据。

# R

评论

Your browser is out-of-date!

Update your browser to view this website correctly. Update my browser now

×