quadmesh icon indicating copy to clipboard operation
quadmesh copied to clipboard

coords for quadmesh()

Open mdsumner opened this issue 5 years ago • 2 comments

Standard like the argument to mesh_plot(), currently can work around with workflow like

## lon, lat, and z raster pseudocode
#r <- raster(file)
#lon <- raster(file, varname = "lon")
#lat <- raster(file, varname = "lat")
#qm <- quadmesh::quadmesh(r)

qm$vb[1, ] <- quadmesh::quadmesh(lon1)$vb[3, ]
qm$vb[2, ] <- quadmesh::quadmesh(lat)$vb[3, ]

with coords we could do

qm <- quadmesh(r, coords = brick(lon, lat))

but also think about corner/centre meshes, and also see #33

mdsumner avatar Mar 13 '20 20:03 mdsumner

The other connection is is angstroms::romsmap() which uses the exact same idea:

https://gist.github.com/mdsumner/61be4632bf5641dc31653d9d7cecdc01#gistcomment-3212031

mdsumner avatar Mar 14 '20 11:03 mdsumner

Should we even allow coords for non-raster input?

library(quadmesh)
library(raster)
#> Loading required package: sp
topo_qm <- quadmesh(etopo)
#> NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
#> NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
#> NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
#> NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
topo_cds <- brick(setValues(etopo, xFromCell(etopo, 1:ncell(etopo))), 
                  setValues(etopo, yFromCell(etopo, 1:ncell(etopo))))
#> NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
#> NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
#> NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files

mp <- maps::map(plot=FALSE)
mesh_plot(topo_qm)

mesh_plot(topo_qm, coords = topo_cds)
points(cbind(mp$x, mp$y), pch = ".")


prj <- "+proj=laea +ellps=sphere"
cds2 <- setValues(topo_cds, reproj::reproj(values(topo_cds), prj))
#> Warning in reproj.matrix(values(topo_cds), prj): 'source' projection not included, but looks like longitude/latitude values:
#>    using '+proj=longlat +datum=WGS84 +no_defs'
mesh_plot(topo_qm, coords = cds2)
points(reproj::reproj(cbind(mp$x, mp$y), prj), pch = ".")
#> Warning in reproj.matrix(cbind(mp$x, mp$y), prj): 'source' projection not included, but looks like longitude/latitude values:
#>    using '+proj=longlat +datum=WGS84 +no_defs'

Created on 2020-03-20 by the reprex package (v0.3.0)

mdsumner avatar Mar 19 '20 20:03 mdsumner