summary: I've tried several ways to draw wrld_simpl over global data
plotted with rasterVis::levelplot, but keep encountering the same
problem:
- only the eastern hemisphere of the map plots
- the eastern hemisphere of the map overlays the western hemisphere of
the data
How to fix? or otherwise provide geographical context?
details:
(I'll attempt to produce a "complete, self-contained example" after I
make the code publicly viewable, but for now, I'm hoping the problem is
sufficiently simple to diagnose.)
I have atmospheric data in a netCDF file specifying gas concentrations
over a 3D space with dimensions longitude, latitude, and (vertical)
level, such that
* the horizontal extents are global
* the vertical coordinates are pressures (specifically hybrid
sigma-pressure)
package=raster provides very usable API for loading the data into a
RasterBrick, and package=rasterVis makes it relatively easy to levelplot
the data in a single "atmospheric column" with the highest pressures at
the bottom. However I also want to provide some geographic context
(e.g., to overlay the data with national boundaries), and have failed to
date.
I have done this with maptools::wrld_simpl successfully in other
contexts, i.e.,
* loading the data with package=ncdf4 (which is considerably more work)
* converting the data to a data.frame
* lattice::levelplot-ing the data
* lattice::xyplot-ing the map
But when I do
data(wrld_simpl)
global.raster <- brick(in.fp, varname=data.var.name)
layers.n <- nbands(global.raster)
global.data.df <- as.data.frame(global.raster)
# some twiddling of names(global.data.df) omitted
global.map <- map('world', plot=FALSE)
global.map.df <- data.frame(lon=global.map$x, lat=global.map$y)
pdf(file=pdf.fp, width=5, height=100)
# using rasterVis::levelplot, not lattice::levelplot
levelplot(global.raster,
margin=FALSE,
names.attr=names(global.data.df),
layout=c(1,layers.n), # all layers in one "atmospheric column"
) + xyplot(lat ~ lon, global.map.df, type='l', lty=1, lwd=1, col='black')
dev.off()
the data plots apparently correctly, but the map plots incorrectly, as
shown by this .png of the top level:
https://docs.google.com/open?id=0BzDAFHgIxRzKOVdlSDFMR29BdWM
As noted above, the map's eastern hemisphere overlays the data's western
hemisphere, and nothing overlays the data's eastern hemisphere.
I tried several variations on the above, but either got the same
results, or outright error. (Unfortunately the error messages are often
printed on the levelplot, but the messages are truncated (so they can't
be understood), and it's annoying to "fail slow" (after waiting for a
long plot to finish).) So after some research I sought to emulate the
means recommended @
http://rastervis.r-forge.r-project.org/#levelplot
data(wrld_simpl)
global.raster <- brick(in.fp, varname=data.var.name)
layers.n <- nbands(global.raster)
global.data.df <- as.data.frame(global.raster)
# some twiddling of names(global.data.df) omitted
global.map.proj <- CRS('+proj=latlon +ellps=WGS84')
global.map <- map('world', plot=FALSE)
global.map.shp <- map2SpatialLines(global.map, proj4string=global.map.proj)
pdf(file=pdf.fp, width=5, height=100)
# using rasterVis::levelplot, not lattice::levelplot
levelplot(global.raster,
margin=FALSE,
names.attr=names(global.data.df),
layout=c(1,layers.n), # all layers in one "atmospheric column"
) + layer(sp.lines(global.map.shp, lwd=0.8, col='darkgray'))
dev.off()
This also fails, and very similarly (though the line colors differ). See
this .png:
https://docs.google.com/open?id=0BzDAFHgIxRzKYzd6cVB2TDE4RE0
I'm guessing I can either do something simple to wrld_simpl's longitude
values, or Something Completely Different that would Just Work Better,
but don't know. Your assistance is appreciated.
TIA, Tom Roche <Tom_Roche at pobox.com>
problems drawing a world map on a levelplot
11 messages · Michael Sumner, Tom Roche, Matthew Landis +1 more
I happened to write this yesterday, which switches the parts of
wrld_simpl from < 0 longitude to the east (Pacific view). I haven't
investigated what you need in detail, so just ignore if I'm off track.
It's not careful about data or anything, since all I needed was the
raw geometry. It might be of use. "xrange" here can be anything from
[-180, 360] though it's really assuming you don't traverse more than
360 in total.
It would be nice to carefully clean up a copy of wrld_simpl like this
and make a wrld_simpl version like the maps package has.
Cheers, Mike.
xrange <- c(140, 250)
##xrange <- c(0, 360)
yrange <- c(-60, 20)
require(maptools)
require(raster)
require(rgeos)
ext <- extent(xrange[1], xrange[2], yrange[1], yrange[2])
data(wrld_simpl, package = "maptools")
if (xrange[2] > 180.0) {
eastrange <- c(xrange[1], 180.0)
westrange <- c(-180.0, xrange[2] - 360)
ew <- extent(westrange[1], westrange[2], yrange[1], yrange[2])
ee <- extent(eastrange[1], eastrange[2], yrange[1], yrange[2])
geome <- as(ee, "SpatialPolygons")
geomw <- as(ew, "SpatialPolygons")
## why does this get dropped?
proj4string(geome) <- CRS(proj4string(wrld_simpl))
proj4string(geomw) <- CRS(proj4string(wrld_simpl))
worlde <- gIntersection(wrld_simpl, geome)
worldw <- gIntersection(wrld_simpl, geomw)
worldw <- elide(worldw, shift = c(360.0, 0.0))
proj4string(worldw) <- CRS(proj4string(wrld_simpl))
dfe <- data.frame(x = seq_len(length(worlde at polygons)))
dfw <- data.frame(x = seq_len(length(worldw at polygons)))
rownames(dfw) <- as.character(as.numeric(rownames(dfe)) + nrow(dfe))
worlde <- SpatialPolygonsDataFrame(worlde, dfe, match.ID = FALSE)
worldw <- SpatialPolygonsDataFrame(worldw, dfw, match.ID = FALSE)
worldw <- spChFIDs(worldw, rownames(dfw))
world <- spRbind(worlde, worldw)
} else {
geom <- as(ext, "SpatialPolygons")
proj4string(geom) <- CRS(proj4string(wrld_simpl))
## TODO: make this spdf if needed
world <- gIntersection(wrld_simpl, geom)
}
plot(world)
On Tue, Dec 18, 2012 at 8:15 PM, Tom Roche <Tom_Roche at pobox.com> wrote:
summary: I've tried several ways to draw wrld_simpl over global data
plotted with rasterVis::levelplot, but keep encountering the same
problem:
- only the eastern hemisphere of the map plots
- the eastern hemisphere of the map overlays the western hemisphere of
the data
How to fix? or otherwise provide geographical context?
details:
(I'll attempt to produce a "complete, self-contained example" after I
make the code publicly viewable, but for now, I'm hoping the problem is
sufficiently simple to diagnose.)
I have atmospheric data in a netCDF file specifying gas concentrations
over a 3D space with dimensions longitude, latitude, and (vertical)
level, such that
* the horizontal extents are global
* the vertical coordinates are pressures (specifically hybrid
sigma-pressure)
package=raster provides very usable API for loading the data into a
RasterBrick, and package=rasterVis makes it relatively easy to levelplot
the data in a single "atmospheric column" with the highest pressures at
the bottom. However I also want to provide some geographic context
(e.g., to overlay the data with national boundaries), and have failed to
date.
I have done this with maptools::wrld_simpl successfully in other
contexts, i.e.,
* loading the data with package=ncdf4 (which is considerably more work)
* converting the data to a data.frame
* lattice::levelplot-ing the data
* lattice::xyplot-ing the map
But when I do
data(wrld_simpl)
global.raster <- brick(in.fp, varname=data.var.name)
layers.n <- nbands(global.raster)
global.data.df <- as.data.frame(global.raster)
# some twiddling of names(global.data.df) omitted
global.map <- map('world', plot=FALSE)
global.map.df <- data.frame(lon=global.map$x, lat=global.map$y)
pdf(file=pdf.fp, width=5, height=100)
# using rasterVis::levelplot, not lattice::levelplot
levelplot(global.raster,
margin=FALSE,
names.attr=names(global.data.df),
layout=c(1,layers.n), # all layers in one "atmospheric column"
) + xyplot(lat ~ lon, global.map.df, type='l', lty=1, lwd=1, col='black')
dev.off()
the data plots apparently correctly, but the map plots incorrectly, as
shown by this .png of the top level:
https://docs.google.com/open?id=0BzDAFHgIxRzKOVdlSDFMR29BdWM
As noted above, the map's eastern hemisphere overlays the data's western
hemisphere, and nothing overlays the data's eastern hemisphere.
I tried several variations on the above, but either got the same
results, or outright error. (Unfortunately the error messages are often
printed on the levelplot, but the messages are truncated (so they can't
be understood), and it's annoying to "fail slow" (after waiting for a
long plot to finish).) So after some research I sought to emulate the
means recommended @
http://rastervis.r-forge.r-project.org/#levelplot
data(wrld_simpl)
global.raster <- brick(in.fp, varname=data.var.name)
layers.n <- nbands(global.raster)
global.data.df <- as.data.frame(global.raster)
# some twiddling of names(global.data.df) omitted
global.map.proj <- CRS('+proj=latlon +ellps=WGS84')
global.map <- map('world', plot=FALSE)
global.map.shp <- map2SpatialLines(global.map, proj4string=global.map.proj)
pdf(file=pdf.fp, width=5, height=100)
# using rasterVis::levelplot, not lattice::levelplot
levelplot(global.raster,
margin=FALSE,
names.attr=names(global.data.df),
layout=c(1,layers.n), # all layers in one "atmospheric column"
) + layer(sp.lines(global.map.shp, lwd=0.8, col='darkgray'))
dev.off()
This also fails, and very similarly (though the line colors differ). See
this .png:
https://docs.google.com/open?id=0BzDAFHgIxRzKYzd6cVB2TDE4RE0
I'm guessing I can either do something simple to wrld_simpl's longitude
values, or Something Completely Different that would Just Work Better,
but don't know. Your assistance is appreciated.
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
Michael Sumner Hobart, Australia e-mail: mdsumner at gmail.com
I meant "wrld_simpl2 version like the maps package has".
On Tue, Dec 18, 2012 at 8:54 PM, Michael Sumner <mdsumner at gmail.com> wrote:
I happened to write this yesterday, which switches the parts of
wrld_simpl from < 0 longitude to the east (Pacific view). I haven't
investigated what you need in detail, so just ignore if I'm off track.
It's not careful about data or anything, since all I needed was the
raw geometry. It might be of use. "xrange" here can be anything from
[-180, 360] though it's really assuming you don't traverse more than
360 in total.
It would be nice to carefully clean up a copy of wrld_simpl like this
and make a wrld_simpl version like the maps package has.
Cheers, Mike.
xrange <- c(140, 250)
##xrange <- c(0, 360)
yrange <- c(-60, 20)
require(maptools)
require(raster)
require(rgeos)
ext <- extent(xrange[1], xrange[2], yrange[1], yrange[2])
data(wrld_simpl, package = "maptools")
if (xrange[2] > 180.0) {
eastrange <- c(xrange[1], 180.0)
westrange <- c(-180.0, xrange[2] - 360)
ew <- extent(westrange[1], westrange[2], yrange[1], yrange[2])
ee <- extent(eastrange[1], eastrange[2], yrange[1], yrange[2])
geome <- as(ee, "SpatialPolygons")
geomw <- as(ew, "SpatialPolygons")
## why does this get dropped?
proj4string(geome) <- CRS(proj4string(wrld_simpl))
proj4string(geomw) <- CRS(proj4string(wrld_simpl))
worlde <- gIntersection(wrld_simpl, geome)
worldw <- gIntersection(wrld_simpl, geomw)
worldw <- elide(worldw, shift = c(360.0, 0.0))
proj4string(worldw) <- CRS(proj4string(wrld_simpl))
dfe <- data.frame(x = seq_len(length(worlde at polygons)))
dfw <- data.frame(x = seq_len(length(worldw at polygons)))
rownames(dfw) <- as.character(as.numeric(rownames(dfe)) + nrow(dfe))
worlde <- SpatialPolygonsDataFrame(worlde, dfe, match.ID = FALSE)
worldw <- SpatialPolygonsDataFrame(worldw, dfw, match.ID = FALSE)
worldw <- spChFIDs(worldw, rownames(dfw))
world <- spRbind(worlde, worldw)
} else {
geom <- as(ext, "SpatialPolygons")
proj4string(geom) <- CRS(proj4string(wrld_simpl))
## TODO: make this spdf if needed
world <- gIntersection(wrld_simpl, geom)
}
plot(world)
On Tue, Dec 18, 2012 at 8:15 PM, Tom Roche <Tom_Roche at pobox.com> wrote:
summary: I've tried several ways to draw wrld_simpl over global data
plotted with rasterVis::levelplot, but keep encountering the same
problem:
- only the eastern hemisphere of the map plots
- the eastern hemisphere of the map overlays the western hemisphere of
the data
How to fix? or otherwise provide geographical context?
details:
(I'll attempt to produce a "complete, self-contained example" after I
make the code publicly viewable, but for now, I'm hoping the problem is
sufficiently simple to diagnose.)
I have atmospheric data in a netCDF file specifying gas concentrations
over a 3D space with dimensions longitude, latitude, and (vertical)
level, such that
* the horizontal extents are global
* the vertical coordinates are pressures (specifically hybrid
sigma-pressure)
package=raster provides very usable API for loading the data into a
RasterBrick, and package=rasterVis makes it relatively easy to levelplot
the data in a single "atmospheric column" with the highest pressures at
the bottom. However I also want to provide some geographic context
(e.g., to overlay the data with national boundaries), and have failed to
date.
I have done this with maptools::wrld_simpl successfully in other
contexts, i.e.,
* loading the data with package=ncdf4 (which is considerably more work)
* converting the data to a data.frame
* lattice::levelplot-ing the data
* lattice::xyplot-ing the map
But when I do
data(wrld_simpl)
global.raster <- brick(in.fp, varname=data.var.name)
layers.n <- nbands(global.raster)
global.data.df <- as.data.frame(global.raster)
# some twiddling of names(global.data.df) omitted
global.map <- map('world', plot=FALSE)
global.map.df <- data.frame(lon=global.map$x, lat=global.map$y)
pdf(file=pdf.fp, width=5, height=100)
# using rasterVis::levelplot, not lattice::levelplot
levelplot(global.raster,
margin=FALSE,
names.attr=names(global.data.df),
layout=c(1,layers.n), # all layers in one "atmospheric column"
) + xyplot(lat ~ lon, global.map.df, type='l', lty=1, lwd=1, col='black')
dev.off()
the data plots apparently correctly, but the map plots incorrectly, as
shown by this .png of the top level:
https://docs.google.com/open?id=0BzDAFHgIxRzKOVdlSDFMR29BdWM
As noted above, the map's eastern hemisphere overlays the data's western
hemisphere, and nothing overlays the data's eastern hemisphere.
I tried several variations on the above, but either got the same
results, or outright error. (Unfortunately the error messages are often
printed on the levelplot, but the messages are truncated (so they can't
be understood), and it's annoying to "fail slow" (after waiting for a
long plot to finish).) So after some research I sought to emulate the
means recommended @
http://rastervis.r-forge.r-project.org/#levelplot
data(wrld_simpl)
global.raster <- brick(in.fp, varname=data.var.name)
layers.n <- nbands(global.raster)
global.data.df <- as.data.frame(global.raster)
# some twiddling of names(global.data.df) omitted
global.map.proj <- CRS('+proj=latlon +ellps=WGS84')
global.map <- map('world', plot=FALSE)
global.map.shp <- map2SpatialLines(global.map, proj4string=global.map.proj)
pdf(file=pdf.fp, width=5, height=100)
# using rasterVis::levelplot, not lattice::levelplot
levelplot(global.raster,
margin=FALSE,
names.attr=names(global.data.df),
layout=c(1,layers.n), # all layers in one "atmospheric column"
) + layer(sp.lines(global.map.shp, lwd=0.8, col='darkgray'))
dev.off()
This also fails, and very similarly (though the line colors differ). See
this .png:
https://docs.google.com/open?id=0BzDAFHgIxRzKYzd6cVB2TDE4RE0
I'm guessing I can either do something simple to wrld_simpl's longitude
values, or Something Completely Different that would Just Work Better,
but don't know. Your assistance is appreciated.
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
-- Michael Sumner Hobart, Australia e-mail: mdsumner at gmail.com
Michael Sumner Hobart, Australia e-mail: mdsumner at gmail.com
Here's a crack at a full example, sorry the clip/merge polygon stuff
is a mess still.
library(maptools)
library(raster)
## first, do the jiggery pokery of wrld_simpl to
## convert from Atlantic to Pacific view
xrange <- c(0, 360)
yrange <- c(-90, 90)
require(maptools)
require(raster)
require(rgeos)
ext <- extent(xrange[1], xrange[2], yrange[1], yrange[2])
data(wrld_simpl, package = "maptools")
## this is a bit more general than needed for this example, since I
## wanted arbitrary crops originally
eastrange <- c(xrange[1], 180.0)
westrange <- c(-180.0, xrange[2] - 360)
ew <- extent(westrange[1], westrange[2], yrange[1], yrange[2])
ee <- extent(eastrange[1], eastrange[2], yrange[1], yrange[2])
geome <- as(ee, "SpatialPolygons")
geomw <- as(ew, "SpatialPolygons")
## why does this get dropped?
proj4string(geome) <- CRS(proj4string(wrld_simpl))
proj4string(geomw) <- CRS(proj4string(wrld_simpl))
worlde <- gIntersection(wrld_simpl, geome)
worldw <- gIntersection(wrld_simpl, geomw)
worldw <- elide(worldw, shift = c(360.0, 0.0))
proj4string(worldw) <- CRS(proj4string(wrld_simpl))
dfe <- data.frame(x = seq_len(length(worlde at polygons)))
dfw <- data.frame(x = seq_len(length(worldw at polygons)))
rownames(dfw) <- as.character(as.numeric(rownames(dfe)) + nrow(dfe))
worlde <- SpatialPolygonsDataFrame(worlde, dfe, match.ID = FALSE)
worldw <- SpatialPolygonsDataFrame(worldw, dfw, match.ID = FALSE)
worldw <- spChFIDs(worldw, rownames(dfw))
world <- spRbind(worlde, worldw)
## so, "world" is longitudes [0, 360], which is pretty much essential
## given the input data shown
r <- raster(nrows = 180, ncols = 360, xmn = 0, xmx = 360, ymn = -90,
ymx = 90, crs = "+proj=longlat +ellps=WGS84 +over")
## fake a raster so we have something equivalent to your data (but simple/r)
rast <- rasterize(world, r)
## fill the ocean parts with something . . .
rast[is.na(rast)] <- 1:sum(is.na(getValues(rast)))
library(rasterVis)
levelplot(rast, margin = FALSE) + layer(sp.polygons(world, lwd=0.8,
col='darkgray'))
On Tue, Dec 18, 2012 at 8:55 PM, Michael Sumner <mdsumner at gmail.com> wrote:
I meant "wrld_simpl2 version like the maps package has". On Tue, Dec 18, 2012 at 8:54 PM, Michael Sumner <mdsumner at gmail.com> wrote:
I happened to write this yesterday, which switches the parts of
wrld_simpl from < 0 longitude to the east (Pacific view). I haven't
investigated what you need in detail, so just ignore if I'm off track.
It's not careful about data or anything, since all I needed was the
raw geometry. It might be of use. "xrange" here can be anything from
[-180, 360] though it's really assuming you don't traverse more than
360 in total.
It would be nice to carefully clean up a copy of wrld_simpl like this
and make a wrld_simpl version like the maps package has.
Cheers, Mike.
xrange <- c(140, 250)
##xrange <- c(0, 360)
yrange <- c(-60, 20)
require(maptools)
require(raster)
require(rgeos)
ext <- extent(xrange[1], xrange[2], yrange[1], yrange[2])
data(wrld_simpl, package = "maptools")
if (xrange[2] > 180.0) {
eastrange <- c(xrange[1], 180.0)
westrange <- c(-180.0, xrange[2] - 360)
ew <- extent(westrange[1], westrange[2], yrange[1], yrange[2])
ee <- extent(eastrange[1], eastrange[2], yrange[1], yrange[2])
geome <- as(ee, "SpatialPolygons")
geomw <- as(ew, "SpatialPolygons")
## why does this get dropped?
proj4string(geome) <- CRS(proj4string(wrld_simpl))
proj4string(geomw) <- CRS(proj4string(wrld_simpl))
worlde <- gIntersection(wrld_simpl, geome)
worldw <- gIntersection(wrld_simpl, geomw)
worldw <- elide(worldw, shift = c(360.0, 0.0))
proj4string(worldw) <- CRS(proj4string(wrld_simpl))
dfe <- data.frame(x = seq_len(length(worlde at polygons)))
dfw <- data.frame(x = seq_len(length(worldw at polygons)))
rownames(dfw) <- as.character(as.numeric(rownames(dfe)) + nrow(dfe))
worlde <- SpatialPolygonsDataFrame(worlde, dfe, match.ID = FALSE)
worldw <- SpatialPolygonsDataFrame(worldw, dfw, match.ID = FALSE)
worldw <- spChFIDs(worldw, rownames(dfw))
world <- spRbind(worlde, worldw)
} else {
geom <- as(ext, "SpatialPolygons")
proj4string(geom) <- CRS(proj4string(wrld_simpl))
## TODO: make this spdf if needed
world <- gIntersection(wrld_simpl, geom)
}
plot(world)
On Tue, Dec 18, 2012 at 8:15 PM, Tom Roche <Tom_Roche at pobox.com> wrote:
summary: I've tried several ways to draw wrld_simpl over global data
plotted with rasterVis::levelplot, but keep encountering the same
problem:
- only the eastern hemisphere of the map plots
- the eastern hemisphere of the map overlays the western hemisphere of
the data
How to fix? or otherwise provide geographical context?
details:
(I'll attempt to produce a "complete, self-contained example" after I
make the code publicly viewable, but for now, I'm hoping the problem is
sufficiently simple to diagnose.)
I have atmospheric data in a netCDF file specifying gas concentrations
over a 3D space with dimensions longitude, latitude, and (vertical)
level, such that
* the horizontal extents are global
* the vertical coordinates are pressures (specifically hybrid
sigma-pressure)
package=raster provides very usable API for loading the data into a
RasterBrick, and package=rasterVis makes it relatively easy to levelplot
the data in a single "atmospheric column" with the highest pressures at
the bottom. However I also want to provide some geographic context
(e.g., to overlay the data with national boundaries), and have failed to
date.
I have done this with maptools::wrld_simpl successfully in other
contexts, i.e.,
* loading the data with package=ncdf4 (which is considerably more work)
* converting the data to a data.frame
* lattice::levelplot-ing the data
* lattice::xyplot-ing the map
But when I do
data(wrld_simpl)
global.raster <- brick(in.fp, varname=data.var.name)
layers.n <- nbands(global.raster)
global.data.df <- as.data.frame(global.raster)
# some twiddling of names(global.data.df) omitted
global.map <- map('world', plot=FALSE)
global.map.df <- data.frame(lon=global.map$x, lat=global.map$y)
pdf(file=pdf.fp, width=5, height=100)
# using rasterVis::levelplot, not lattice::levelplot
levelplot(global.raster,
margin=FALSE,
names.attr=names(global.data.df),
layout=c(1,layers.n), # all layers in one "atmospheric column"
) + xyplot(lat ~ lon, global.map.df, type='l', lty=1, lwd=1, col='black')
dev.off()
the data plots apparently correctly, but the map plots incorrectly, as
shown by this .png of the top level:
https://docs.google.com/open?id=0BzDAFHgIxRzKOVdlSDFMR29BdWM
As noted above, the map's eastern hemisphere overlays the data's western
hemisphere, and nothing overlays the data's eastern hemisphere.
I tried several variations on the above, but either got the same
results, or outright error. (Unfortunately the error messages are often
printed on the levelplot, but the messages are truncated (so they can't
be understood), and it's annoying to "fail slow" (after waiting for a
long plot to finish).) So after some research I sought to emulate the
means recommended @
http://rastervis.r-forge.r-project.org/#levelplot
data(wrld_simpl)
global.raster <- brick(in.fp, varname=data.var.name)
layers.n <- nbands(global.raster)
global.data.df <- as.data.frame(global.raster)
# some twiddling of names(global.data.df) omitted
global.map.proj <- CRS('+proj=latlon +ellps=WGS84')
global.map <- map('world', plot=FALSE)
global.map.shp <- map2SpatialLines(global.map, proj4string=global.map.proj)
pdf(file=pdf.fp, width=5, height=100)
# using rasterVis::levelplot, not lattice::levelplot
levelplot(global.raster,
margin=FALSE,
names.attr=names(global.data.df),
layout=c(1,layers.n), # all layers in one "atmospheric column"
) + layer(sp.lines(global.map.shp, lwd=0.8, col='darkgray'))
dev.off()
This also fails, and very similarly (though the line colors differ). See
this .png:
https://docs.google.com/open?id=0BzDAFHgIxRzKYzd6cVB2TDE4RE0
I'm guessing I can either do something simple to wrld_simpl's longitude
values, or Something Completely Different that would Just Work Better,
but don't know. Your assistance is appreciated.
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
-- Michael Sumner Hobart, Australia e-mail: mdsumner at gmail.com
-- Michael Sumner Hobart, Australia e-mail: mdsumner at gmail.com
Michael Sumner Hobart, Australia e-mail: mdsumner at gmail.com
An embedded and charset-unspecified text was scrubbed... Name: not available URL: <https://stat.ethz.ch/pipermail/r-sig-geo/attachments/20121219/91829958/attachment.pl>
https://stat.ethz.ch/pipermail/r-sig-geo/2012-December/017014.html (rearranged)
following your code I think you are not setting the CRS of your RasterBrick.
True, but since it was also global data, brick(...) "did the right thing." I have now set its CRS just to be sure.
Instead of:
global.map.proj <- CRS('+proj=latlon +ellps=WGS84')
you should try:
global.map.proj <- CRS('+proj=longlat +ellps=WGS84')
Thanks again, but please note that I was following (see starred line) http://rastervis.r-forge.r-project.org/#levelplot
let's add the administrative borders. This information is available at the GADM service.
library(maptools)
** proj <- CRS('+proj=latlon +ellps=WGS84')
##Change to your folder
mapaSHP <- readShapeLines('~/Datos/ESP_adm/ESP_adm2.shp', proj4string=proj)
p <- levelplot(SISmm, layers=1, FUN.margin=median)
p + layer(sp.lines(mapaSHP, lwd=0.8, col='darkgray'))
...
Date: 2012-08-10 Author: Oscar Perpi??n Lamigueiro
Perhaps you could fix that? TIA, Tom Roche <Tom_Roche at pobox.com>
summary: How to best/easiest range-shift lon-lat data? details: https://stat.ethz.ch/pipermail/r-sig-geo/2012-December/017011.html
I happened to write [the following code, omitted] yesterday, which switches the parts of wrld_simpl from < 0 longitude to the east (Pacific view).
<self-administered dopeslap/> Of course! the problem is my input! see starred line: R session
global.proj <- CRS('+proj=longlat +ellps=WGS84')
global.raster <- brick(in.fp, varname=data.var.name, crs=global.proj)
global.raster
class : RasterBrick dimensions : 96, 144, 13824, 56 (nrow, ncol, ncell, nlayers) resolution : 2.5, 1.894737 (x, y)
* extent : -1.25, 358.75, -90.94737, 90.94737 (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
which is wrong: - latitudes outside the range [-90, +90]? - (as you note) 0-based longitudes are likely troublesome (I'm pretty sure the other data with which this must integrate is [-180, +180]) - ... and these longitudes aren't even 0-based! I didn't even check :-( since the data's source is much more eminent than a mere student like myself. But now I note $ ncdump -h ./2008N2O_restart.nc | head
netcdf \2008N2O_restart {
dimensions:
lev = 56 ;
lat = 96 ;
lon = 144 ;
ilev = 57 ;
variables:
double N2O(lev, lat, lon) ;
N2O:long_name = "N2O" ;
N2O:units = "ppmV" ;
$ ncdump -v lat ./2008N2O_restart.nc | tail -n 28
lat = -90, -88.1052631578947, -86.2105263157895, -84.3157894736842,
-82.4210526315789, -80.5263157894737, -78.6315789473684,
-76.7368421052632, -74.8421052631579, -72.9473684210526,
-71.0526315789474, -69.1578947368421, -67.2631578947368,
-65.3684210526316, -63.4736842105263, -61.5789473684211,
-59.6842105263158, -57.7894736842105, -55.8947368421053, -54,
-52.1052631578947, -50.2105263157895, -48.3157894736842,
-46.4210526315789, -44.5263157894737, -42.6315789473684,
-40.7368421052632, -38.8421052631579, -36.9473684210526,
-35.0526315789474, -33.1578947368421, -31.2631578947368,
-29.3684210526316, -27.4736842105263, -25.5789473684211,
-23.6842105263158, -21.7894736842105, -19.8947368421053, -18,
-16.1052631578947, -14.2105263157895, -12.3157894736842,
-10.4210526315789, -8.52631578947369, -6.63157894736842,
-4.73684210526316, -2.84210526315789, -0.947368421052634,
0.947368421052634, 2.84210526315788, 4.73684210526315,
6.63157894736842, 8.52631578947369, 10.4210526315789,
12.3157894736842, 14.2105263157895, 16.1052631578947, 18,
19.8947368421053, 21.7894736842105, 23.6842105263158,
25.5789473684211, 27.4736842105263, 29.3684210526316,
31.2631578947368, 33.1578947368421, 35.0526315789474,
36.9473684210526, 38.8421052631579, 40.7368421052632,
42.6315789473684, 44.5263157894737, 46.4210526315789,
48.3157894736842, 50.2105263157895, 52.1052631578947, 54,
55.8947368421053, 57.7894736842105, 59.6842105263158,
61.578947368421, 63.4736842105263, 65.3684210526316,
67.2631578947368, 69.1578947368421, 71.0526315789474,
72.9473684210526, 74.8421052631579, 76.7368421052632,
78.6315789473684, 80.5263157894737, 82.4210526315789,
84.3157894736842, 86.2105263157895, 88.1052631578947, 90 ; }
I'm not sure why raster is reporting y range=[-90.94737, 90.94737], except that 2 * 0.94737 = 1.89474 == the latitude resolution. $ ncdump -v lon ./2008N2O_restart.nc | tail -n 15
lon = 0, 2.5, 5, 7.5, 10, 12.5, 15, 17.5, 20, 22.5, 25, 27.5, 30,
32.5, 35, 37.5, 40, 42.5, 45, 47.5, 50, 52.5, 55, 57.5, 60, 62.5,
65, 67.5, 70, 72.5, 75, 77.5, 80, 82.5, 85, 87.5, 90, 92.5, 95,
97.5, 100, 102.5, 105, 107.5, 110, 112.5, 115, 117.5, 120, 122.5,
125, 127.5, 130, 132.5, 135, 137.5, 140, 142.5, 145, 147.5, 150,
152.5, 155, 157.5, 160, 162.5, 165, 167.5, 170, 172.5, 175, 177.5,
180, 182.5, 185, 187.5, 190, 192.5, 195, 197.5, 200, 202.5, 205,
207.5, 210, 212.5, 215, 217.5, 220, 222.5, 225, 227.5, 230, 232.5,
235, 237.5, 240, 242.5, 245, 247.5, 250, 252.5, 255, 257.5, 260,
262.5, 265, 267.5, 270, 272.5, 275, 277.5, 280, 282.5, 285, 287.5,
290, 292.5, 295, 297.5, 300, 302.5, 305, 307.5, 310, 312.5, 315,
317.5, 320, 322.5, 325, 327.5, 330, 332.5, 335, 337.5, 340, 342.5,
345, 347.5, 350, 352.5, 355, 357.5 ; }
Again, it seems raster is reporting range(x)=[-1.25, 358.75] only because resolution(x)=2.5?. I'll defer to its superior wisdom. Given my need to integrate this data with other more conventional, longitudes=[-180, +180] datasets, it seems better to fix this data than to fix wrld_simpl. So I'm wondering, how best to "shift" my data from longitudes=[0, 360] to longitudes=[-180, +180]? TIA, Tom Roche <Tom_Roche at pobox.com>
An embedded and charset-unspecified text was scrubbed... Name: not available URL: <https://stat.ethz.ch/pipermail/r-sig-geo/attachments/20121219/4ec151f2/attachment.pl>
An embedded and charset-unspecified text was scrubbed... Name: not available URL: <https://stat.ethz.ch/pipermail/r-sig-geo/attachments/20121219/35265d5e/attachment.pl>
Tom Roche Wed, 19 Dec 2012 13:52:38 -0500
How to best/easiest range-shift lon-lat data?
Matt Landis Wed, 19 Dec 2012 14:09:57 -0500
[raster::shift] will work, but [raster::rotate] is even easier
Thanks! Now this code-------------------------------------------------
global.proj <- CRS('+proj=longlat +ellps=WGS84')
global.raster <- brick(in.fp, varname=data.var.name, crs=global.proj)
global.raster
# class : RasterBrick
# dimensions : 96, 144, 13824, 56 (nrow, ncol, ncell, nlayers)
# resolution : 2.5, 1.894737 (x, y)
# extent : -1.25, 358.75, -90.94737, 90.94737 (xmin, xmax, ymin, ymax)
# coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
layers.n <- nbands(global.raster)
# make longitudes NOT zero-based
global.raster <-
rotate(global.raster,
overwrite=TRUE) # else levelplot does one layer per page?
global.raster
# class : RasterBrick
# dimensions : 96, 144, 13824, 56 (nrow, ncol, ncell, nlayers)
# resolution : 2.5, 1.894737 (x, y)
# extent : -181.25, 178.75, -90.94737, 90.94737 (xmin, xmax, ymin, ymax)
# ...
# TODO: if using rasterVis::levelplot, get names with `layerNames`
global.df <- as.data.frame(global.raster)
names(global.df) <- formatC(as.numeric(sub(x=names(global.df), pattern="X", replacement="", fixed=TRUE)), format="g", digits=sigdigs)
# debugging
summary(unlist(global.df))
# Min. 1st Qu. Median Mean 3rd Qu. Max.
# 0.01459 0.22440 0.48350 0.36550 0.48530 0.50210
# matches known stats
global.map <- map('world', plot=FALSE)
global.map.shp <- map2SpatialLines(global.map, proj4string=global.proj)
# start plot driver
cat(sprintf(
'%s: plotting %s (may take awhile! TODO: add progress control)\n',
this.fn, pdf.fp))
pdf(file=pdf.fp, width=5, height=100)
# using rasterVis::levelplot, not lattice::levelplot
levelplot(global.raster,
margin=FALSE,
names.attr=names(global.df),
layout=c(1,layers.n), # all layers in one "atmospheric" column
) + layer(sp.lines(global.map.shp, lwd=0.8, col='darkgray'))
dev.off()
--------------------------------------------------------------produces
https://docs.google.com/open?id=0BzDAFHgIxRzKaGhMOE5BblpEeVE
with both map hemispheres overlaying the data.
thanks all! Tom Roche <Tom_Roche at pobox.com>
An embedded and charset-unspecified text was scrubbed... Name: not available URL: <https://stat.ethz.ch/pipermail/r-sig-geo/attachments/20121219/32e0395b/attachment.pl>