提供了编程的基础技术教程

zl程序教程

您现在的位置是：首页 >  Java

当前栏目

详谈R语言构建地理投影系统绘制高端地图

2023-02-18 16:35:31 时间

❝本节来详细介绍如何使用R语言来构建地理投影系统绘制世界地图，细节挺多的小编做了详细的注释；结果仅供参考❞

加载R包

library(tidyverse)
library(sf)
library(camcorder)

导入数据

world <- read_sf("countries.geojson") %>% janitor::clean_names() %>% 
  rmapshaper::ms_simplify(keep = 0.2)

tomato_prod <- read_csv("tomato-production.txt") %>% janitor::clean_names() %>% 
  rename("tonnes" = last_col()) %>% filter(year == 2020)

数据联接

tomato_world <- world %>% 
  left_join(tomato_prod, by = c("iso_a3" = "code"))

构建投影系统

crs_wintri <- "+proj=wintri +datum=WGS84 +no_defs +over"

❝crs_wintri变量，它包含了一个投影系统的信息。 投影系统用来确定如何在平面上把地球表面的地理空间信息进行投影，以便更好地展示和分析数据。 这个投影系统叫做 "Winkel Tripel"，它是一种广泛用于地图制作的投影方式，可以保证较为均衡地表示地球上的经纬度信息。 ❞

数据转换

tomato_world_wintri <- lwgeom::st_transform_proj(tomato_world, crs = crs_wintri)

❝st_transform_proj 函数并将另一个变量 tomato_world 中的地理空间数据转换到了新的投影系统（即前面定义的 crs_wintri）中。 该函数接受两个参数：tomato_world：这是一个地理空间数据框架，包含了某些地理空间数据（例如地图上的点、线或多边形）。 crs参数，用来指定目标投影系统。在这里，我们将 tomato_world 中的数据转换到了 crs_wintri 指定的投影系统中。 ❞

❝在这个例子中，我们假设 tomato_world 是一个包含世界地图信息的地理空间数据框架，而 crs_wintri 则指定了一种新的投影方式，即 Winkel Tripel 投影。 这样，通过运行上面的代码，我们就能得到一个新的地理空间数据框架 tomato_world_wintri，它将原来的地图信息转换到了新的投影系统中。 ❞

构建经纬线

grat_wintri <- st_graticule(lat = c(-89.9, seq(-80, 80, 20), 89.9)) %>%
  lwgeom::st_transform_proj(crs = crs_wintri)

❝grat_wintri数据包含了一组新的经纬线。使用 st_graticule 函数来生成经纬线。该函数接受一个参数 lat，用来指定经纬线的纬度间隔。 在此将纬度间隔设为 -89.9，seq(-80, 80, 20) 和 89.9。这意味着我们将会在纬度为 -89.9、-80、-60、-40、-20、0、20、40、60、80 和 89.9 的位置绘制经纬线。 st_transform_proj将刚刚生成的经纬线转换到了新的投影系统中，即前面定义的 crs_wintri。这样就得到一个新的地理空间数据框架 grat_wintri，它包含了经过投影转换后的经纬线信息。 ❞

定义经纬度值

lats <- c(90:-90, -90:90, 90)
longs <- c(rep(c(180, -180), each = 181), 180)

❝lats 中包含了从 90 到 -90 的纬度值，并且会有重复的纬度值。这意味着我们会在纬度为 90、-90 和 90 的位置绘制经纬线。 longs 中包含了 180、-180 和 180 三组经度值。由于每一组经度值都有 181 个，因此实际上会在经度为 180、-180 和 180 的位置绘制经纬线。通过这些纬度和经度值，可以在地图上绘制经纬线。 例如，我们可以把每一对纬度和经度值看作一个点，并将这些点连接起来，从而得到一组经纬线。 ❞

构建polygon对象

wintri_outline <- 
  list(cbind(longs, lats)) %>%
  st_polygon() %>%
  st_sfc(crs = "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs") %>% 
  st_sf() %>%
  lwgeom::st_transform_proj(crs = crs_wintri)

❝使用list(cbind(longs, lats))将经度long和纬度lats绑定在一起，得到一个列表对象。 st_polygon函数将这个列表对象转换为一个空间几何(sfc)对象，表示一个多边形。 st_sfc函数将这个空间几何对象包装在一个简单空间几何集合sfc对象中，并为其指定坐标参考系统 st_sf函数将这个简单空间几何集合对象转换为一个空间数据框对象，并使用st_transform_proj函数将其转换为另一个坐标参考系统。 ❞

数据可视化

pal <- RColorBrewer::brewer.pal(8, "RdBu")[1:8]

ggplot(tomato_world_wintri) +
  geom_sf(data = wintri_outline, fill = "#5BBCD6", color = NA) +
  geom_sf(data = grat_wintri, color = "grey20", linewidth = 0.15) +
  geom_sf(aes(fill = tonnes)) +
  scale_fill_stepsn(colors = pal, na.value = "grey97", breaks = c(0, 1e5, 5e5, 1e6, 2.5e6, 5e6, 65e6),
                    labels = scales::label_number(suffix = " M", scale = 1e-6), trans = "log") +
  coord_sf(datum = NULL) +
  guides(fill = guide_colorsteps(title.position = "top", title.hjust = 0.5)) +
  theme_void(base_family = "Outfit") +
  theme(legend.position = "top",
    legend.title = element_text(size = 16),
    legend.key.width = unit(3.5, "lines"),
    legend.key.height = unit(0.5, "lines"),
    plot.background = element_rect(fill = "white", color = NA),
    plot.margin = margin(0, 0, 0, 0))

相关文章