画地图这件事,永远是看上去越简单、做起来越麻烦。我最早接触R语言绘图,就是想复现一张漂亮的南美洲地图:十二个国家、绵长的安第斯山脉、亚马逊雨林里密布的采样点,全都要落到一张图上。那时候以为调用一个函数就完事,结果光是数据格式就折腾了一个下午。后来做多了才发现,R语言绘制南美洲地图的整个流程——从地理数据获取、投影变换、点数据标记,到配色、版面、主题润色——其实是一条环环相扣的技术链路,踩坑和惊喜一样多。这篇文章我把整套方法拆开揉碎写给你,包括R语言安装环境的基础准备、sf与ggplot2的组合使用、nc文件读取和生信指标的叠加,以及把SARIMA模型预测结果画成时空地图的进阶玩法。无论你是刚入门的学生,还是已经在用R做生信、气象或环境数据分析的从业者,照着这篇文章走一遍,都能画出一张能放到论文里的南美洲地图。
1. 为什么是R语言和南美洲地图
1.1 一张地图背后的数据链路
很多人以为画地图是GIS的专属工作,ArcGIS和QGIS确实能干,但它们的学习成本不低,而且和数据分析的衔接非常别扭。R语言不一样,它本身就是为数据分析而生的,ggplot2把可视化语法统一起来,sf包又让空间对象能以数据框的形式参与一切常规操作。也就是说,你可以在一套工作流里完成三件事:处理业务数据(比如计算α多样性指数、对空气污染时间序列建模)、组织地理边界(国界、州界、河流湖泊)、让两者在地图上碰撞出直观的结果。
南美洲这个区域尤其适合用来练手。它跨赤道南北,地形东西差异巨大,西边是狭长的安第斯山脉,东边是巴西高原和亚马逊雨林,国界形状又极其不规则。这种地块对投影方式特别敏感,画出来丑不丑、变形严重不严重,一眼就能看出来。南方医科大学那边每年都有同学做大作业时被这类地理可视化任务卡住,我见过不少拿着R语言教程翻来覆去就是画不出完整南美洲地图的案例。问题多半不在函数不会用,而是没有理解整条数据链路。
1.2 南美洲地图数据从哪里来
先解决一个根本问题:南美洲的边界数据从哪里来。这里我把常用的几个数据源做一次梳理,你根据自己的精度需求和场景选就行。
| 数据源 | 获取方式 | 精度 | 适合场景 |
|---|---|---|---|
| rnaturalearth | ne_countries() | 中低精度,体积小 | 快速出图、教学演示 |
| rnaturalearthhires | ne_download() | 高精度海岸线 | 论文出版级底图 |
| GADM | raster::getData()或在线下载 | 精细到省/州 | 需要划分行政边界时 |
| geodata | geodata::gadm() | 精细行政 | 生物地理学、物种分布 |
| OpenStreetMap数据 | osmdata包 | 道路河流极细 | 局部放大图 |
如果你只想要一张干净的南美洲轮廓,rnaturalearth是最省事的。它在后台把Natural Earth的矢量数据整理成sf对象,一个函数就能拿到国界、海岸线和国名的组合体。我的习惯是先用它快速验证配色和布局,最后出终稿时才换高精度的GADM或rnaturalearthhires。
还有一类数据不是矢量边界,而是网格化的栅格数据,最常见的就是.nf结尾的NetCDF气候文件。后面我会专门讲怎么把这类数据读进R并与南美洲地图对齐。一句话总结:地图数据没有所谓“最好”,只有“够用”和“适合当前场景”。
2. 环境准备:从零搭好R绘图工作台
2.1 R与RStudio的安装细节
这一步看起来太基础了,但恰恰是新手翻车最多的地方。R语言官网下载安装包时,建议选择R 4.x以上版本,不要在老旧版本上死磕,因为sf、rnaturalearth这些包的新版本对R的最低版本有要求。下载后安装时,Windows用户记得把“添加到系统PATH”选上,虽然RStudio可以自动发现R,但后续用命令行工具调Rscript时会少很多麻烦。
RStudio的版本也值得较真,最新版的Positron和RStudio都支持sf的可视化预览,但如果你手头项目依赖旧版shiny或tidyverse系列,最好保持稳定版RStudio不要手贱点“更新全部包”,这是我踩了几次坑换来的教训。
更关键的一步是Windows用户必须安装Rtools。很多R包(包括sf)在Windows上需要本地编译,没有Rtools会直接报错“无法编译”。macOS用户则需要确保系统里存在Xcode Command Line Tools,Linux用户建议先安装libgdal、libgeos、libproj这些底层地理库。sf包的核心依赖就是这几个C++库,它们不装好,后面所有空间操作都会罢工。
2.2 核心包安装顺序与依赖处理
绘制南美洲地图的标准工具箱并不复杂,下面这行命令可以一次性装齐:
install.packages(c( "sf", "ggplot2", "rnaturalearth", "rnaturalearthdata", "ggspatial", "ggrepel", "RColorBrewer", "viridis", "scales", "ncdf4", "terra", "vegan", "forecast", "patchwork" ), repos = "https://mirrors.tuna.tsinghua.edu.cn/CRAN/")从国内环境安装时,建议把repos参数设置为清华或中科大镜像,速度会快几个量级。装rnaturalearthdata时它会被一并拉下来,因为它本质上是rnaturalearth的低精度数据配套包。
我在实际项目里还会用到qgisprocess调用QGIS做复杂拓扑运算,但纯粹画南美洲地图的话完全没有必要。ggplot2负责画图,sf负责空间数据结构和变换,ggspatial提供比例尺与指北针,ggrepel解决采样点标签覆盖,viridis提供色盲友好的渐变配色。这七个包搭配起来,已经能覆盖从入门到出版级的大部分需求。
如果你在安装sf时遇到“namespace 'sf' is being loaded”老版本包冲突,不要硬撑,果断重启R会话,然后分开逐个安装:先安装sf,再安装ggplot2,最后装其他依赖。这个顺序可以避免tidyverse系列与sf版本之间的潜在冲突。
3. 获取与清洗南美洲的地理数据
3.1 用rnaturalearth获取国界与海岸线
获取南美洲边界数据非常直接:
library(sf) library(rnaturalearth) sa <- ne_countries(scale = "medium", continent = "South America", returnclass = "sf") plot(sa["name"])跑完这段代码你会得到一张带国名的简单地图。但注意,ne_countries返回的对象里其实不止南美,它还包含法属圭亚那、福克兰群岛等。如果你只想画大陆主体,可以手动过滤:
sa_clean <- sa[!sa$name %in% c("Falkland Is.", "Falkland Islands"), ]这个过滤逻辑经常被忽略,导致画出来的地图在多出几个小岛之后构图变得失衡。还有一个小技巧:用sf::st_geometry_type()检查一下数据是否都是MultiPolygon,如果有POLYGON混在里面,后续做合并和裁剪时会被坑。
rnaturalearth最大的便利在于它把国名和ISO代码作为属性字段存储在sf对象里,这让后续按国家分组填色变得极其自然,不用再费劲去匹配外键。
3.2 投影方式如何选:从经纬度到等面积投影
投影是地图绘制中最影响观感的一步,也是概念最多的一步。简单说,地球是个球面,屏幕是平面,把球面压成平面必然变形。R里默认的EPSG:4326就是经纬度坐标,它其实不算一个适合展示的投影,画出来的南美洲会显得又扁又歪。
南美洲东西跨度不大但南北跨度极大,从北纬12度一直到南纬55度,因此非常推荐使用Albers等面积投影(Albers Equal Area Conic)。等面积的意思是图上每一块区域的比例关系都正确,这对生态学、流行病学上有实际意义的统计图特别重要。很多气候模型在做区域平均时也是基于这类投影。
具体转换代码如下:
sa_aea <- st_transform(sa_clean, crs = "+proj=aea +lat_1=5 +lat_2=-42 +lat_0=-32 +lon_0=-60 +datum=WGS84")这几行参数的意思是标准纬线设置在5度和-42度,中心经度-60度,正好覆盖南美洲主体。转完之后你会发现南美洲的形状一下子“挺拔”了,安第斯山脉和海岸线的弯曲关系也更接近地图集里的感觉。
如果你的图里要叠加海洋专题数据(比如海表温度),就需要另一个投影体系了,比如Plate Carrée或墨卡托。投影是关键决策,必须在画图之前定好,不要等全部图层画完再转变换,那样叠加的采样点和边界会反复横跳,消耗的调试时间足以让你怀疑人生。
3.3 读入nc文件和点数据的坐标标准化
再说说很多环境专业同学必然遇到的一类文件——NetCDF,后缀一般就是.nc。打开.nc文件的第一步检查内容:
library(ncdf4) nc <- nc_open("era5_temp_2023.nc") print(nc)print输出会列出所有维度、变量和属性,你要找到纬度和经度对应的变量名。大多数时候名字就叫latitude和longitude,但也有人叫lat和lon,甚至是y和x。拿到之后提取:
lat <- ncvar_get(nc, "latitude") lon <- ncvar_get(nc, "longitude") temp <- ncvar_get(nc, "temperature")这里有个排序的坑:lat可能是降序(从北到南),lon也可能按0到360排列而不是-180到180。南美洲的经度范围是-90到-30,如果数据显示成0到360,你需要做一次换算:lon[lon > 180] <- lon[lon > 180] - 360。
点数据标记在操作中经常是CSV或Excel表格,里面是采样点的经纬度。读取后第一步不是画图,而是先检查坐标合法性:
pts <- read.csv("sites.csv") stopifnot( all(!is.na(pts$longitude)), all(pts$longitude >= -180 & pts$longitude <= 180), all(pts$latitude >= -90 & pts$latitude <= 90) )如果检查不通过,最常见的错是把经纬度写反了,或者数值里有文本格式残留。用stopifnot强制在进入绘图前暴露问题,能省下大量排查时间。最后用st_as_sf把普通数据框转成空间点对象,指定crs=4326,再按需投影到aea。
4. 基础地图:把南美洲画出来
4.1 第一张地图:国界线与陆地区域
进入绘图环节,先用一行代码画最基础的南美大陆:
library(ggplot2) p_base <- ggplot(sa_aea) + geom_sf(fill = "#f5f0e6", color = "#444444", linewidth = 0.4) + theme_void() ggsave("south_america_base.png", p_base, width = 8, height = 10, dpi = 300)这个版本干净但很素。如果想要更精细的层次感,可以在地形上做文章:用terra包读一张全球DEM高程数据,裁剪到南美范围,再用geom_raster绘制高程渐变,最后把国界用geom_sf叠加在上面。这种做法能把安第斯山脉清晰地表现出来,观感立刻提升一个档次。
不过做这个操作时要小心内存占用,全球DEM在30弧秒精度下有上亿个像元,直接读取可能让电脑死机。稳妥做法是先读入全球数据,再用terra::crop和terra::mask裁剪到南美,如果还嫌卡就把分辨率聚合到0.1度再转成数据框。
4.2 叠加河流、湖泊和城市点数据标记
“r语言点数据标记”是很多人搜索的高频需求,这里我给出一套完整实用的做法。假设你有一个采样点文件sites_sf,里面存了物种观测点的经纬度,你想把这些点标在地图上,并且用颜色区分属性:
ggplot(sa_aea) + geom_sf(fill = "#dfe6df", color = "#888888", linewidth = 0.3) + geom_sf(data = sites_projected, aes(color = diversity_group), size = 2.2, alpha = 0.85) + scale_color_brewer(palette = "Set2", name = "多样性分组") + theme_void() + theme(legend.position = "bottom")这里的关键在于sites_projected也必须使用与底图相同的投影。如果你在底图上用的sa_aea是aea投影,而站点对象还是原始经纬度,直接用geom_sf叠加时R会自动做投影转换,理论上没问题,但数据量很大时会导致重投影耗时,而且某些边缘点会因拓扑错误而消失。我的经验是先统一转换再画,也就是:
sites_projected <- st_transform(sites_sf, crs = st_crs(sa_aea))点与点之间的标签重叠也是高频问题,用ggrepel的geom_text_repel可以自动避让:
library(ggrepel) geom_text_repel( data = sites_projected, aes(x = st_coordinates(sites_projected)[,1], y = st_coordinates(sites_projected)[,2], label = site_name), family = "STSong", size = 3, max.overlaps = 20 )记住不要把标签直接放在原始data.frame上,因为sf对象在投影后坐标会变,最好从变换后的对象里提取坐标,否则标签会漂移。
4.3 把业务数据映射到地图坐标
地图本身只是载体,真正有价值的是叠加在空间位置上的业务指标。以生态学中常用的α多样性为例,计算香农多样性指数是vegan包的主场:
library(vegan) shannon <- diversity(species_matrix, index = "shannon") diversity_df <- data.frame(site = rownames(species_matrix), shannon = shannon)拿到指数后,把它和经纬度表merge到一起,再转sf投影,然后画成大小渐变或颜色渐变的点。我常用的映射方式是点大小为样本量、颜色连续映射α多样性指数,这样一张图能同时读出来三个信息维度。
如果是转录组测序这类生信场景,方法也是一样的。基因表达量矩阵要先做标准化才能跨样本比较,很多同学拿到FPKM就直接画图,这是不对的。标准做法是先把FPKM换算成TPM,公式如下:
fpkm_to_tpm <- function(fpkm) { tpm <- fpkm / sum(fpkm) * 1e6 return(tpm) } tpm_matrix <- apply(fpkm_matrix, 2, fpkm_to_tpm)这个换算的口径是简化版,因为严谨的FPKM转TPM还会考虑转录本长度,但实际生信分析中多数场景直接按列归一化到百万分之一就能满足比较需求。换算之后按样品分组求均值,再把均值关联到空间坐标上,配合南美地图,就能直观展示不同生态区域或采样点的表达量空间分布。
说句题外话,类似“生物学年龄”这样的生信指标也是同一个套路:先算指标,再经纬度映射,最后叠加地图。地图的本质就是坐标轴,其他统计量和它无关,但一旦和坐标对齐,就产生了空间故事性。
5. 艺术的本质:配色、布局与主题定制
5.1 配色逻辑:从离散到连续,如何选得高级
画地图最容易暴露审美短板的就是配色。默认的rainbow七彩色基本是灾难现场,饱和度过高、颜色之间没有层次过渡,打印出来更是惨不忍睹。我自从开始画地图之后,默认优先使用viridis家族的配色,它对色觉障碍人群友好,而且渐变非常有质感:
scale_fill_viridis_c(option = "magma", direction = -1)离散型分类变量的配色推荐用RColorBrewer的Set2、Dark2或Accent,这些配色饱和度适中,多个图例并排时不打架。如果你要表达的是“空气污染程度由低到高”这类顺序变量,用单色渐变通常是更高级的选择,比如从浅黄色到暗红色:
scale_fill_gradientn(colors = c("#fef0d9", "#fcbba1", "#fc9272", "#de2d26", "#a50f15"))我更偏好在项目里用scico包里的berlin、lajolla等地质期刊常用配色,这些配色是为数值可视化专门设计的,比RColorBrewer的连续色板更稳定。
关于“r语言层次分割分析方法贡献率”这种描述性统计诉求,如果贡献率要落到地图上,最稳妥的方式是分面多张子图或在地图旁配堆积条形图图例。把贡献率本身作为地图填充色很容易误导读者,因为它本质上是比例概念,不是空间梯度概念。这个我试过,视觉效果会很乱。
5.2 主题与版面:地图也要讲究呼吸感
很多人画地图的毛病是元素堆满全图,既不放比例尺也不管图例位置,出图后又觉得“哪里不对”。地图的版面设计其实讲究呼吸感,该留白的地方必须留白。
我推荐使用theme_void()作为基础主题,然后手动补足必要的图形元素:
library(ggspatial) p_final <- p_base + annotation_scale(location = "bl", width_hint = 0.3) + annotation_north_arrow(location = "tl", which_north = "true", style = north_arrow_fancy_orienteering()) + theme_void() + theme( legend.position = "bottom", legend.title = element_text(family = "STSong", size = 10), legend.text = element_text(family = "STSong", size = 8), plot.margin = margin(10, 10, 10, 10) )比例尺放在左下角,指北针放左上角,图例统一到下方并加标题。这套模板我用了好几年,基本适配所有区域地图项目。如果你的地图用于论文投稿,还要在caption里加入数据来源,用labs(caption = "Data: Natural Earth")即可。字体要格外注意,中文字体在R里的默认family很容易变成豆腐块,常见做法是把字体显式指定为STSong、SimHei或STKaiti。
5.3 组装图与动画联动:把SARIMA预测结果画到地图上
地图不只可以画静态的一个时间点,还可以承载时间序列预测结果。这里我以“空气污染r语言”为背景展开:假设你有南美洲各城市过去三年的PM2.5月均浓度,想用SARIMA模型预测未来一年,并把这12个月的预测结果展现在地图上。
按城市为单位建立ts对象,按月周期做季节差分:
library(forecast) fit <- auto.arima(ts(ts_data$pm25, frequency = 12), seasonal = TRUE, stepwise = FALSE) forecast_12 <- forecast(fit, h = 12)$mean给每个城市都跑一遍模型,然后把预测均值按月份堆成长表:
monthly_df <- do.call(rbind, lapply(cities, function(city){ data.frame(city = city, month = 1:12, pred = get_forecast(city)) }))把monthly_df和城市经纬度合并,再转sf,用facet_wrap(~month)就得到12个月的小倍数地图:
ggplot(sa_aea) + geom_sf(fill = "#f5f5f5", color = "#aaaaaa") + geom_sf(data = pred_sf, aes(color = pred), size = 2) + facet_wrap(~month, ncol = 4) + scale_color_viridis_c(option = "turbo") + theme_void()这组图放在一起,南美洲各城市空气污染的季节性一下子就能看出来了:哪些城市冬季爆表、哪些季节性不明显,一目了然。SARIMA模型不是万能药,自动定阶出来的模型偶尔会有离谱的预测波动,我建议在出图前把预测值可视化检查一遍,把增长率超过500%的异常点纠偏后再进入地图流程。
如果你想要动画效果,可以改用gganimate把月份映射到frame变量,输出gif动画,那视觉冲击力会更强。不过做动画前先确认静态图足够精美,动画只是锦上添花。
6. 常见问题和排查技巧实录
6.1 sf对象Geometry报错和Invalid geometry
使用sf进行空间操作时最常见的报错是“Geometry is invalid”或“TopologyException: Input geom 0 is invalid”。根源通常是底图数据本身存在自相交或重复顶点,这在Natural Earth的低分辨率数据里偶尔会出现。排查方法:
st_is_valid(sa_aea)如果有伪造的要素,用st_make_valid修复:
sa_valid <- st_make_valid(sa_aea)这个函数会把自相交的多边形拆成合法几何,代价是可能让属性表多出几行。修复之后再去做裁剪、相交和点重合判断,才能得到稳定结果。处理高精度的国境线数据时务必记得这一步,我见过不止一次因为没修复导致全省的点都标不上的情况。
6.2 地图上的中文标签乱码
在Windows上直接用ggplot2输出PDF地图,中文标签经常显示成乱码或者空白。解决思路是把绘图设备换成Cairo,或显式指定中文字体:
grDevices::cairo_pdf("map.pdf", width = 8, height = 10) ggplot(...) + theme(text = element_text(family = "STSong")) dev.off()导出PNG时用png(type = "cairo", family = "STSong")也能解决问题。一个小经验:如果你在RStudio里预览图形是正常中文,但ggsave导出后乱码,八九不离十是设备不支持,统一在ggsave里加device = cairo_pdf或type = "cairo"就能修复。
6.3 nc文件读取遇到的时间维与经纬度顺序
处理.nc文件时,维度顺序可能是(time, latitude, longitude),也可能是(longitude, latitude, time),直接用ncvar_get得到的数组维度要和维度变量对齐,否则提取出来的数据完全是乱的。建议:先用dim和dimnames检查,然后用aperm调整:
temp <- ncvar_get(nc, "temperature") # dims: time x lat x lon temp_perm <- aperm(temp, c(3, 2, 1)) # lat x lon x time这样每一层就对应一段时间。时间维对应的变量名常叫time或valid_time,其单位可能不是年月日而是“hours since 1900-01-01”,需要用as.POSIXct换算成年份。换算之后要把网格数据展平成data.frame,再裁剪到南美范围并用ggplot的geom_contour或geom_raster画等值面。这个过程代码量不大,但一旦经纬度顺序搞错,画出来的图会像雪花一样混乱。
6.4 FPKM转TPM、merge与坐标NA的处理
生信数据合并到地图坐标时我踩过最狠的坑是merge之后坐标列出现NA。原因是merge的by参数设错,匹配键的类型不统一,比如左边是字符右边是因子,merge就会静默产生大量不匹配记录。解决方式是在merge前统一类型:
meta$site <- as.character(meta$site) diversity_df$site <- as.character(diversity_df$site) merged <- merge(meta, diversity_df, by = "site", all.x = TRUE)合并后立刻执行complete.cases检查,把缺失坐标的样本过滤掉。不要在地图绘制的时候才发现有NA,那时候整个图层会崩溃。还有一个细节:如果你的坐标里有重复采样点,多个属性值对应同一经纬度,地图上会重叠成小黑圈。建议先按采样点分组汇总:group_by(site) %>% summarise(value = mean(value)),画出平均值并在地图里标注样本量。
以上这些坑每一个我都实际踩过,排查过程快则半小时,慢则一晚上。把它们整理成下面这个速查表,可以帮你快速定位:
| 症状 | 大概率原因 | 处理方式 |
|---|---|---|
| 地图变形严重 | 投影选择错误 | 改用Albers等面积投影 |
| 点标记全部不显示 | 点对象CRS与底图不一致 | st_transform统一CRS |
| 中文标签乱码 | 绘图设备缺少中文字体 | 使用cairo设备并指定family |
| 报错Geometry invalid | 源数据拓扑错误 | st_make_valid |
| 网格数据位置对不上 | nc维度顺序搞错 | aperm调整维度 |
| merge后坐标NA | 匹配键类型不一致 | as.character统一 |
| 预测值画出来离谱 | SARIMA自动定阶过于激进 | 检查并过滤异常预测 |
结尾:几句实际的心里话
画了半年的南美洲地图,我体会最深的一点是:地图的美感,本质上来自投影、配色和留白之间的平衡。R语言给了我们一套强大的工具链,但真正的“艺术”不在函数里,而在于你对数据空间结构的理解——你知道什么时候该选等面积投影,什么时候该用单色渐变而不是彩虹色。最后分享一个小技巧:每次画南美洲底图之前,我会先用st_bbox把几个关键城市(比如利马、巴西利亚、布宜诺斯艾利斯)的坐标打印出来,确认它们是否落在边界框内。这个小步骤能提前发现90%的数据对齐问题,比你画完图再回头检查坐标轴省事得多。希望这篇文章能帮你少走一些弯路,画出自己满意的南美洲地图。