minimal L1B in R

[mdsumner]

Posted here just to archive it, I'm cleaning up

## download a MODIS file
f <- 
"ftp://ladsweb.nascom.nasa.gov/allData/6/MOD021KM/2015/114/MOD021KM.A2015114.0200.006.2015114135358.hdf"
library(raster)
library(rgdal)
download.file(f, file.path("/rdsi/PRIVATE/dev/modis_ses/test", basename(f)), mode = "wb")
 f <- file.path("/rdsi/PRIVATE/dev/modis_ses/test", basename(f))
 
 ## get the subdataset names (we want lon, lat, and emissivity)
 tx <- system(sprintf("gdalinfo %s", f), intern = TRUE)
 patt <- "^  SUBDATASET_.*_NAME="
 gnames <- gsub(patt, "", grep(patt, tx, value = TRUE))
 lon <- raster(gnames[21])
 lat <- raster(gnames[20])
 em <- raster(gnames[24])
 
## set up a local projection
proj <- sprintf("+proj=laea +lat_0=%i +lon_0=%i", as.integer(mean(values(lat))), as.integer(mean(values(lon))))
 
 ## resample the lon/lat mesh to the data dimensions (should do this in planar coords probably)
## and this could be sped up by doing the resampling once only for lon and apply the index to lat
mlon <- resample(lon, setExtent(em, extent(lon)), method = "bilinear")
mlat <- resample(lat, setExtent(em, extent(lat)), method = "bilinear")
xy <- project(cbind(values(mlon), values(mlat)), proj)
scl <- function(x) (x - min(x, na.rm = TRUE))/diff(range(x, na.rm = TRUE))
plot(xy, pch = ".", col = grey(seq(0,1, length = 255))[scl(values(em)) * 254 + 1], asp = 1, axes = FALSE, xlab = "", ylab = "")

library(rworldxtra)
data(countriesHigh)
w <- spTransform(subset(countriesHigh, SOVEREIGNT == "Antarctica"), CRS(proj))
plot(w, add = TRUE, border = "lightgrey")
## dummy object for plotting the graticule
dummy <- spTransform(rasterToPoints(raster(extent(range(values(lon)), range(values(lat))),  nrow = 50, ncol = 50, crs = "+proj=longlat +ellps=WGS84"), spatial = TRUE), CRS(proj))
op <- par(xpd = NA); llgridlines(dummy, ndiscr = 80, side = "EN"); par(op)