在 R 中从 netcdf 创建栅格的最准确方法是什么？
What is the most accurate way of creating a raster from netcdf in R?
多年来我一直在处理 netCDF 数据。 netCDF 用于空气污染物数据，纬度和经度作为 单独的 变量提供，而不是原始网格的一部分。
LINK 迄今为止：Sample Netcdf
这些 netCDF 文件提供 2 级二氧化氮数据，它们是从 NASA Earthdata 门户网站下载的。卫星为Sentinel-5P，仪器为TROPOMI
所以在处理这个数据的时候，你要为NO2，纬度和经度创建变量。我正在尝试创建栅格图层，然后将它们保存为 GeoTIFF 文件以供我研究。
这里的问题与我不知道如何最好地创建这些栅格有关。纬度和经度在整个数据集中的分布不均，我还没有找到准确创建这些图像的方法。我使用 netCDF 文件提供的行数和列数创建了一个模型网格。在变量列表中，这称为扫描线 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, "DETAILED_RESULTS/nitrogendioxide_total_column",
"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,lonvec,latvec)
colnames(dfsat) <- c('z', 'x', 'y')
colnames(dfqa) <- c('z', 'x', 'y')
df <- rbind(df, dfsat)
dfqa <- rbind(df,dfqa)
rm(lat,lon,no2tc,qa,latvec,lonvec,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[, 2:3], r, df[, 1], fun = mean)
crs(xx) <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"
crs(qa_raster) <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,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, "PRODUCT/nitrogendioxide_tropospheric_column","_FillValue")
my_unit = ncdf4::ncatt_get(nc, "PRODUCT/nitrogendioxide_tropospheric_column","units")
my_product_name = ncdf4::ncatt_get(nc, "PRODUCT/nitrogendioxide_tropospheric_column", "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,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, lat_mid), fun = geosphere::distHaversine)
extent_distance_horizontal <- one_degree_horizontal_distance *
horizontal_distance
} else {
extent_distance_horizontal <-
geosphere::distm(c(raster::extent(pts)[1],
lat_mid),
c(raster::extent(pts)[2], lat_mid),
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)
多年来我一直在处理 netCDF 数据。 netCDF 用于空气污染物数据，纬度和经度作为 单独的 变量提供，而不是原始网格的一部分。
LINK 迄今为止：Sample Netcdf
这些 netCDF 文件提供 2 级二氧化氮数据，它们是从 NASA Earthdata 门户网站下载的。卫星为Sentinel-5P，仪器为TROPOMI
所以在处理这个数据的时候，你要为NO2，纬度和经度创建变量。我正在尝试创建栅格图层，然后将它们保存为 GeoTIFF 文件以供我研究。
这里的问题与我不知道如何最好地创建这些栅格有关。纬度和经度在整个数据集中的分布不均，我还没有找到准确创建这些图像的方法。我使用 netCDF 文件提供的行数和列数创建了一个模型网格。在变量列表中，这称为扫描线 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, "DETAILED_RESULTS/nitrogendioxide_total_column",
"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,lonvec,latvec)
colnames(dfsat) <- c('z', 'x', 'y')
colnames(dfqa) <- c('z', 'x', 'y')
df <- rbind(df, dfsat)
dfqa <- rbind(df,dfqa)
rm(lat,lon,no2tc,qa,latvec,lonvec,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[, 2:3], r, df[, 1], fun = mean)
crs(xx) <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"
crs(qa_raster) <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,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, "PRODUCT/nitrogendioxide_tropospheric_column","_FillValue")
my_unit = ncdf4::ncatt_get(nc, "PRODUCT/nitrogendioxide_tropospheric_column","units")
my_product_name = ncdf4::ncatt_get(nc, "PRODUCT/nitrogendioxide_tropospheric_column", "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,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, lat_mid), fun = geosphere::distHaversine)
extent_distance_horizontal <- one_degree_horizontal_distance *
horizontal_distance
} else {
extent_distance_horizontal <-
geosphere::distm(c(raster::extent(pts)[1],
lat_mid),
c(raster::extent(pts)[2], lat_mid),
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)