如何在世界地图上叠加 netcdf4 光栅图像？

Question

我正在尝试绘制从 Sentinel 5 Precursor/TROPOMI 2 级数据中提取的气溶胶高度数据的世界地图，该数据来自数据中心 here 在 dropbox 上。

我特别关注夏威夷（北纬 5-30 度，西经 130-180 度），并试图绘制 netcdf4 (.nc) 格式的气溶胶高度。

library(ncdf4)
library(raster)
library(ggplot2)
library(sp)
library(ggmap)
library(rworldmap)
library(maps)
library(usmap)
library(sf)
library(rasterVis)
library(lattice)

ncname <- "~/Desktop/Summer 2020/Tropomi/Aerosol Height/S5P_RPRO_L2__AER_LH_20180715T201241_20180715T215609_03908_01_010301_20190530T173144.nc"
nc <- nc_open(ncname)

lon <- ncvar_get(nc, varid = "PRODUCT/longitude")
lat <- ncvar_get(nc, varid = "PRODUCT/latitude")
time <- ncvar_get(nc, varid = "PRODUCT/time_utc")
ah <- ncvar_get(nc, varid = "PRODUCT/aerosol_mid_height")

test <- raster("S5P_RPRO_L2__AER_LH_20180715T201241_20180715T215609_03908_01_010301_20190530T173144.nc", 
               varname = "PRODUCT/aerosol_mid_height", 
               xmn=min(lon), xmx=max(lon), ymn=min(lat), ymx=max(lat), 
               crs=CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs+ towgs=0,0,0"))

plot(test)
plot(wrld_simpl, add = TRUE)

结果图片distorted image 1

> wrld_simpl
class       : SpatialPolygonsDataFrame 
features    : 246 
extent      : -180, 180, -90, 83.57027  (xmin, xmax, ymin, ymax)
crs         :  +proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0 
variables   : 11
names       : FIPS, ISO2, ISO3,  UN,           NAME,    AREA,    POP2005, REGION, SUBREGION,      LON,     LAT 
min values  :     ,   AD,  ABW,   4, Aaland Islands,       0,          0,      0,         0, -178.131, -80.446 
max values  :   ZI,   ZW,  ZWE, 894,       Zimbabwe, 1638094, 1312978855,    150,       155,  179.219,   78.83 

> test
class      : RasterLayer 
dimensions : 3245, 448, 1453760  (nrow, ncol, ncell)
resolution : 1, 1  (x, y)
**extent     : -0.5, 447.5, -0.5, 3244.5  (xmin, xmax, ymin, ymax)**
crs        : +proj=longlat +ellps=WGS84 +datum=WGS84 +towgs84=0,0,0 
source     : /Users/mattcohen/Desktop/Summer 2020/Tropomi/Aerosol Height/S5P_RPRO_L2__AER_LH_20180715T201241_20180715T215609_03908_01_010301_20190530T173144.nc 
names      : Height.at.center.of.aerosol.layer.relative.to.geoid 
z-value    : 269308800 
zvar       : PRODUCT/aerosol_mid_height 

另外，当我尝试时：

r <- brick(ncname, var = "PRODUCT/aerosol_mid_height")
plot(r)
plot(wrld_simpl, add = TRUE)

生成的图像仍然失真 distorted image 2

show(r)
class      : RasterBrick 
dimensions : 3245, 448, 1453760, 1  (nrow, ncol, ncell, nlayers)
resolution : 1, 1  (x, y)
extent     : -0.5, 447.5, -0.5, 3244.5  (xmin, xmax, ymin, ymax)
crs        : NA 
source     : /Users/mattcohen/Desktop/Summer 2020/Tropomi/Aerosol Height/S5P_RPRO_L2__AER_LH_20180715T201241_20180715T215609_03908_01_010301_20190530T173144.nc 
names      : X269308800 
PRODUCT/time (seconds since 2010-01-01 00:00:00): 269308800 
varname    : PRODUCT/aerosol_mid_height 

我相信wrld_simpl 和test 的预测之间存在脱节。 test 的 extent 类似乎也很奇特。有谁知道为什么会出现这个问题？提前谢谢你。

更新

s1 <- data.frame(as.vector(lon), as.vector(lat), as.vector(ah))
crsLatLon <- "+proj=longlat +datum=WGS84"
ex <- extent(c(-180,180,-90,90))
pmraster <- raster(ncol=360*10, nrow=180*10, crs=crsLatLon,ext=ex)
pmraster <- rasterize(s1[,1:2], pmraster, s1[,3], fun=mean, na.rm=T)

exHI <- extent(c(-180,-140,10,30))
levelplot(crop(pmraster,exHI))

尝试绘图时，我收到此错误：

Error: $ operator is invalid for atomic vectors
In addition: Warning messages:
1: In min(x) : no non-missing arguments to min; returning Inf
2: In max(x) : no non-missing arguments to max; returning -Inf
3: In min(x) : no non-missing arguments to min; returning Inf
4: In max(x) : no non-missing arguments to max; returning -Inf

为什么我无法绘制栅格？

Answer 1

从带有光栅的 ncdf 文件中读取的标准方法是：

library(raster)
ncname <- "~/Desktop/Summer 2020/Tropomi/Aerosol Height/S5P_RPRO_L2__AER_LH_20180715T201241_20180715T215609_03908_01_010301_20190530T173144.nc"
r <- brick(ncname, var = "PRODUCT/aerosol_mid_height")

plot(r)
plot(wrld_simpl, add = TRUE)

如果这看起来不太好，show(r) 会显示什么？请提供一个小示例文件。

查看文件，我发现它不符合光栅（或 terra/GDAL）所期望的 CF 标准。它可能需要某种形式的转换。坐标位于 (PRODUCT/ground_pixel, PRODUCT/scanline)。也就是说，它们指的是像素，而不是坐标，我看不到它们在哪里，它们是均匀分布的等等。该文件还说纬度/经度边界是全局的。人们可以抱有最好的希望并执行以下操作，但未经验证我不会相信这一点。

 extent(r) <- c(-180,180,-90,90)

要查看您可以执行的所有元数据

 print(r)