Re: [R-sig-Geo] Raster package writeRaster format=NCF
Elizabeth, Creating a raster object from a nc file is probably easiest done like this: r - raster(filename, band=1) Or to get the all: b - brick(filename) then: ba - aggregate(b, ) Normally, you should be able to do this to save it as ncdf: ba - writeRaster(ba, filename='out.nc') Or in one step: ba - aggregate(b, filename='out.nc') However, you have uncovered a bug in writing to a ncdf file when values are not in memory, I'll think I have fixed that in the development version (available from R-Forge in 24 hrs). For now, a work around is to do 'readAll' before you use writeRaster (this will fail if the file is too large, but that is probably not the case) r - readAll(r) r - writeRaster(... Robert On Tue, Oct 5, 2010 at 6:20 AM, Elizabeth Crisfield ea...@psu.edu wrote: I have read a 10 band .nc file into R using *open.ncdf* and *get.var.ncdf*and I've read each band into a raster format e.g. *raster(data[,,1])*. I used *aggregate* to upscale the data, and now I want to write it back out to netcdf format. For some reason, I can write each band out separately (e.g. *writeRaster(newdataband1, newdataband1, format=CDF)* but if I *stack* the data first, I get an error when I try to write the stack out to netcdf. The error says *CDF is not a supported file format. See writeFormats()*. Sure enough, *writeFormats()* gives me a long list of formats that are enabled, and CDF is not in the list. ? I am running R (R 2.11.1 GUI 1.34 Leopard build 64-bit (5589)) on Mac OS X (10.6.4) with up-to-date packages of ncdf and Raster. I've also installed netcdf through MacPorts and the libraries are in /opt/local rather than /usr/local (I'm wondering if I have a library or path problem?) How can I get the writeRaster to add CDF to the list of writeFormats? Is there another way to reconstruct this netcdf file from the 10 raster layers in R? (*create.ncdf*?) Alternatively, since I think I have them written out to separate netcdfs I could maybe bind them back together using another utility outside of R? Thanks in advance. Elizabeth Elizabeth Crisfield Research Assistant Geography Department Penn State University phone 814-777-3395 ea...@psu.edu www.elizabethcrisfield.com [[alternative HTML version deleted]] ___ R-sig-Geo mailing list R-sig-Geo@stat.math.ethz.ch https://stat.ethz.ch/mailman/listinfo/r-sig-geo ___ R-sig-Geo mailing list R-sig-Geo@stat.math.ethz.ch https://stat.ethz.ch/mailman/listinfo/r-sig-geo
[R-sig-Geo] accessing raster::lineValues
Recently, Robert kindly added a function to raster package (raster_1.5-14 from RForge) called lineValues. After loading the new package, I am unable to access (or use) the function. Any ideas what I might be doing wrong? rst.cln - raster::lineValues(lns = rst.lines, rst.cellnum) Error: 'lineValues' is not an exported object from 'namespace:raster' lineValues Error: object 'lineValues' not found raster::lineValues Error: 'lineValues' is not an exported object from 'namespace:raster' sessionInfo() R version 2.11.1 (2010-05-31) i386-pc-mingw32 locale: [1] LC_COLLATE=Slovenian_Slovenia.1250 LC_CTYPE=Slovenian_Slovenia.1250 [3] LC_MONETARY=Slovenian_Slovenia.1250 LC_NUMERIC=C [5] LC_TIME=Slovenian_Slovenia.1250 attached base packages: [1] stats graphics grDevices utils datasets methods base other attached packages: [1] raster_1.5-14spatstat_1.20-2 deldir_0.0-12mgcv_1.6-2 [5] adehabitatMA_0.1 gpclib_1.5-1 sp_0.9-66 loaded via a namespace (and not attached): [1] grid_2.11.1lattice_0.18-8 Matrix_0.999375-39 nlme_3.1-96 [5] tools_2.11.1 Cheers, Roman -- In God we trust, all others bring data. [[alternative HTML version deleted]] ___ R-sig-Geo mailing list R-sig-Geo@stat.math.ethz.ch https://stat.ethz.ch/mailman/listinfo/r-sig-geo
[R-sig-Geo] increment and distance in correlog
Hello again,
okay, i found the ncf package and the correlog function.
# correlog(x, y, z, w = NULL, increment, resamp = 100, latlon =
FALSE, na.rm = FALSE, quiet = FALSE)
my short question is:
if latlon = T and the positions are "degree, decimal-minutes"
then increment has the unit of "km" and distance bins are in "km"
as well - is that correct ?
I found it in the sourcecode :
#then generating geographic distances
if(latlon){
#these are geographic distances from lat-lon
coordinates
dmat - matrix(0, nrow = n, ncol = n)
for(i in 1:(n-1)) {
for(j in (i+1):n) {
dmat[j, i] - gcdist(x[i], y[i],
x[j], y[j])
dmat[i, j] - dmat[j, i]
with gcdist =
r - 360/(2 * pi)
lon1 - x1 / r
lat1 - y1 / r
lon2 - x2 / r
lat2 - y2 / r
dlon - lon2 - lon1
dlat - lat2 - lat1
a - (sin(dlat/2))^2 + cos(lat1) * cos(lat2) *
(sin(dlon/2))^2
c - 2 * atan2( sqrt(a), sqrt(1-a) )
return(6370 * c)
}
but i cannot the answer from thisdoes gcdist deliver km units ?
Cheers!
jens
Am 05.10.2010 12:57, schrieb Jens Floeter:
Dear
list members,
i need a hint - i am stuck on my search for a function.
I would like to calculate cross-correlograms between two variables
(predator and prey densities) and use this as an index of spatial
overlap.
I used the spdep pakcage and use e.g.
sp.correlogram(xy.data.nb, prey, order = 5 , method = "I",
zero.policy=TRUE, style = "W")
sp.correlogram(xy.data.nb, predator, order = 5 , method = "I",
zero.policy=TRUE, style = "W")
to calculate the (auto-) correlograms for predator and prey
individually.
But how do i combine them to calculate the cross-correlogram ?
Any hint would be highly appreciated...
Best wishes !
jens
___
R-sig-Geo mailing list
R-sig-Geo@stat.math.ethz.ch
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
___
R-sig-Geo mailing list
R-sig-Geo@stat.math.ethz.ch
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] increment and distance in correlog
Jens,
in your code snippet below, the number 6370 rings a bell. According to
http://en.wikipedia.org/wiki/Earth_radius it is close to the average
earth radius. In km. Given that there's just one such number, it seems
the code assumes a sphere is close enough to approximate the earth.
On 10/06/2010 12:37 PM, Jens Floeter wrote:
Hello again,
okay, i found the ncf package and the correlog function.
# correlog(x, y, z, w = NULL, increment, resamp = 100, latlon = FALSE, na.rm
=
FALSE, quiet = FALSE)
my short question is:
*if latlon = T and the positions are degree, decimal-minutes then increment
has the unit of km and distance bins are in km as well - is that correct
?*
I found it in the sourcecode :
#then generating geographic distances
if(latlon){
#these are geographic distances from lat-lon coordinates
dmat - matrix(0, nrow = n, ncol = n)
for(i in 1:(n-1)) {
for(j in (i+1):n) {
dmat[j, i] - gcdist(x[i], y[i], x[j], y[j])
dmat[i, j] - dmat[j, i]
with gcdist =
r - 360/(2 * pi)
lon1 - x1 / r
lat1 - y1 / r
lon2 - x2 / r
lat2 - y2 / r
dlon - lon2 - lon1
dlat - lat2 - lat1
a - (sin(dlat/2))^2 + cos(lat1) * cos(lat2) * (sin(dlon/2))^2
c - 2 * atan2( sqrt(a), sqrt(1-a) )
return(6370 * c)
}
but i cannot the answer from thisdoes gcdist deliver km units ?
Cheers!
jens
Am 05.10.2010 12:57, schrieb Jens Floeter:
Dear list members,
i need a hint - i am stuck on my search for a function.
I would like to calculate cross-correlograms between two variables (predator
and prey densities) and use this as an index of spatial overlap.
I used the spdep pakcage and use e.g.
sp.correlogram(xy.data.nb, prey, order = 5 , method = I, zero.policy=TRUE,
style = W)
sp.correlogram(xy.data.nb, predator, order = 5 , method = I,
zero.policy=TRUE, style = W)
to calculate the (auto-) correlograms for predator and prey individually.
But how do i combine them to calculate the cross-correlogram ?
Any hint would be highly appreciated...
Best wishes !
jens
___
R-sig-Geo mailing list
R-sig-Geo@stat.math.ethz.ch
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
___
R-sig-Geo mailing list
R-sig-Geo@stat.math.ethz.ch
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
--
Edzer Pebesma
Institute for Geoinformatics (ifgi), University of Münster
Weseler Straße 253, 48151 Münster, Germany. Phone: +49 251
8333081, Fax: +49 251 8339763 http://ifgi.uni-muenster.de
http://www.52north.org/geostatistics e.pebe...@wwu.de
___
R-sig-Geo mailing list
R-sig-Geo@stat.math.ethz.ch
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
[R-sig-Geo] Basic Info. requried
I want to read Iknos image and then want to see display in R. Just read and display nothing else. I try using raster package rasterfile-system.file(C:/iknos/edwards_airbase_usa_1m_tc.tif,package=raster) Then plot with the help of this command plot(rasterfile) Error in plot.window(...) : need finite 'ylim' values In addition: Warning messages: 1: In min(x) : no non-missing arguments to min; returning Inf 2: In max(x) : no non-missing arguments to max; returning -Inf Need help from where I can learn basic things, I already read raster package manual but unable to see the image. with Best Regards, Malik Shahzad Visiting Researcher National Institute of Informatics (NII) Tokyo, Japan Doctoral Student Asian Institute of Technology (AIT) Bangkok, Thailand +66-8-7676-5616 [[alternative HTML version deleted]] ___ R-sig-Geo mailing list R-sig-Geo@stat.math.ethz.ch https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] Basic Info. requried
Hi Malik, To a look at the rgdal package to read your tiff file and spplot from the sp package to plot it. cheers, Paul On 10/06/2010 01:36 PM, Malik Shahzad wrote: I want to read Iknos image and then want to see display in R. Just read and display nothing else. I try using raster package rasterfile-system.file(C:/iknos/edwards_airbase_usa_1m_tc.tif,package=raster) Then plot with the help of this command plot(rasterfile) Error in plot.window(...) : need finite 'ylim' values In addition: Warning messages: 1: In min(x) : no non-missing arguments to min; returning Inf 2: In max(x) : no non-missing arguments to max; returning -Inf Need help from where I can learn basic things, I already read raster package manual but unable to see the image. with Best Regards, Malik Shahzad Visiting Researcher National Institute of Informatics (NII) Tokyo, Japan Doctoral Student Asian Institute of Technology (AIT) Bangkok, Thailand +66-8-7676-5616 [[alternative HTML version deleted]] ___ R-sig-Geo mailing list R-sig-Geo@stat.math.ethz.ch https://stat.ethz.ch/mailman/listinfo/r-sig-geo -- Drs. Paul Hiemstra Department of Physical Geography Faculty of Geosciences University of Utrecht Heidelberglaan 2 P.O. Box 80.115 3508 TC Utrecht Phone: +3130 253 5773 http://intamap.geo.uu.nl/~paul http://nl.linkedin.com/pub/paul-hiemstra/20/30b/770 ___ R-sig-Geo mailing list R-sig-Geo@stat.math.ethz.ch https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] Basic Info. requried
I believe there is a mistake in your command. You skipped the raster command itself. From the raster help files here is the correct way: r - raster(system.file(external/test.grd, package=raster)) On Wed, Oct 6, 2010 at 7:36 AM, Malik Shahzad geoma...@live.com wrote: I want to read Iknos image and then want to see display in R. Just read and display nothing else. I try using raster package rasterfile-system.file(C:/iknos/edwards_airbase_usa_1m_tc.tif,package=raster) Then plot with the help of this command plot(rasterfile) Error in plot.window(...) : need finite 'ylim' values In addition: Warning messages: 1: In min(x) : no non-missing arguments to min; returning Inf 2: In max(x) : no non-missing arguments to max; returning -Inf Need help from where I can learn basic things, I already read raster package manual but unable to see the image. with Best Regards, Malik Shahzad Visiting Researcher National Institute of Informatics (NII) Tokyo, Japan Doctoral Student Asian Institute of Technology (AIT) Bangkok, Thailand +66-8-7676-5616 [[alternative HTML version deleted]] ___ R-sig-Geo mailing list R-sig-Geo@stat.math.ethz.ch https://stat.ethz.ch/mailman/listinfo/r-sig-geo [[alternative HTML version deleted]] ___ R-sig-Geo mailing list R-sig-Geo@stat.math.ethz.ch https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] Basic Info. requried
Indeed, and as Malik is not using a system file (files packaged with R to serve as examples), you can should leave that bit out and do: r -raster(C:/iknos/edwards_airbase_usa_1m_tc.tif') plot(r) Robert On Wed, Oct 6, 2010 at 4:44 PM, kapo coulibaly kmcou...@gmail.com wrote: I believe there is a mistake in your command. You skipped the raster command itself. From the raster help files here is the correct way: r - raster(system.file(external/test.grd, package=raster)) On Wed, Oct 6, 2010 at 7:36 AM, Malik Shahzad geoma...@live.com wrote: I want to read Iknos image and then want to see display in R. Just read and display nothing else. I try using raster package rasterfile-system.file(C:/iknos/edwards_airbase_usa_1m_tc.tif,package=raster) Then plot with the help of this command plot(rasterfile) Error in plot.window(...) : need finite 'ylim' values In addition: Warning messages: 1: In min(x) : no non-missing arguments to min; returning Inf 2: In max(x) : no non-missing arguments to max; returning -Inf Need help from where I can learn basic things, I already read raster package manual but unable to see the image. with Best Regards, Malik Shahzad Visiting Researcher National Institute of Informatics (NII) Tokyo, Japan Doctoral Student Asian Institute of Technology (AIT) Bangkok, Thailand +66-8-7676-5616 [[alternative HTML version deleted]] ___ R-sig-Geo mailing list R-sig-Geo@stat.math.ethz.ch https://stat.ethz.ch/mailman/listinfo/r-sig-geo [[alternative HTML version deleted]] ___ R-sig-Geo mailing list R-sig-Geo@stat.math.ethz.ch https://stat.ethz.ch/mailman/listinfo/r-sig-geo ___ R-sig-Geo mailing list R-sig-Geo@stat.math.ethz.ch https://stat.ethz.ch/mailman/listinfo/r-sig-geo
[R-sig-Geo] R Basic Info. requried
Dear thanks for help, In first line there is no problem rasterfile-system.file(C:/iknos/edwards_airbase_usa_1m_tc.tif,package=raster) This works fine without any error, It mean image has been loaded into variable raster file, but now how i can display raster file in R. I mean I want to see my image that is stored in this variable. When I use plot(rasterfile) it shows error. Before this I work with rImage package to load image(*.jpg) file and then show with the help of plot(image). it works there fine and I can see image but now I am unable to see my output. I hope you got my point. with Best Regards, Malik Shahzad Visiting Researcher National Institute of Informatics (NII) Tokyo, Japan Doctoral Student Asian Institute of Technology (AIT) Bangkok, Thailand +66-8-7676-5616 Date: Wed, 6 Oct 2010 17:49:35 +0300 Subject: Re: [R-sig-Geo] Basic Info. requried From: r.hijm...@gmail.com To: kmcou...@gmail.com CC: geoma...@live.com; r-sig-geo@stat.math.ethz.ch Indeed, and as Malik is not using a system file (files packaged with R to serve as examples), you can should leave that bit out and do: r -raster(C:/iknos/edwards_airbase_usa_1m_tc.tif') plot(r) Robert On Wed, Oct 6, 2010 at 4:44 PM, kapo coulibaly kmcou...@gmail.com wrote: I believe there is a mistake in your command. You skipped the raster command itself. From the raster help files here is the correct way: r - raster(system.file(external/test.grd, package=raster)) On Wed, Oct 6, 2010 at 7:36 AM, Malik Shahzad geoma...@live.com wrote: I want to read Iknos image and then want to see display in R. Just read and display nothing else. I try using raster package rasterfile-system.file(C:/iknos/edwards_airbase_usa_1m_tc.tif,package=raster) Then plot with the help of this command plot(rasterfile) Error in plot.window(...) : need finite 'ylim' values In addition: Warning messages: 1: In min(x) : no non-missing arguments to min; returning Inf 2: In max(x) : no non-missing arguments to max; returning -Inf Need help from where I can learn basic things, I already read raster package manual but unable to see the image. with Best Regards, Malik Shahzad Visiting Researcher National Institute of Informatics (NII) Tokyo, Japan Doctoral Student Asian Institute of Technology (AIT) Bangkok, Thailand +66-8-7676-5616 [[alternative HTML version deleted]] ___ R-sig-Geo mailing list R-sig-Geo@stat.math.ethz.ch https://stat.ethz.ch/mailman/listinfo/r-sig-geo [[alternative HTML version deleted]] ___ R-sig-Geo mailing list R-sig-Geo@stat.math.ethz.ch https://stat.ethz.ch/mailman/listinfo/r-sig-geo [[alternative HTML version deleted]] ___ R-sig-Geo mailing list R-sig-Geo@stat.math.ethz.ch https://stat.ethz.ch/mailman/listinfo/r-sig-geo
[R-sig-Geo] spatial correlogram
Hello I am trying to use the correlog function to estimate a spatial correlogram for the residuals of a logistic regression and I have run accross the following error. summary(binom1 - glm(Use~X20mslop+X20mdem+soilsst, family=binomial, + data=M60m2000NE_1.df)) correlog1.1 - correlog(M60m2000NE_1.df$X, M60m2000NE_1.df$Y, residuals(binom1),w=NULL, + na.rm=T, increment=1, resamp=0) Error in split.default(moran, dkl) : Group length is 0 but data length 0 Can anyone help me to understand what this error is referring to and suggest how I can fix it? Thanks Tara PS Please see the file uploaded below :) http://r-sig-geo.2731867.n2.nabble.com/file/n5607783/Modified60m2000TrainingSetANoEcotour_05_25B.csv Modified60m2000TrainingSetANoEcotour_05_25B.csv -- View this message in context: http://r-sig-geo.2731867.n2.nabble.com/spatial-correlogram-tp5607783p5607783.html Sent from the R-sig-geo mailing list archive at Nabble.com. ___ R-sig-Geo mailing list R-sig-Geo@stat.math.ethz.ch https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] Why can you not calculate Moran's i for glm() residuals?
Hello Jenn Unfortunately I can not offer any advice but I am facing very similar challenges and actively searching for answers. I have presence/absence data for a colonial species and I am trying to compare the performance of a logistic model with spatial models: 1) a glm including an autocovariate and 2) a glmm. Currently I have resorted to using spatial correlograms for visual inspection of residual spatial autocorrelation (interesting that we both came up with the same posible option). Not ideal and I have not found any information as to whether this is valid but I am stumped. I was wondering if you have found any more information on the topics you described and if so would you mind sharing what you found with me. I will certainly return the favour. Thanks so much and Good Luck! Sincerely Tara MSc Candidate -- View this message in context: http://r-sig-geo.2731867.n2.nabble.com/Why-can-you-not-calculate-Moran-s-i-for-glm-residuals-tp5568942p5608678.html Sent from the R-sig-geo mailing list archive at Nabble.com. ___ R-sig-Geo mailing list R-sig-Geo@stat.math.ethz.ch https://stat.ethz.ch/mailman/listinfo/r-sig-geo
[R-sig-Geo] spplot with two SpatialPolygonsDataFrame: 1) filling background polygon; 2) left justifying text
Greetings, I have two questions/problems related to spplot: 1. I am trying to plot one SpatialPolygonsDataFrame over another, and fill the background polygon with grey. I can't seem to fill the background polygon without having it cover over my supposed foreground polygon, even if I specify first = T (which fails): link1 - http://sites.google.com/site/ldemisc/spplot-question/demo.mdshp.RData; link2 - http://sites.google.com/site/ldemisc/spplot-question/demo.sashp.RData; load(url(link1)); load(url(link2)) # Set up plot background shape and additional text #sa.lst = list(sp.polygons, sa.shp, fill = grey, col = black) # Covers over md.shp #sa.lst = list(sp.polygons, sa.shp, fill = grey, first = F, col = black) # Covers over md.shp #sa.lst = list(sp.polygons, sa.shp, fill = grey, first = T, col = black) # sa.shp doesn't plot sa.lst = list(sp.polygons, sa.shp, col = black) # Shape of South Africa as background polygon - # works, but no fill # Additional text - overall mean yield. text1 - list(sp.text, c(36, -321), paste(Overall mean yield =, round(mean(md@data[-c(22,69), 4]), 0), kg/ha), col =red, cex = 1) # Scale bar limits and title scmax - round(max(md@data[, 4]) / 50, 0) * 50 scmin - round(min(md@data[, 4]) / 50, 0) * 50 byval - 25 titlevec - Mean district maize yields (kg/ha) 1981-99: Cultivar PAN6528 # Spplot dev.new() spplot(md.shp[-c(22,69),], z = 4, sp.layout = list(sa.lst, text1), col.regions = topo.colors(500), at = seq(scmin, scmax, by = byval), main = titlevec) Any advice on how to correctly set sa.lst so that sa.shp fills underneat md.shp will be much appreciated. I haven't been able to find a ready answer in the help files, except for this: The order of items in sp.layout matters; objects are drawn in the order they appear. Plot order and prevalence of sp.layout items: for points and lines, sp.layout items are drawn before the points (to allow for grids and polygons); for grids and polygons sp.layout is drawn afterwards (so the item will not be overdrawn by the grid and/or polygon). Although a matter of taste, transparency may help when combining things. But, beyond what I have tried so far, I can't figure out how to write the code so that sa.shp is drawn first. The second question: 2. Is there an option to pass to: text1 - list(sp.text, c(36, -321), paste(Overall mean yield =, round(mean(md@data[-c(22,69), 4]), 0), kg/ha), col =red, cex = 1) so that the text is left-justified on the coordinates c(36, -321)? At the moment it seems that the text is centered at this location. Thanks in advance for any help you can provide. Best, Lyndon ___ R-sig-Geo mailing list R-sig-Geo@stat.math.ethz.ch https://stat.ethz.ch/mailman/listinfo/r-sig-geo
