如何解决从 R 中的 netcdf 创建栅格的最准确方法是什么？
我处理了多年的 netCDF 数据。 netCDF 用于空气污染物数据，纬度和经度作为单独变量提供，而不是作为原始网格的一部分。
这些 netCDF 文件提供 2 级二氧化氮数据，它们是从 NASA Earthdata 门户下载的。卫星为Sentinel-5P，仪器为TROPOMI。
因此，在处理这些数据时，您必须为 NO2、纬度和经度创建变量。我正在尝试创建栅格图层，然后将它们保存为 GeoTIFF 文件以供我研究。
这里的问题与我不知道如何最好地创建这些栅格有关。整个数据集中的纬度和经度不是等距的，我还没有想出一种方法来准确地创建这些图像。我使用 netCDF 文件提供的行数和列数创建了一个模型网格。在变量列表中，这称为 scanline 和 ground_pixel，但是当我绘制它时，最终图像中的单元格看起来不正确。
这是我上传数据的方式：
## Open the netcdf
ncname <- no2files$filename[m]
ncfname <- paste(ncname,sep = "")
nc <- nc_open(ncfname)
## Get the necessary variables.
no2tc <-ncvar_get(nc,"PRODUCT/nitrogendioxide_tropospheric_column")
lat <- ncvar_get(nc,"PRODUCT/latitude")
lon <- ncvar_get(nc,"PRODUCT/longitude")
qa <- ncvar_get(nc,"PRODUCT/qa_value")
fillvalue = ncatt_get(nc,"DETAILED_RESULTS/nitrogendioxide_total_column","_FillValue")
mfactor <- ncatt_get(nc,"multiplication_factor_to_convert_to_molecules_percm2")
fillvalue_qa = ncatt_get(nc,"PRODUCT/qa_value","_FillValue")
no2tc[no2tc == fillvalue$value] <- NA
no2tc <- no2tc * mfactor$value
qa[qa == fillvalue_qa$value] <- NA
nc_close(nc)
# rm(ncfname)
no2vec <- as.vector(no2tc)
latvec <- as.vector(lat)
lonvec <- as.vector(lon)
qavec <- as.vector(qa)
dfsat <- data.frame(no2vec,lonvec,latvec)
dfqa <- data.frame(qavec,latvec)
colnames(dfsat) <- c('z','x','y')
colnames(dfqa) <- c('z','y')
df <- rbind(df,dfsat)
dfqa <- rbind(df,dfqa)
rm(lat,lon,no2tc,qa,latvec,no2vec,qavec)
这是我目前创建栅格的方式：
## Create the raster. The ncol = 3245 and Now = 450 are from the scanline and ground_pixel variables.
e <- extent(-180,180,-90,90)
r <- raster(e,ncol = 3245,nrow = 450)
xx <- rasterize(df[,2:3],r,df[,1],fun = mean)
qa_raster <- rasterize(dfqa[,fun = mean)
crs(xx) <- "+proj=longlat +datum=wgs84 +no_defs +ellps=wgs84 +towgs84=0,0"
crs(qa_raster) <- "+proj=longlat +datum=wgs84 +no_defs +ellps=wgs84 +towgs84=0,0"
## Crop and plot the raster
## change shapefile coordinate system
# border <- spTransform(ontario,crs(xx))
aoi <- spTransform(ontario_buffer,crs(xx))
## Mask values with qa < 0.5 (this is the recommended value)
xx[qa_raster < 0.5 & xx < 0] <- NA
## This is the final plot
plot_tif <- crop(xx,extent(aoi))
### Use this if you want to view the plot.
mask_tif <- mask(plot_tif,aoi)
# plot(mask_tif)
# final <- plot(border,add=TRUE)
## Plot the raster
filename <- paste(i,".tif",sep="")
writeraster(mask_tif,filename = filename,"GTiff",overwrite=TRUE)
最终结果如下：
然后我尝试了我在网上找到的另一种方法，但是您必须设置分辨率。我可以这样做，但我只想按原样绘制单元格，无需任何修改。
ncfname <- "S5P_OFFL_L2__NO2____20200107T173517_20200107T191647_11582_01_010302_20200109T103930.nc"
nc <- ncdf4::nc_open(ncfname)
mfactor = ncdf4::ncatt_get(nc,"PRODUCT/nitrogendioxide_tropospheric_column","multiplication_factor_to_convert_to_molecules_percm2")
fillvalue = ncdf4::ncatt_get(nc,"_FillValue")
my_unit = ncdf4::ncatt_get(nc,"units")
my_product_name = ncdf4::ncatt_get(nc,"long_name")
mfactor <- mfactor$value
fillvalue <- fillvalue$value
vals <- ncdf4::ncvar_get(nc,"PRODUCT/nitrogendioxide_tropospheric_column")
lat <- ncdf4::ncvar_get(nc,"PRODUCT/latitude")
lon <- ncdf4::ncvar_get(nc,"PRODUCT/longitude")
vals[vals == fillvalue] <- NA
vals_df = NULL
vals_df <- rbind(vals_df,data.frame(lat = as.vector(lat),lon = as.vector(lon),vals = as.vector(vals)))
pts <- vals_df
sp::coordinates(pts) <- ~lon + lat
my_projection <- "+proj=longlat +datum=wgs84 +no_defs +ellps=wgs84 +towgs84=0,0"
sp::proj4string(pts) <- sp::CRS(my_projection)
my_aoi <- ontario
crs_test <- raster::compareCRS(pts,my_aoi)
my_aoi <- sp::spTransform(my_aoi,CRS = as.character(raster::crs(pts)))
p <- methods::as(raster::extent(my_aoi),"Spatialpolygons")
sp::proj4string(p) <- sp::CRS(my_projection)
pts <- raster::crop(pts,p)
extent_distance_vertical <- geosphere::distm(c(raster::extent(pts)[1],raster::extent(pts)[3]),c(raster::extent(pts)[1],raster::extent(pts)[4]),fun = geosphere::disthaversine)
vertical_mid_distance <- (raster::extent(pts)[4] - raster::extent(pts)[3])/2
lat_mid <- raster::extent(pts)[3] + vertical_mid_distance
horizontal_distance <- raster::extent(pts)[2] - raster::extent(pts)[1]
if (horizontal_distance > 180) {
one_degree_horizontal_distance <- geosphere::distm(c(1,lat_mid),c(2,fun = geosphere::disthaversine)
extent_distance_horizontal <- one_degree_horizontal_distance *
horizontal_distance
} else {
extent_distance_horizontal <-
geosphere::distm(c(raster::extent(pts)[1],c(raster::extent(pts)[2],fun = geosphere::disthaversine)
}
my_res <- 20000
ncol_rast <- as.integer(extent_distance_horizontal/my_res)
nrow_rast <- as.integer(extent_distance_vertical/my_res)
print(paste0("Create raster file from points"))
rast <- raster::raster(nrows = nrow_rast,ncols = ncol_rast,crs = as.character(raster::crs(pts)),ext = raster::extent(pts),vals = NULL)
final <- raster::rasterize(pts,rast,pts$vals,fun = mean)
final <- raster::mask(final,my_aoi)
sp::plot(final)
如何准确地创建这些栅格图层？谢谢！
解决方法
使用示例文件
f <- "S5P_OFFL_L2__NO2____20200107T173517_20200107T191647_11582_01_010302_20200109T103930.nc"
你可以做到
library(terra)
r <- rast(f,paste0("/PRODUCT/",c("longitude","latitude","nitrogendioxide_tropospheric_column")))
r
#class : SpatRaster
#dimensions : 4172,450,3 (nrow,ncol,nlyr)
#resolution : 1,1 (x,y)
#extent : -0.5,449.5,-0.5,4171.5 (xmin,xmax,ymin,ymax)
#coord. ref. :
#sources : longitude
# latitude
# nitrogendioxide_tropospheric_column
#varnames : longitude (pixel center longitude)
# latitude (pixel center latitude)
# nitrogendioxide_tropospheric_column (Tropospheric vertical column of nitrogen dioxide)
#names : longitude,latitude,nitrogendi~ric_column
#unit : degrees_east,degrees_north,mol m-2
#time : 2020-01-07
plot(r,nr=1)
该图说明数据不是按常规栅格数据组织的（事实上，纬度会有 N-S 梯度，经度会有 E-W 梯度）。另见plot(r$longitude,r$latitude)。您可以将它们视为积分：
dp <- as.data.frame(r)
p <- vect(dp,geom=c("longitude","latitude"))
绘图需要一段时间，因为有 > 150 万个点，所以我取了一个样本
plot(p[sample(nrow(p),10000)],"nitrogendioxide_tropospheric_column")
如果您希望将数据组织为常规栅格，可以使用 rasterize
x <- rast(res=1/6)
x <- rasterize(p,x,"nitrogendioxide_tropospheric_column",fun=mean)
plot(x > 0)
你可以像这样得到安大略
can <- vect(raster::getData("GADM",country="CAN",level=1))
ontario <- can[can$NAME_1=="Ontario",]
x <- crop(x,ontario)
x <- mask(x,ontario)
plot(x)
版权声明：本文内容由互联网用户自发贡献，该文观点与技术仅代表作者本人。本站仅提供信息存储空间服务，不拥有所有权，不承担相关法律责任。如发现本站有涉嫌侵权/违法违规的内容， 请发送邮件至 dio@foxmail.com 举报，一经查实，本站将立刻删除。
