【R代码分析】甲烷排放绘制代码-参考论文

论文-J2021-【墨西哥】Unravelling a large methane emission discrepancy in Mexico using satellite observations-Remote Sensing of Environment

相应绘图代码下载-Shen_2020_Mexico（完整R绘图代码）

Figure1：

处理和可视化墨西哥东部区域的甲烷（CH₄）排放和卫星观测数据（TROPOMI XCH₄）

墨西哥东部区域的甲烷（CH₄）排放如下图所示：

R代码解释

1. 初始化与加载依赖包

rm(list=ls())
library(fields)
library(maps)
library(ncdf4)
library(abind)
library(geosphere)
library(mapdata)
library(Hmisc)

• rm(list=ls())：清除当前工作环境中的所有变量。
• 加载多个处理地理数据、绘图、netCDF 文件等的 R 包。

---

2. 设置工作目录与加载自定义函数

setwd("~/Documents/CH4")
source('Function/get_geo.R')
source('Function/get_met.R')

• 设置工作目录。
• 加载两个自定义函数脚本，可能包含获取地理或气象数据的函数。

---

3. 自定义函数：plot.box

plot.box=function(loc,lwd=lwd){arrows(x0=loc[1],y0=loc[3],x1=loc[2],code=0,col=1,lwd=lwd)arrows(x0=loc[1],y0=loc[4],x1=loc[2],code=0,col=1,lwd=lwd)arrows(x0=loc[1],y0=loc[3],y1=loc[4],code=0,col=1,lwd=lwd)arrows(x0=loc[2],y0=loc[3],y1=loc[4],code=0,col=1,lwd=lwd)
}

• 这个函数绘制矩形框，用于地图上标出感兴趣区域。
• loc 是一个4元素向量：c(xmin, xmax, ymin, ymax)。

---

4. 自定义函数：plot.field

这是一个复杂而强大的绘图函数，用于绘制二维空间数据（如遥感图像、气体浓度场）。

该函数：

• 使用 image.plot（来自 fields 包）绘图。
• 支持不同类型数据（正负值、百分比、绝对值等）。
• 支持自定义颜色、图例、坐标范围等。
• 支持太平洋中心地图或默认地图。

主要功能包括：

• 数据预处理（翻转纬度、设置色标范围等）
• 自动或自定义颜色尺度
• 添加地图边框（map()）

---

5. 加载排放数据并修改某区域

ss=load("/Users/lu/Documents/CH4_Mexico/Data/Emis_Tia/emis_per_area.Rdata")
ind11=(lon.in>=-93.0 & lon.in <=-91.7)
ind21=(lat.in>=18.6 & lat.in <=19.8)
emis_per_area[ind11,ind21]=emis_per_area[ind11,ind21]*0.1

• 加载每平方公里 CH₄ 排放数据。
• 将一个特定区域的排放量缩小为原来的10%（可能是数据校正或假设）。

---

6. 加载卫星数据：TROPOMI XCH₄

ss=load("/Users/lu/Documents/CH4_Mexico/Figures/Figure 1/TROPOMI_XCH4_01x01.Rdata")

• 加载卫星观测数据，包括 XCH₄（柱浓度甲烷），分辨率 0.1° × 0.1°。

---

7. 绘图：两幅图组成的图像

pdf(file="CH4_posteriri_Mexico.pdf",width=6.5,height=2.5)
par(mfrow=c(1,2))

• 输出 PDF 文件，包含 2 幅子图（1 行 2 列）。

(a) TROPOMI XCH₄ 图像

ap= avg_corrected; ap[count<=1]=NA
plot.field(ap[ind1,ind2], ...)
title("(a) TROPOMI XCH4, elevation corrected", ...)

• 使用 avg_corrected 数据（即海拔校正后的 XCH₄ 平均值）。
• 将观测次数小于等于 1 的像元设为 NA。
• 使用 plot.field 绘图。
• 添加多个矩形框（plot.box）来标出研究区域。

(b) 排放数据图像

ap= emis_per_area
ap[ap==0]=NA
plot.field(log10(ap[ind1,ind2]), ...)
title("(b) Prior oil/gas emissions in 2015", ...)

• 加载排放数据，取对数（log10）以便更好地视觉呈现。
• 零值设为 NA。
• 添加标注（R1-R7）对应不同的感兴趣区域。

---

8. 绘制一个感兴趣区域（如 R5）的放大图

pdf(file="R5.pdf",width=4.5,height=3)
xlim=c(-94.5,-91.9); ylim=c(17.3,18.6)
...
points(lonlat[,1],lonlat[,2],cex=0.5)

• 画出研究区域 R5 的放大图。
• points() 添加具体的设施或测站点位置（经纬度）。
• 继续使用 plot.field 显示 TROPOMI 数据。

---

Figure2：

参考