Re: [R-sig-Geo] problems drawing a world map on a levelplot
Hello,
Perhaps the problem could be related with the CRS setting you are using.
Instead of:
global.map.proj - CRS('+proj=latlon +ellps=WGS84')
you should try:
global.map.proj - CRS('+proj=longlat +ellps=WGS84')
Besides, following your code I think you are not setting the CRS of your
RasterBrick.
I have tried successfully the next toy example:
global.map.proj - CRS('+proj=longlat +ellps=WGS84')
global.map - map('world', plot=FALSE)
global.map.shp - map2SpatialLines(global.map, proj4string=global.map.proj)
r - raster()
r - init(r, fun=function(x)1)
projection(r) - global.map.proj
levelplot(r) + layer(sp.lines(global.map.shp))
Besides, you can modify the theme to fill the sea (NA values) with color:
myTheme - modifyList(RdBuTheme(),
list(panel.background=list(col='skyblue3')))
## you won't see any differences here because r does not contains NA cells.
levelplot(r, par.settings=myTheme) + layer(sp.lines(global.map.shp))
Finally, you may find useful these sources for the countries boundaries:
http://www.diva-gis.org/Data
http://www.naturalearthdata.com/downloads/110m-cultural-vectors/
Best,
Oscar.
2012/12/19 Michael Sumner mdsum...@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@polygons)))
dfw - data.frame(x = seq_len(length(worldw@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 mdsum...@gmail.com
wrote:
I meant wrld_simpl2 version like the maps package has.
On Tue, Dec 18, 2012 at 8:54 PM, Michael Sumner mdsum...@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)
[R-sig-Geo] Upcoming training courses: ISRIC's Spring School 2013 and GEOSTAT Quebec City 2013
Dear R-sig-geo, Registrations for the two upcoming 5-day training courses are now open: --- Hands-on Global Soil Information Facilities, 22-26 April 2013, Wageningen University The purpose of the Spring School is to introduce participants to software for soil data analysis and visualisation, digital soil mapping and soil-web services. The focus is on using the Free and Open Source Software for statistical and geographical computing for soil analysis and mapping: R (the main programming environment), SAGA GIS, Python and Google Earth. Upon completing the course, participants will be able to use the software and data sets for generating soil GIS layers following the international standards (e.g. GlobalSoilMap specifications) and for land use planning. Registration deadline: 20th of January 2013. URL: http://www.isric.org/content/isrics-spring-school-2013 --- GEOSTAT Quebec City 2013, 26 May - 1 of June, 2013, Quebec City The GEOSTAT Quebec city 2013 summer school is 8th in a series of summer schools organized by R and OS GIS developers and enthusiasts. GEOSTAT aims at PhD students and R-sig-geo enthusiasts in a range of environmental and GIS sciences, with a special focus on analyzing spatio-temporal gridded data in R and connected OS GIS software. The participants will learn how to import and organize space-time data i.e. time series of rasters and how to program statistical and geographical analysis using a combination of R and OS GIS functionality. The school will be hosted at the Faculty of Forestry, Geomatics and Geography, on the campus of Laval University, a few kilometers from the heart of Quebec City, Canada. Registration deadline: 1st of February 2013. URL: http://geostat-course.org/Quebec_2013 PS: If you are aware of any such (future or past!) training course or event that is related to the R-sig-geo topics, please consider adding it to the GEOSTAT timeline at http://geostat-course.org/timeline/nodes T. Hengl http://www.wewur.wur.nl/popups/vcard.aspx?id=HENGL001 ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
[R-sig-Geo] reading PostGIS table into sp data frame
I wanted to create a global map with squares in lat-lon. I have PostGIS tables
to define these squares – but I haven’t been able to figure out an efficient
way of reading those tables into R. The code I am using now is:
crs - CRS(+proj=longlat +ellps=WGS84)
s - paste(select id, st_astext(geom) as geom from geo.cs10d;, sep=)
r - dbGetQuery(con, s)
p - readWKT(r$geom[1],id=r$id[1],p4s=crs)
for(i in 2:length(r$id)){
p - rbind(p, readWKT(r$geom[i], id=r$id[i], p4s=crs))
}
where geo.cs10d is the table with squares, id the primary key of the table, and
geom the binary geometry field.
The code above works fine for the larger squares, such as 10 degrees, of which
I only need 648 to cover the globe. For finer resolutions, the above takes just
too long – I assume because the rbind function rewrites the whole sp object
each time it executes. I’ve seen other R scripts that initiate an empty data
frame of the correct length to go round similar problems with the rbind
function; I haven’t been able to find an equivalent for spatial polygons. How
can I initiate an empty data frame with the right structure, and the right
length?
A preferable solution would be if there would be a single function to load a
complete PostGIS table, rather than having to load the polygons one by one in a
loop. Is there such a function?
I’m using PostgreSQL 8.4, PostGIS 1.5, R 2.15.2, platform x86_64-w64-mingw32;
IDE is StatET 3.0.1 plugin for Eclipse 3.7.2.
Any help would be much appreciated.
Edward
___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] reading PostGIS table into sp data frame
2012/12/19 Edward Vanden Berghe evber...@gmail.com
I wanted to create a global map with squares in lat-lon. I have PostGIS
tables to define these squares but I havent been able to figure out an
efficient way of reading those tables into R. The code I am using now is:
crs - CRS(+proj=longlat +ellps=WGS84)
s - paste(select id, st_astext(geom) as geom from geo.cs10d;,
sep=)
r - dbGetQuery(con, s)
p - readWKT(r$geom[1],id=r$id[1],p4s=crs)
for(i in 2:length(r$id)){
p - rbind(p, readWKT(r$geom[i], id=r$id[i], p4s=crs))
}
where geo.cs10d is the table with squares, id the primary key of the
table, and geom the binary geometry field.
The code above works fine for the larger squares, such as 10 degrees, of
which I only need 648 to cover the globe. For finer resolutions, the above
takes just too long I assume because the rbind function rewrites the
whole sp object each time it executes. Ive seen other R scripts that
initiate an empty data frame of the correct length to go round similar
problems with the rbind function; I havent been able to find an equivalent
for spatial polygons. How can I initiate an empty data frame with the right
structure, and the right length?
A preferable solution would be if there would be a single function to load
a complete PostGIS table, rather than having to load the polygons one by
one in a loop. Is there such a function?
Im using PostgreSQL 8.4, PostGIS 1.5, R 2.15.2, platform
x86_64-w64-mingw32; IDE is StatET 3.0.1 plugin for Eclipse 3.7.2.
Any help would be much appreciated.
Edward
Use postgis st_transfom to convert your table in epsg:4326, I would suggest
to use a view for that.
Then use readOGR (in package rgdal) to read the table/view (geom and
attributes) in a sp object with :
sp.object - readOGR(PG:dbname=your_db host=your_host user=username
password=xxx, geo.cs10d )
regards
-r
--
*L'unica speranza di catarsi, ammesso che ne esista una, resta affidata
all'istinto di ribellione, alla rivolta non isterilita in progetti, alla
protesta violenta e viscerale. (V. Evangelisti)
*
[[alternative HTML version deleted]]
___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] problems drawing a world map on a levelplot
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_ro...@pobox.com
___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
[R-sig-Geo] help for spdep
Dear all. I try to fit conditional auto regressive models (CAR and SAR) in package spdep. Also, I have fitted some other models like GLM, Empirical Bayes and ... My program to CAR and SAR is: esar1f1 - spautolm(IMR.m~ 0+PCFP+Sup+B.F,data =data1, listw=nb2listw(nb6, style=W), family=SAR, method=full, verbose=TRUE) summary(esar1f1) esar1f2 - spautolm(IMR.m~ 0+PCFP+Sup+B.F,data =data1, listw=nb2listw(nb6, style=W), family=CAR, method=full, verbose=TRUE) esar1f1 Call: spautolm(formula = IMR.m ~ 0 + PCFP + Sup + B.F, data = data1, listw = nb2listw(nb6, style = W), family = SAR, method = full, verbose = TRUE) Coefficients: PCFPSupB.F lambda 0.1711382 -0.2262700 0.2414336 -1.7946991 Log likelihood: -56.24405 What is lambda?? I couldn't find a good refrences to understand tese model!!! How can I predict this model for new data?? How can I cross validate results of CAR and SAR?? -- Saman Monfared Msc, Department of Statistics, Shiraz University, Shiraz 71454, Iran Email: samanmonfar...@gmail.com Tel: +98 917 5305167 ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
[R-sig-Geo] Runtime error when opening netCDF
Dear all, I have a similar error reported by Robert Buitenwerf here: http://grokbase.com/t/r/r-sig-geo/11a3b5ms2c/netcdf-to-raster-error I am trying to read/open a netCDF file using the ncdf package, but it keeps giving me the Microsoft Visual C++ Runtime Library error: This application has requested the Runtime to terminate it in an unusual way... I am trying to open the CRU TS3.1 climate data (open.ncdf(nc file)). I am able to open a netCDF file smaller than the 2.5GB CRU data so is this a size issue? I am running R on Windows 7 64 bit. Sincerely, Prajjwal [[alternative HTML version deleted]] ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] Runtime error when opening netCDF
2 suggestions: Panday, Prajjwal Wed Dec 19 19:19:30 CET 2012 I am trying to read/open a netCDF file using the ncdf package, Use package=ncdf4, not package=ncdf: the latter is no longer supported, per the developer (of both packages). but it keeps giving me the Microsoft Visual C++ Runtime Library error: Use linux. HTH, 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
Re: [R-sig-Geo] problems drawing a world map on a levelplot
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_ro...@pobox.com
___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
Re: [R-sig-Geo] Runtime error when opening netCDF
Prajjwal, You can also read these files using ncdf with the the 32-bit version of R (on windows 64). Robert On Wed, Dec 19, 2012 at 10:39 AM, Tom Roche tom_ro...@pobox.com wrote: 2 suggestions: Panday, Prajjwal Wed Dec 19 19:19:30 CET 2012 I am trying to read/open a netCDF file using the ncdf package, Use package=ncdf4, not package=ncdf: the latter is no longer supported, per the developer (of both packages). but it keeps giving me the Microsoft Visual C++ Runtime Library error: Use linux. HTH, 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 [[alternative HTML version deleted]] ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] problems drawing a world map on a levelplot
library(raster)
?shift
On Wed, Dec 19, 2012 at 1:52 PM, Tom Roche tom_ro...@pobox.com wrote:
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
Re: [R-sig-Geo] rasterize func: polygon values to raster cells - error msg
Gustaf, Thanks, I think I found the bug. It happens when there is no overlap at all between polygons and raster. I think it is fixed in raster 2.0-39 (you can install from here https://r-forge.r-project.org/R/?group_id=294). Robert On Tue, Dec 18, 2012 at 10:35 AM, Gustaf Granath gustaf.gran...@ebc.uu.sewrote: Robert, Thanks for looking at this. Unfortunately, updating R and Raster did not help. Here is info you requested. Hope that helps! Gustaf tt-rasterize(SpPDF,tes,field=idd) Error in writeValues(out, fun(n), tr$row[i]) : error in evaluating the argument 'v' in selecting a method for function 'writeValues': Error in fun(n) : unused argument(s) (n) traceback() 6: writeValues(out, fun(n), tr$row[i]) 5: init(rstr, function() NA) 4: .polygonsToRaster(x, y, field = field, fun = fun, background = background, mask = mask, update = update, updateValue = updateValue, filename = filename, getCover = getCover, silent = silent, ...) 3: .local(x, y, ...) 2: rasterize(SpPDF, tes, field = idd) 1: rasterize(SpPDF, tes, field = idd) show(SpPDF) class : SpatialPolygonsDataFrame nfeatures : 58824 extent : -11.08333, 31.91753, 33.91667, 71.91743 (xmin, xmax, ymin, ymax) coord. ref. : +proj=longlat +datum=WGS84 nvariables : 92 names : X.Lon, Lat, JulT00, P1900, NH1900, N1900, JulT10, P1910, NH1910, N1910, JulT20, P1920, NH1920, N1920, JulT30, ... min values : -10.5, 34.8333,0.3, 260, 0.103, 0.0163, 0, 261, 0.103, 0.0184,0.1, 245, 0.118, 0.0215,0.3, ... max values : 31.8333, 71, 27.2, 3267, 14.837, 1.5572, 27.2, 3226, 11.927, 1.2981, 27.8, 3303, 12.945, 1.4166, 27.5, ... Warning messages: 1: closing unused connection 4 (C:/Users/veGustaf/AppData/Local/Temp/R_raster_tmp/veGustaf/raster_tmp_2012-12-18_192108_37929.gri) 2: closing unused connection 3 (C:/Users/veGustaf/AppData/Local/Temp/R_raster_tmp/veGustaf/raster_tmp_2012-12-18_191936_11753.gri) show(tes) class : RasterLayer dimensions : 12768, 11050, 141086400 (nrow, ncol, ncell) resolution : 250, 250 (x, y) extent : 2782250, 5544750, 2220250, 5412250 (xmin, xmax, ymin, ymax) coord. ref. : +proj=longlat +datum=WGS84 sessionInfo() R version 2.15.2 (2012-10-26) Platform: i386-w64-mingw32/i386 (32-bit) locale: [1] LC_COLLATE=Swedish_Sweden.1252 LC_CTYPE=Swedish_Sweden.1252 [3] LC_MONETARY=Swedish_Sweden.1252 LC_NUMERIC=C [5] LC_TIME=Swedish_Sweden.1252 attached base packages: [1] stats graphics grDevices utils datasets methods base other attached packages: [1] raster_2.0-31 sp_1.0-2 loaded via a namespace (and not attached): [1] grid_2.15.2 lattice_0.20-10 tools_2.15.2 On 2012-12-17 07:32, Robert J. Hijmans wrote: Gustaf, When reporting that something does not work, please include the result of your sessionInfo() Just after the error occurs. traceback() can also be helpful. In this case it would appear that that you are not using the current version of raster. If that is the case, please update.packages() first (perhaps after getting the most recent version of R). Some more info on the SpatialPolygonsDataFrame that causes the error would also be helpful (use the show function after raster is loaded). Robert On Wed, Dec 12, 2012 at 11:55 PM, Gustaf Granath gustaf.gran...@ebc.uu.se wrote: Hi, I have used the rasterize function to get values from a spatialPolyDataFrame variable (@data$variable) to raster: rasterize(SPDFobj, rasterObj, field=variable) For a specific case a get a, to me, strange error message . Error in writeValues(out, fun(n), tr$row[i]) : error in evaluating the argument 'v' in selecting a method for function 'writeValues': Error in fun(n) : unused argument(s) (n) It works with one SpatialPolyDataFrame that Im working with, but not with another one and I cannot figure out why. Any ideas of what is going on? Unfortunately, I have not been able to reproduce the error with simulated data (and my data files are huge) but this is essentially what I have done (but the code below works - actually the raster layer in this code works with my SPDFobj): p1 - rbind(c(-180,-20), c(-140,55), c(10, 0), c(-140,-60), c(-180,-20)) p2 - rbind(c(-10,0), c(140,60), c(160,0), c(140,-55), c(-10,0)) p3 - rbind(c(-125,0), c(0,60), c(40,5), c(15,-45), c(-125,0)) pols - SpatialPolygons( list( Polygons(list(Polygon(p1), Polygon(hole)), 1), Polygons(list(Polygon(p2)), 2), Polygons(list(Polygon(p3)), 3))) pols -SpatialPolygonsDataFrame(pols, data.frame(W=1:3)) r - raster(ncol=80, nrow=90) test-rasterize(pols, r, field=W) Here is an overview of the objects I am using. summary(SPDFobj) Object of class SpatialPolygonsDataFrame Coordinates: min max x -11.08333 31.91753 y 33.91667 71.91743 Is projected: FALSE proj4string : [+proj=longlat +datum=WGS84] Data attributes:
Re: [R-sig-Geo] Runtime error when opening netCDF
Robert Yes I was able to read the file using ncdf. However, I am not able to load the entire variable in R-32bit using get.var.ncdf command. I get the error cannot allocate vector of size 2.5 Gb. Also, would this be a problem down the road if I am working with large datasets in R-32bit? Thanks, Prajjwal -- View this message in context: http://r-sig-geo.2731867.n2.nabble.com/Runtime-error-when-opening-netCDF-tp7582014p7582021.html Sent from the R-sig-geo mailing list archive at Nabble.com. ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] Runtime error when opening netCDF
Prajjwal, You could use 32bit R and the raster package to write the data to new files; and then continue with 64bit R. Robert On Wed, Dec 19, 2012 at 11:16 AM, ppanday ppan...@clarku.edu wrote: Robert Yes I was able to read the file using ncdf. However, I am not able to load the entire variable in R-32bit using get.var.ncdf command. I get the error cannot allocate vector of size 2.5 Gb. Also, would this be a problem down the road if I am working with large datasets in R-32bit? Thanks, Prajjwal -- View this message in context: http://r-sig-geo.2731867.n2.nabble.com/Runtime-error-when-opening-netCDF-tp7582014p7582021.html Sent from the R-sig-geo mailing list archive at Nabble.com. ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo [[alternative HTML version deleted]] ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
[R-sig-Geo] Make raster from oddly organized ASCII
Dear all, Today I received some ASCII files that need to be converted to raster. They are organized in a rather odd way, and I cannot think of a method to read and organize the data to make a raster. Please take a look at a sample file: sample-readLines(url(http://dl.dropbox.com/u/27700634/gghswmp.lis;)) Useful information that I know about the data: - global coverage at 1 degree; - each line's 2nd value corresponds to the latitude coordinate; - the records are ordered from north (lat = 89) to south (lat = 90) within each hemisphere, the eastern hemisphere (lon=0) appearing first and the western (lon=180) second; - lon is the longitude and lat is the latitude of the southwest corner of the first 1x1 degree cell. Initially I thought of creating a XYZ data.frame and then use rasterfromXYZ to produce the raster. Could anyone help me to create a XYZ data.frame from the supplied file? Many thanks in advance, Thiago. ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] problems drawing a world map on a levelplot
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. MedianMean 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_ro...@pobox.com
___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] Runtime error when opening netCDF
Robert, from what I understand, use raster package in 32-bit to open and write as a new netCDF file (using writeRaster) and then open in 64-bit using ncdf. This way I should be able to bypass size limit issue if I come across large datasets. So far, I have read the netCDF file as raster using brick() command. In my next steps, I am trying to compute monthly anomalies of the netcdf data and extract monthly means using another raster layer. Any suggestions regarding these steps would be greatly appreciated. Thanks, Prajjwal -- View this message in context: http://r-sig-geo.2731867.n2.nabble.com/Runtime-error-when-opening-netCDF-tp7582014p7582027.html Sent from the R-sig-geo mailing list archive at Nabble.com. ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
[R-sig-Geo] Raster crosstab output: -1 values?
Hi everyone, I'm computing a very large crosstab analysis (using the raster package) for around 700 tiles of 25m resolution vegetation raster data. Essentially, I want to figure out how much forest is contained within each land use, land tenure and jurisdictional boundary (eight Australian States) over a 39 year time period (20 years with data). The land tenure, land use and State layers have all been cut up into the same extent, resolution and projection as the forest cover data. I then create a raster stack of the forest cover, State and land tenure-use (the latter was combined into one raster layer) tiles, and subsequently a crosstab. This process is repeated for each tile and year (~700 loops, taking a rather long time). The forest cover data (in this particular year - 1972) has values of only 0 and 1 (non-forest/forest). However, I'm finding in some cases the crosstab output contains data in a -1 category, despite this value not being part of the original forest cover data. I'm at a loss to explain this; all datasets have the same number of cells, resoluton, and no missing data. Would really appreciate any insights! Cheers, Megan Data details below: Forest Cover (ttile): class : RasterLayer dimensions : 16000, 28360, 45376 (nrow, ncol, ncell) resolution : 0.00025, 0.00025 (x, y) extent : 112.91, 120, -28, -24 (xmin, xmax, ymin, ymax) coord. ref. : +proj=longlat +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +no_defs data source : D:/Analysis/R/Data/forestextent_reg/forest_v8_72_11/sg50/forst72.grd names : layer values : 0, 1 (min, max) Raster Attribute Table fields : ID COUNT min : 0 30065009 max : 1 423694991 State (st_tex): class : RasterLayer dimensions : 16000, 28360, 45376 (nrow, ncol, ncell) resolution : 0.00025, 0.00025 (x, y) extent : 112.91, 120, -28, -24 (xmin, xmax, ymin, ymax) coord. ref. : +proj=longlat +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +no_defs data source : D:/Analysis/R/Data/sttiles/forest_v8_72_11/stsg50.grd names : states_p values : 8, 8 (min, max) Land Use/Tenure (lu_ten_tex): class : RasterLayer dimensions : 16000, 28360, 45376 (nrow, ncol, ncell) resolution : 0.00025, 0.00025 (x, y) extent : 112.91, 120, -28, -24 (xmin, xmax, ymin, ymax) coord. ref. : +proj=longlat +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +no_defs data source : D:/Analysis/R/Data/lutentiles/forest_v8_72_11/lutensg50.grd names : lu_ten values : 1, 73 (min, max) Raster Attribute Table fields : ID COUNT min : 1 -1816924096 max : 73 1681761600 Relevant code: # Create raster stack between forestextent, lu and state x -stack(ttile, lu_ten_tex, st_tex) #Crosstab will calculate the frequency of each level (0=non-forest, 1=primforest, 2=regforest) in each land use/tenure/state vegfreq -crosstab(x, digits=0, long=FALSE, progress=text) vegfreq , , states_p = 8 lu_ten layer 0 1 311136163 33 0 400975 2059222 2355223 28336430 13144612 2565663 316999241 42855 1 15814 1098 19824 6497248 3958995761517 17902798 345 -111 1280 153 3522 39320 15961 0 lu_ten layer43732353715112 22 0 8447 16777 64401 50205 19375 17463334230 664199 1 1153 7223 66799 995 31825 137396970 126201 -1 0 0 0 0 0 0 0 0 lu_ten layer62 232 0 310439 1586 13062 1 20156114 2938 -1 0 0 0 -- View this message in context: http://r-sig-geo.2731867.n2.nabble.com/Raster-crosstab-output-1-values-tp7582028.html Sent from the R-sig-geo mailing list archive at Nabble.com. ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] R 3.0.0 and spatial classes
Is there interest in something like this code for generating SpatialLines?
I use a variation of this fairly regularly:
SimpleSpatialLines - function(coords.1, coords.2, proj4string =
CRS(as.character(NA))){
coords.1 - coordinates(coords.1)
coords.2 - coordinates(coords.2)
if (nrow(coords.1) != nrow(coords.2)) stop(coords.1 must be of same
length as coords.2)
coords.both - cbind(coords.1, coords.2)
coords.both - cbind(coords.both, 1:nrow(coords.both))
lines - apply(X=coords.both,
MARGIN=1,
FUN=function(x){
Lines(list(Line(cbind(c(x[1], x[3]),
c(x[2], x[4],
ID=x[5])})
SpatialLines(LinesList=lines, proj4string=proj4string)
}
Regards,
Francis Markham
Australian National University
francis.mark...@anu.edu.au
On 17 December 2012 06:58, Barry Rowlingson b.rowling...@lancaster.ac.ukwrote:
On Sun, Dec 16, 2012 at 6:53 PM, Chris English sgl...@hotmail.com wrote:
Edzer:
On Tue, Jul 17, 2012 at 2:34 PM, Agustin Lobo ~ from Barry Rowlingson
wrote:
To convert to SpatialLines, get the coordinates and build in the
usual convoluted manner:
s=data.frame(x=runif(10),y=1:10,z=rnorm(10))
coordinates(s)=~x+y
L = SpatialLines(list(Lines(list(Line(coordinates(s))),X)))
plot(L)
Convoluted is not the same as orphaned, certainly, but one gets the
sense that
'Line' owes its existence to matters of plotting rather than line as
line, independent of
drawing it, and this may have some import upon line analysis and the
possibility of
arriving at topology and dispensing with shared lines and the like.
Line is a non-spatial Line, made up from an ordered set of (x,y)
coordinate pairs. It cannot have a coordinate system assigned to it.
Lines is a list of Lines, making a non-spatial set of Line
segments. It too cannot have a coordinate system assigned.
SpatialLines is a spatial set of Lines, for when you have a number
of features each of which may be composed of several disconnected
segments. It can have a coordinate system assigned.
You are making a SpatialLines object where a single feature has a
single line segment. The slightly annoying
list(Lines(list(Line(... dance is unavoidable because the lists
are necessary since the lists can have more than one element, but you
could easy create a 'SimpleSpatialLine' function that did all that. It
might even make sp for R 3.0.0! (SimpleSpatialPolygon might be handy
too...)
Another possibility might be to write methods for SpatialLines that
takes a matrix and skips the complexity for simple cases...
Barry
___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
[[alternative HTML version deleted]]
___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
