如何在世界地图上叠加 netcdf4 栅格图像？
How to overlay a netcdf4 raster image on world map?
我正在尝试绘制从 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 的预测之间存在脱节。 extent class for test 似乎也很奇怪。有谁知道为什么会出现这个问题？提前谢谢你。
更新
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
为什么我无法绘制光栅？
从带有光栅的 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) 中。也就是说，它们指的是像素，而不是坐标，我看不到它们在哪里，它们是均匀分布的等等。该文件还说 lat/lon 边界是全局的。人们可以抱有最好的希望并执行以下操作，但未经验证我不会相信这一点。
extent(r) <- c(-180,180,-90,90)
要查看您可以执行的所有元数据
print(r)
我正在尝试绘制从 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 的预测之间存在脱节。 extent class for test 似乎也很奇怪。有谁知道为什么会出现这个问题？提前谢谢你。
更新
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
为什么我无法绘制光栅？
从带有光栅的 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) 中。也就是说，它们指的是像素，而不是坐标，我看不到它们在哪里，它们是均匀分布的等等。该文件还说 lat/lon 边界是全局的。人们可以抱有最好的希望并执行以下操作，但未经验证我不会相信这一点。
extent(r) <- c(-180,180,-90,90)
要查看您可以执行的所有元数据
print(r)