suppressPackageStartupMessages(library(stars))
s5p = system.file("sentinel5p/S5P_NRTI_L2__NO2____20180717T120113_20180717T120613_03932_01_010002_20180717T125231.nc", package = "starsdata")
if (s5p != "") {
lat_ds = paste0("HDF5:\"", s5p, "\"://PRODUCT/latitude")
lon_ds = paste0("HDF5:\"", s5p, "\"://PRODUCT/longitude")
nit_ds = paste0("HDF5:\"", s5p, "\"://PRODUCT/SUPPORT_DATA/DETAILED_RESULTS/nitrogendioxide_summed_total_column")
lat = read_stars(lat_ds)
lon = read_stars(lon_ds)
nit = read_stars(nit_ds)
nit[[1]][nit[[1]] > 9e+36] = NA
ll = setNames(c(lon, lat), c("x", "y"))
nit.c0 = st_as_stars(nit, curvilinear = ll)
# more direct method:
nit.c = read_stars(s5p, sub = "//PRODUCT/SUPPORT_DATA/DETAILED_RESULTS/nitrogendioxide_summed_total_column",
curvilinear = c("//PRODUCT/longitude", "//PRODUCT/latitude"), driver = NULL)
if (inherits(nit.c[[1]], "units"))
nit.c[[1]] = units::drop_units(nit.c[[1]])
nit.c[[1]][nit.c[[1]] > 9e+36] = NA
all.equal(nit.c0, nit.c)
st_crs(nit.c) = 4326
print(nit.c)
if (FALSE && capabilities()["png"]) {
png("nit1.png", 800, 800)
plot(nit.c, breaks = "equal", reset = FALSE, axes = TRUE, as_points = TRUE, pch = 16)
maps::map('world', add = TRUE, col = 'red')
dev.off()
png("nit2.png", 800, 800)
plot(nit.c, breaks = "equal", reset = FALSE, axes = TRUE, as_points = FALSE, border = NA)
maps::map('world', add = TRUE, col = 'red')
dev.off()
nit.c = stars:::st_downsample(nit.c, 8)
print(nit.c)
png("nit3.png", 800, 800)
plot(nit.c, breaks = "equal", reset = FALSE, axes = TRUE, as_points = TRUE, pch = 16)
maps::map('world', add = TRUE, col = 'red')
dev.off()
png("nit4.png", 800, 800)
plot(nit.c, breaks = "equal", reset = FALSE, axes = TRUE, as_points = FALSE, border = NA)
maps::map('world', add = TRUE, col = 'red')
dev.off()
}
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.