I am attempting to plot a world map of aerosol height data pulled from the Sentinel 5 Precursor/TROPOMI level 2 data from the data hub available here on dropbox.
I am specifically looking at Hawaii (5-30 degree N, 130-180 degree W) and am attempting to plot the aerosol height which is in netcdf4 (.nc) format.
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)
Resulting image 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
Additionally, when I try:
r <- brick(ncname, var = "PRODUCT/aerosol_mid_height")
plot(r)
plot(wrld_simpl, add = TRUE)
The resulting image is still distorted 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
I believe there is a disconnect between the projections of wrld_simpl and test. The extent class for test also seems peculiar. Does anyone know why this problem is occurring? Thank you in advance.
Update
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))
When attempting to plot, I receive this error:
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
Why am I unable to plot the raster?
The standard way to read from a ncdf file with raster is:
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)
If that does not look good, what does show(r) show? And please make a small example file available.
Looking at the file I see that it is does not follow that CF standard that raster (or terra/GDAL) expects. And it may need transformation of some sort. The coordinates are in (PRODUCT/ground_pixel, PRODUCT/scanline). That is, they refer to pixels, not to coordinates, and I cannot see where they are, of they are evenly spaced etc. The file also says that the lat/lon bounds are global. One could hope for the best and do the below but I would not trust this without verification.
extent(r) <- c(-180,180,-90,90)
To see all the metadata you can do
print(r)