是否有更好的方法从R中的Sentinel 5p产品加载hdf5文件并生成光栅文件.tif？

问题描述

我正在尝试从有效载荷Sentinel 5p（ESA）的卫星导入R HDF5文件。目标是创建一个GeoTIFF文件，以供Python中的另一个应用程序使用。 到目前为止，我已经开发了以下代码，但是我不确定光栅的分辨率。有更好的方法吗？ 谢谢！

$ mypy foo.py 
Success: no issues found in 1 source file

library(RNetCDF)
library(rgdal)
library(ncdf4)
library(raster)
library(stringr)
library(gdalUtils)
library(rgdal)
library(gstat)

### define a rectangular extent I want to use to crop my HDF5 file
e <- extent(c(11.92586,15.65291,35.49295,38.81766))
plot(e)

# make a spatial polygon from the extent
p <- as(e,"SpatialPolygons")
plot(p)
proj4string(p) <- CRS("+proj=longlat +datum=WGS84")


# save shp file for the rectangular DOMAIN
setwd("...dir/sicily_shapefile")
shapefile(p,"rectangular_domain.shp",overwrite=TRUE)

# reload and plot rectangular domani domain
shp_rect <- readOGR(dsn = dir,layer = "rectangular_domain")
# ----- Transform to EPSG 4326 - WGS84 
shp_rect <- spTransform(shp_rect,CRS("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0"))
shp_rect@data$name <- 1:nrow(shp_rect)
plot(shp_rect)

## setup directory containing an HDF5 Sentinel 5p file (for NO2)
setwd("....dir/sentinel5p")


## load all HDF5 file
filenames <- "S5P_NRTI_L2__NO2____20200901T122909_20200901T123409_14955_01_010302_20200901T132207.nc"

sentinel_file <- nc_open(filenames)
name_vari <- names(sentinel_file $var)


## get longitude,latitude and NO2 data (these are arrays)
LON <- ncvar_get(sentinel_file,"PRODUCT/longitude")
LAT <- ncvar_get(sentinel_file,"PRODUCT/latitude")
NO2 <- ncvar_get(sentinel_file,"PRODUCT/nitrogendioxide_tropospheric_column")

## transform arrays into vectors
LON <- as.vector(LON)
LAT <- as.vector(LAT)
NO2 <- as.vector((NO2))

## make a spatial matrix
LAT <- matrix(LAT,nrow=450,ncol=372,byrow=TRUE)
LON <- matrix(LON,byrow=TRUE)
NO2 <- matrix(NO2,byrow=TRUE)   
    
xmn= min(LON)
xmx=max(LON)
ymn=min(LAT)
ymx=max(LAT)
    

LAT <- c(LAT)
LON <- c(LON)
M <- c(NO2)

### make an unique vector with lat,lon and data
data <- cbind(LON,LAT,M)
data <- as.data.frame(data)
colnames(data) <- c("Lon","Lat","M")

xmn= min(data$Lon)
xmx=max(data$Lon)
ymn=min(data$Lat)
ymx=max(data$Lat)
 
## rasterize the vector (lon,lat,M)
## setup the resolution of the raster at ~ 10km
resl_ras <- 0.1
colnames(data) <- c('x','y','z')

x.range <- as.numeric(c(xmn,xmx))  
y.range <- as.numeric(c(ymn,ymx)) 

## make the base grid
grd <- expand.grid(x = seq(from = x.range[1],to = x.range[2],by = resl_ras),y = seq(from = y.range[1]*0.95,to = y.range[2]*0.95,by = resl_ras))  
    
grd_1<- dplyr::filter(grd,grd$x == grd$x[1])
grd_2<- dplyr::filter(grd,grd$y == grd$y[1])


r <- raster(xmn=min(data$x),xmx=max(data$x),ymn=min(data$y),ymx=max(data$y),ncol=nrow(grd_2),nrow= nrow(grd_1))

r <- rasterize(data[,1:2],r,data[,3],fun=mean)
projection(r) <- CRS("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0")
plot(r)

### crop all over the rectangular domani
r <- crop(r,extent(shp_rect))
r <- mask(r,shp_rect) 

plot(r)
## save raster as .tif file (GeoTIFF)
writeRaster(r,"NO2_Sentinel5p.tif",options= "INTERLEAVE=BAND",overwrite=T)

#### visualize raster in a Leafleft map  

library(leaflet)

map <- leaflet() %>%
  addTiles() %>%
  addRasterImage(r,opacity = 0.5,group = "SEVIRI")
map

解决方法

暂无找到可以解决该程序问题的有效方法，小编努力寻找整理中！

如果你已经找到好的解决方法，欢迎将解决方案带上本链接一起发送给小编。

小编邮箱:dio#foxmail.com (将#修改为@）