AustralianAntarcticDivision / angstroms

Tools for model grid data (with raster in R)
http://australianantarcticdivision.github.io/angstroms/
4 stars 1 forks source link

better polar example #21

Open mdsumner opened 4 years ago

mdsumner commented 4 years ago
library(angstroms)
#> Loading required package: raster
#> Loading required package: sp
library(anglr)

## 1.) remap to a regular lonlat grid
x <- romsdata3d("ai.1979.nc")
coords <- brick(romsdata2d("lon0.nc"), 
                romsdata2d("lat0.nc"))

## points we want
r <- raster(extent(0, 360, -85, -50), 
            res = c(1.5, 0.5), crs = "+proj=longlat +datum=WGS84")
ll <- coordinates(r)
## map them to the native grid
xy <- romsmap(ll, coords = coords)

plot(x[[1]])  ## first time step
contour(coords[[1]], add = TRUE, col = "dodgerblue") ## longitude lines
contour(coords[[2]], add = TRUE, col = "firebrick") ## latitude lines
points(xy, pch = ".")


## now extract those point values
r1 <- setValues(r, extract(x[[1]], xy, method = "bilinear"))
plot(r1)
maps::map("world2", add = TRUE)


## 2.) leverage mesh_plot to plot a polar projection
## I was a bit surprised this works, gives me something to go on
prj <- "+proj=laea +lat_0=-90 +lon_0=147 +datum=WGS84"
wmap <- sp::spTransform(anglr::simpleworld, 
                        prj)
## remap coords for plot
## (could use a similar trick as above to create this as a regular grid)
coords_xy <- setValues(coords, reproj::reproj(values(coords), 
                                              source = "+proj=longlat +datum=WGS84", 
                                              target = prj)[,1:2])

mesh_plot(x[[1]], coords = coords_xy)
plot(wmap, add = TRUE)

Created on 2020-04-29 by the reprex package (v0.3.0)