curvilinear support (rasters with coordinate arrays)

[mdsumner]

From issue these files contain "lonlat" arrays which are coordinates stored as data (instead of a regular transform). These can be plotted with the anglr package function mesh_plot(var, coords = brick(lon, lat)) where 'var' is one of the raster layers, and lon,lat the first (2,1) arrays shown in the figure above.

I've long wanted something in raster/terra to "know about" this kind of curvilinear coordinates, and I put all my helpers into angstroms and the anglr package (with long plans to somehow formalize those ideas better). The crux perspective is that of a "quad mesh", like the mesh3d type in rgl - which stores each cell corner uniquely with a mesh index.

Finally, what is seen as 'curvilinear' is simply sometimes a map projection with untransformed lonlat arrays, but that's not always the case.

[mdsumner]

Here's a perspective using the anglr::mesh_plot() function, no warping required or intense GDAL grid registration, just plot as a mesh.

  library(raster)
#> Loading required package: sp

the_ras <- paste(tempdir(), "/vwnd.10m.2015.nc", sep = "")

utils::download.file("ftp://ftp.cdc.noaa.gov/Datasets/NARR/Dailies/monolevel/vwnd.10m.2015.nc",
                     the_ras, 
                     method = "auto", 
                     mode = "wb", 
                     cacheOK = TRUE, 
                     extra = getOption("download.file.extra"))

lonlat <- brick(raster(the_ras, varname = "lon"), 
                raster(the_ras, varname = "lat"))
#> Loading required namespace: ncdf4

r <- raster(the_ras)

library(anglr)  ## not on CRAN atm remotes::install_github("hypertidy/anglr")

## looks a bit crazy
mesh_plot(r, coords = lonlat)

## because the longitudes are crazy
plot(lonlat)

## choose a projection that's more suited to the region
mesh_plot(r, coords = lonlat, crs = "+proj=laea +lat_0=90")

## and "fix" the longitude values in how they wrap

lon <- lonlat[[1]]
msk <- lon > 50
lon[msk] <- lon[msk] - 360
#plot(lon)

prj <- "+proj=laea +lat_0=90"
mesh_plot(r, coords = brick(lon, lonlat[[2]]), crs = prj, asp = 1)
mm <- maps::map(plot = FALSE)
points(rgdal::project(cbind(mm$x, mm$y), prj), , pch = ".")

## what is the projection probably? (something something North America, im guessing wildly)
alb <- "+proj=aea +lat_1=20 +lat_2=60 +lat_0=40 +lon_0=-96 +x_0=0 +y_0=0 +ellps=GRS80 +datum=NAD83 +units=m +no_defs"
mesh_plot(r, coords = brick(lon, lonlat[[2]]), crs = alb, asp = 1)
points(rgdal::project(cbind(mm$x, mm$y), alb), , pch = ".")

^{Created on 2020-12-01 by the reprex package (v0.3.0)}

[mdsumner]

(I see that projection is properly LCC, but that's not really the point here - this data set can just be treated as a proper raster, but the general case still holds)