summary: How to convert from a simple (x,y) matrix (as produced by
M3::get.map.lines.M3.proj) to a SpatialLines, or to coordinates
consumable by latticeExtra::layer?
details:
Perhaps I should have asked the question above first, instead of
putting it at the end of
http://stackoverflow.com/questions/14865507/how-to-display-a-projected-map-on-an-rlatticelayerplot
where it is obviously overlooked :-( The problem is also explained
there:
http://stackoverflow.com/questions/14865507/how-to-display-a-projected-map-on-an-rlatticelayerplot
(formatted for mail)
## Read in variable with daily ozone.
o3.raster - raster::raster(
x=lcc.file, varname=O3, crs=lcc.crs, level=1)
o3.raster@crs - lcc.crs # why does the above assignment not take?
# start debugging
o3.raster
summary(coordinates(o3.raster)) # thanks, Felix Andrews!
M3::get.grid.info.M3(lcc.file)
# end debugging
# o3.raster
# class : RasterLayer
# band: 1
# dimensions : 112, 148, 16576 (nrow, ncol, ncell)
# resolution : 1, 1 (x, y)
# extent : 0.5, 148.5, 0.5, 112.5 (xmin, xmax, ymin, ymax)
# coord. ref. :
+proj=lcc +lat_1=33 +lat_2=45 +lat_0=40 +lon_0=-97
+a=637 +b=637
# data source : .../ozone_lcc.nc
# names : O3
# z-value : 1
# zvar: O3
# level : 1
# summary(coordinates(o3.raster))
#xy
# Min. : 1.00 Min. : 1.00
# 1st Qu.: 37.75 1st Qu.: 28.75
# Median : 74.50 Median : 56.50
# Mean : 74.50 Mean : 56.50
# 3rd Qu.:111.25 3rd Qu.: 84.25
# Max. :148.00 Max. :112.00
# M3::get.grid.info.M3(lcc.file)
# $x.orig
# [1] -2736000
# $y.orig
# [1] -2088000
# $x.cell.width
# [1] 36000
# $y.cell.width
# [1] 36000
# $hz.units
# [1] m
# $ncols
# [1] 148
# $nrows
# [1] 112
# $nlays
# [1] 1
# The grid (`coordinates(o3.raster)`) is integers, because this
# data is stored using the IOAPI convention. IOAPI makes the grid
# Cartesian by using an (approximately) equiareal projection (LCC)
# and abstracting grid positioning to metadata stored in netCDF
# global attributes. Some of this spatial metadata can be accessed
# by `M3::get.grid.info.M3` (above), and some can be accessed via
# the coordinate reference system (`CRS`, see `lcc.proj4` above)
However I can
## Get state boundaries in projection units.
map.lines - M3::get.map.lines.M3.proj(
file=lcc.file, database=state, units=m)
# start debugging
class(map.lines)
summary(map.lines)
summary(map.lines$coords)
# end debugging
# class(map.lines)
# [1] list
# summary(map.lines)
#Length Class Mode
# coords 23374 -none- numeric
# units 1 -none- character
# summary(map.lines$coords)
#x y
# Min. :-2272238 Min. :-1567156
# 1st Qu.: 94244 1st Qu.: -673953
# Median : 913204 Median : -26948
# Mean : 689997 Mean : -84644
# 3rd Qu.: 1744969 3rd Qu.: 524531
# Max. : 2322260 Max. : 1265778
# NA's :168NA's :168
Note how these grid coordinates derive linearly from the IOAPI
metadata above, and to the raster bounds (farther above): e.g.,
map.lines.XY.min.raw -
c(-2272238,-1567156) # from summary(map.lines$coords)
x.orig - -2736000
x.cell.width - 36000
y.orig - -2088000
y.cell.width - 36000
map.lines.XY.min.IOAPI - c(
(map.lines.XY.min.raw[1] - x.orig)/x.cell.width,
(map.lines.XY.min.raw[2] - y.orig)/y.cell.width
)
map.lines.XY.min.IOAPI
# [1] 12.88228 14.46789
Looking @
http://i.stack.imgur.com/eijZI.png (example 2)
eyeball the intersection of the latitude of south tip of Texas with
the westernmost longitude of California. Then look at the
corresponding location in
http://i.stack.imgur.com/z7OA0.png (example 1)
you will see that corresponding example1 coordinates are approximately
(12.88228,14.46789) as computed above. Hardly a proof :-) but the this
certainly shows that the transform from `map.lines`-coordinates to
IOAPI-coordinates is straightforward. So what I need to know, in order
to feed output from M3::get.map.lines.M3.proj(...) to
latticeExtra::layer(...), is how to set the coordinates of a
SpatialLines as produced by, e.g.,
state.map.shp -
maptools::map2SpatialLines(state.map, proj4string=lcc.crs)
or
latticeExtra::layer(
sp::sp.lines(state.map.shp, lwd=0.8, col='darkgray'))
or otherwise how to set the (x,y) coordinates being read by
latticeExtra::layer(...)
TIA, Tom Roche tom_ro...@pobox.com
___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
