Please advise me regarding overlay of an LCC-projected map on `raster` data plotted by `rasterVis::layerplot`. What I mean, why I ask: As discussed in detail @ http://stackoverflow.com/questions/14865507/how-to-display-a-projected-map-on-an-rlatticelayerplot (more detail than is feasible for an email post) I'm trying to overlay LCC-projected data from a netCDF file with an LCC-projected map. The netCDF uses the IOAPI convention (described in link), and also manages to hit yet again the `raster` filename-extension problem. (Dunno why, but that's twice in one week for me--bad karma?) The data's horizontal extents are CONUS, so I'm looking for an LCC CONUS map. In one example (in the link above), I use CRAN package=M3 to get (basically) a matrix of coordinates map.lines <- M3::get.map.lines.M3.proj( file=lcc.file, database="state", units="m") and feed that to old-style graphics. But I'd much prefer to make this work like the other example, which resembles other code I have that is using `rasterVis::layerplot`. Unfortunately that other code is using unprojected/lon-lat data, and I must also handle projected (probably all LCC) data. TIA, Tom Roche <Tom_Roche at pobox.com>
how to display a projected map on an rasterVis::layerplot?
4 messages · Tom Roche, Oscar Perpiñan
Hello, Yesterday I read your question in stackoverflow. I tried to display the ozone file with the state boundaries overlaid without success. In my opinion, the problem is that the extent of the RasterLayer as read by raster is quite different from the state boundaries. If I am not wrong, this is the same problem you asked and answered in another stackoverflow question: http://stackoverflow.com/questions/13005181/how-to-set-the-display-bounds-of-a-rasterplot I download the file from https://bitbucket.org/tlroche/mozart-global-to-aqmeii-na/downloads/ozone_lcc.nc and then: library(raster) library(maps) library(maptools) lcc.crs <- "+proj=lcc +lat_1=33 +lat_2=45 +lat_0=40 +lon_0=-97 +a=6370000 +b=6370000" o3.raster <- raster('ozone_lcc.nc', varname="O3", crs=lcc.crs, band=1)
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=longlat +datum=WGS84 data source : /home/oscar/temp/ozone_lcc.nc names : O3 z-value : 1 zvar : O3 level : 1 o3.raster state.map <- map(database="state", plot=FALSE, projection='lambert', parameters=c(40, 33)) state.map.shp <- map2SpatialLines(state.map, proj4string=CRS(lcc.crs))
state.map.shp
class : SpatialLines nfeatures : 169 extent : -0.2226344, 0.2112617, -0.9123794, -0.6458026 (xmin, xmax, ymin, ymax) coord. ref. : +proj=lcc +lat_1=33 +lat_2=45 +lat_0=40 +lon_0=-97 +a=6370000 +b=6370000 Best, Oscar. Grupo de Sistemas Fotovoltaicos (IES-UPM) Dpto. Ingenier?a El?ctrica (EUITI-UPM) URL: http://procomun.wordpress.com Twitter: @oscarperpinan 2013/2/15 Tom Roche <Tom_Roche at pobox.com>:
Please advise me regarding overlay of an LCC-projected map on `raster` data plotted by `rasterVis::layerplot`. What I mean, why I ask: As discussed in detail @ http://stackoverflow.com/questions/14865507/how-to-display-a-projected-map-on-an-rlatticelayerplot (more detail than is feasible for an email post) I'm trying to overlay LCC-projected data from a netCDF file with an LCC-projected map. The netCDF uses the IOAPI convention (described in link), and also manages to hit yet again the `raster` filename-extension problem. (Dunno why, but that's twice in one week for me--bad karma?) The data's horizontal extents are CONUS, so I'm looking for an LCC CONUS map. In one example (in the link above), I use CRAN package=M3 to get (basically) a matrix of coordinates map.lines <- M3::get.map.lines.M3.proj( file=lcc.file, database="state", units="m") and feed that to old-style graphics. But I'd much prefer to make this work like the other example, which resembles other code I have that is using `rasterVis::layerplot`. Unfortunately that other code is using unprojected/lon-lat data, and I must also handle projected (probably all LCC) data. TIA, Tom Roche <Tom_Roche at pobox.com>
_______________________________________________ R-sig-Geo mailing list R-sig-Geo at r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
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 at 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=6370000 +b=6370000 # data source : .../ozone_lcc.nc # names : O3 # z-value : 1 # zvar : O3 # level : 1
# > summary(coordinates(o3.raster)) # x y # 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 :168 NA'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_Roche at pobox.com>
https://stat.ethz.ch/pipermail/r-sig-geo/2013-February/017531.html
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?
One way to do this, though ugly, is the current answer @ http://stackoverflow.com/questions/14865507/how-to-display-a-projected-map-on-an-rlatticelayerplot I.e., "fix" the coordinates in the `state.map` obtained from `maps::map`, using the linear IOAPI transform, and process that into a SpatialLines as previously--see code and plot @ link above. Is there a cleaner way to do this?