Re: [R-sig-Geo] memory Usage setting
Edzer, This is also my biggest frustration with R at the moment. I can not load grids that are bigger than 2M pixels and run geostatistics on them (if I run the same analysis in a GIS software such as SAGA GIS, then I do not get this problems). I often get messages such as: Reached total allocation of 1536Mb: see help(memory.size) Error: cannot allocate vector of size 101.1Mb I still did not figure out how to avoid this obstacle. FYI, I just came back from the geocomputation conference in Ireland where I met some Colleagues from the Centre for e-Science in Lancaster (http://e-science.lancs.ac.uk). They are just about to release an R package called MultiR that should be able to significantly speed up R calculations by employing the grid computing facilities. Daniel Grose (http://e-science.lancs.ac.uk/personnel.html#daniel) mentioned that they plan to organize a workshop on MultiR in November 2007. Grose, D. et al., 2006. sabreR: Grid-Enabling the Analysis of MultiProcess Random Effect Response Data in R. In: P. Halfpenny (Editor), Second International Conference on e-Social Science. National Centre for e-Social Science, Manchester, UK, pp. 12. http://epubs.cclrc.ac.uk/work-details?w=36493 Tom Hengl http://spatial-analyst.net -Original Message- From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] On Behalf Of Edzer J. Pebesma Sent: Thursday, September 13, 2007 12:15 PM To: [EMAIL PROTECTED] Cc: r-sig-geo@stat.math.ethz.ch Subject: Re: [R-sig-Geo] memory Usage setting I think R will never do it's own memory swapping, as that is a typical OS task. There are however several developments (provided in add-on packages) that will not load all data in memory at start-up, but instead call some data base whenever a data element is needed. You might search r-help for rsqlite or biglm, and there are others; also look at the award winners at useR this year. Here, we've run pretty successful R sessions needing 10-11 Gb of memory on a 8Gb RAM 64 bits linux machine with lots of swap space. Needs some patience, and still R might crash other parts of the system when memory usage becomes too excessive. Best regards, -- Edzer Didier Leibovici wrote: Thanks Roger I feel we've got a low RAM machine which would need a bit of an uplift (recent server though)! The linux machine is unfortunately also with 4Gb of RAM But I persist to say it would be interesting to have within R a way of automatically performing swapping memory if needed ... Didier Roger Bivand wrote: On Tue, 11 Sep 2007, [EMAIL PROTECTED] wrote: These days in GIS on may have to manipulate big datasets or arrays. Here I am on WINDOWS I have a 4Gb my aim was to have an array of dim 298249 12 10 22 but that's 2.9Gb Assuming double precision (no single precision in R), 5.8Gb. It used to be (maybe still is?) the case that a single process could only 'claim' a chunk of max size 2GB on Windows. Also remember to compute overhead for R objects... 58 bytes per object, I think it is. It is also strange that once a dd needed 300.4Mb and then 600.7Mb (?) as also I made some room in removing ZZ? Approximately double size - many things the interpreter does involve making an additional copy of the data and then working with *that*. This might be happening here, though I didn't read your code carefully enough to be able to be certain. which I don't really know if it took into account as the limit is greater than the physical RAM of 4GB. ...? :) would it be easier using Linux ? possibly a little bit - on a linux machine you can at least run a PAE kernel (giving you a lot more address space to work with) and have the ability to turn on a bit more virtual memory. usually with data of the size you're trying to work with, i try to find a way to preprocess the data a bit more before i apply R's tools to it. sometimes we stick it into a database (postgres) and select out the bits we want our inferences to be sourced from. ;) it might be simplest to just hunt up a machine with 8 or 16GB of memory in it, and run those bits of the analysis that really need memory on that machine... Yes, if there is no other way, a 64bit machine with lots of RAM would not be so contrained, but maybe this is a matter of first deciding why doing statistics on that much data is worth the effort? It may be, but just trying to read large amounts of data into memory is perhaps not justified in itself. Can you tile or subset the data, accumulating intermediate results? This is the approach the biglm package takes, and the R/GDAL interface also supports subsetting from an external file. Depending on the input format of the data, you should be able to do all you need provided that you do not try to keep all the data in memory. Using a database may be a good idea, or if the data are multiple remote sensing images, subsetting
Re: [R-sig-Geo] Spdep: help needed calculating Moran's I
You can also fit a variogram to the residuals using the gstat package and then record the nugget and sill variation (see http://www.gstat.org/manual/node7.html). I do not know how to test that the nugget variation is statistically significant from the sill variation. You could split the variogram cloud into closer points and more distant points and then run some statistical tests on the two sub-sets (e.g. the t-test). But where to put a boundary between closer/more distant points? Tom Hengl http://spatial-analyst.net -Original Message- From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] On Behalf Of Geertje Van der Heijden Sent: vrijdag 26 oktober 2007 18:12 To: r-sig-geo@stat.math.ethz.ch Subject: [R-sig-Geo] Spdep: help needed calculating Moran's I Hi, I have just posted the same question on the general R help mailing list, but thought that this list might be more appropriate. I am a new user of R. Here is my problem: I have 58 sites from across South America. I have done a regression analysis to relate environmental and biogeographical variables to species richness and want to test whether my residuals are autocorrelated. As far as I understand the Moran's I, I have to take all possible combinations between all points into account to test this. So I have used dnearneigh() with the lower boundary set to 0 and the upper boundary set arbitrarily high to make sure all connections are included. coords - as.matrix(cbind(lowland$long, lowland$lat)) coord.nb - dnearneigh(coords, 0, 1, longlat=TRUE) coord.list - nb2listw(coord.nb, style=W) lianasp.lm - lm(lianasprich ~ log(averdist) + dsl + lianadens + wooddens) lm.morantest(lianasp.lm, coord.list, alternative=two.sided) However, this gives me a Moran's I which is exactly the same as the expected Moran's I (and hence a p-value of 1). If I change the lower or upper boundary slightly so that not all possible links are taken into account, the value is different, but still really near to the expected Moran's I. I don't understand why these values are or the same or nearly so. I am new to spatial statistics, so this might me a really basic question and my appologies if it is, but I am generally a bit at a loss now about the Moran's I and I am wondering if I have calculated it right. Have used to right method to convert my coordinates into neighbourhood distances (and if not, which method should I have used) and am I understanding and calculation the Moran's I correctly? Any help would be greatly appreciated. Many thanks, Geertje Geertje van der Heijden PhD student Tropical Ecology School of Geography University of Leeds Leeds LS2 9JT Tel: (+44)(0)113 3433345 Email: [EMAIL PROTECTED] [[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
Re: [R-sig-Geo] regression kriging in gstat with skewed distributions
Dear Giovanni, Logit transformation can be automatically applied to any variables which has a lower and upper physical limits (e.g. 0-100%). In R, you can transform a variable to logits by e.g.: points = read.dbf(points.dbf) points$SANDt = log((points$SAND/100)/(1-(points$SAND/100))) After you interpolate your variable, you can back-transform the values by using: SAND.rk = krige(fsand$call$formula, points[sel,], SPC, sand.rvgm) SAND.rk$pred=exp(SAND.rk$var1.pred)/(1+exp(SAND.rk$var1.pred))*100 The prediction variance can not be back-transformed, but you can use the normalized prediction variance by dividing it with the sampled variance. See also section 4.2.1 of my lecture notes (http://geostat.pedometrics.org/). There are many transformations that can be applied to force a normality of your target variable (see e.g. http://en.wikipedia.org/wiki/Data_transformation_(statistics) ). The most generic transformation is to work with the probability density function values (see e.g. http://dx.doi.org/10.1016/j.jneumeth.2006.11.004 ), this way you do not have to think about how the histogram looks at all. But then the interpretation of the regression plots becomes rather difficult. In any case, you should apply the transformation already to the target variable because also a requirement for linear regression is that the residuals are normally distributed around the regression line. see also: FITTING DISTRIBUTIONS WITH R (by Vito Ricci) http://cran.r-project.org/doc/contrib/Ricci-distributions-en.pdf Tom Hengl http://spatial-analyst.net -Original Message- From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] On Behalf Of G. Allegri Sent: dinsdag 15 januari 2008 15:28 To: r-sig-geo@stat.math.ethz.ch Subject: [R-sig-Geo] regression kriging in gstat with skewed distributions I'm trying to realize e regression kriging with gstat package on my soil samples data. The response variable (ECe measuere) and covariates appear positvely skewed. As Tomislav Hengl suggests in its framework for RK [1], a logistic transformation is proposed as a generic way to reduce the skeweness by using the physical limits of the data. Is it really a transformation that can be applied in the generic case of skewed datas? I mean,in my case I have non-normal residuals (from original data regression), and I'm trying to transform the residuals (and not the original values) to do SK on them . Is this approach correct? A related question is how to do normal score transformations (for my residuals) in R and gstat. I know gstat doesn't manage transformations and back-transformations, so it should be done previously in R... but I can't find any package that permit it in a straisghtforward way. I've found something with qqnorm(ppoints(data)) and the approx() function. Is that all? Giovanni [1] A generic framework for spatial prediction of soil variables based on regressionkriging Geoderma 122 (1–2), 75–93. ___ 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] kernel density estimation to google earth
# Modified by T. Hengl (http://spatial-analyst.net) # 27.01.2008 library(rgdal) library(maptools) library(splancs) # Import the points and study area: data.shp - readOGR(C:/, layer=events) str(data.shp) poly.shp - readOGR(C:/, layer=hull) str(poly.shp) poly - getPolygonCoordsSlot(getPolygonsPolygonsSlot(getSpPpolygonsSlot(poly.shp)[[1]])[[1]]) grd - GridTopology(cellcentre.offset=c(round([EMAIL PROTECTED]r1,min],0), round([EMAIL PROTECTED]r2,min],0)), cellsize=c(2000,2000), cells.dim=c(30,25)) # Run the 2D kernel smoother: kbw2000 - spkernel2d(data.shp, poly, h0=2000, grd) hist(kbw2000) # Pack and plot a SpatialGridDataFrame: kbw2000.grd - SpatialGridDataFrame(grd, data=data.frame(kbw2000)) proj4string(kbw2000.grd) - [EMAIL PROTECTED] spplot(kbw2000.grd, col.regions=terrain.colors(16), sp.layout=list(sp.points,pch=+,cex=1.2,col=black,data.shp)) # Export to KML # First, reproject the grid to longlat: kbw2000.ll - spTransform(kbw2000.grd, CRS(+proj=longlat +datum=WGS84)) str(kbw2000.ll) # The cell size you need to determine yourself!! width = ([EMAIL PROTECTED]coords.x1,max[EMAIL PROTECTED]coords.x1,min])/2000 height = ([EMAIL PROTECTED]coords.x2,max[EMAIL PROTECTED]coords.x2,min])/2000 geogrd.cell = ([EMAIL PROTECTED]s1,max[EMAIL PROTECTED]s1,min])/width # Define a new grid: geogrd = spsample(kbw2000.ll, type=regular, cellsize=c(geogrd.cell,geogrd.cell)) gridded(geogrd) = TRUE gridparameters(geogrd) # cellcentre.offset cellsize cells.dim # x1 15.90165 0.0263668530 # x2 47.95541 0.0263668516 # This is an empty grid without any topology (only grid nodes are defined) and coordinate # system definition. To create topology, we coerce a dummy variable (1s), then # specify that the layer has a full topology: nogrids = [EMAIL PROTECTED]@cells.dim[x1[EMAIL PROTECTED]@cells.dim[x2] geogrd = SpatialGridDataFrame([EMAIL PROTECTED], data=data.frame(rep(1, nogrids)), [EMAIL PROTECTED]) # and estimate the values of the reprojected map at new grid locations using the bilinear resampling: # this can be time-consuming for large grids!!! library(gstat) kbw2000.llgrd = krige(kbw2000~1, kbw2000.ll, geogrd, nmax=4) # Optional, convert the original shape to latlong coordinates: data.ll - spTransform(data.shp, CRS(+proj=longlat +datum=WGS84)) spplot(kbw2000.llgrd[var1.pred], col.regions=terrain.colors(16), scales=list(draw=TRUE), sp.layout=list(sp.points,pch=+,cex=1.2,col=black,data.ll)) # The final grid map can be exported to KML format using the maptools package and kmlOverlay method: kbw2000.kml = GE_SpatialGrid(kbw2000.llgrd) tf - tempfile() png(file=paste(tf, .png, sep=), width=kbw2000.kml$width, height=kbw2000.kml$height, bg=transparent) par(mar=c(0,0,0,0), xaxs=i, yaxs=i) image(as.image.SpatialGridDataFrame(kbw2000.llgrd[1]), col=bpy.colors(), xlim=kbw2000.kml$xlim, ylim=kbw2000.kml$ylim) plot(data.ll, pch=+, cex=1.2, add=TRUE, bg=transparent) kmlOverlay(kbw2000.kml, paste(tf, .kml, sep=), paste(tf, .png, sep=)) dev.off() # see also: # Hengl, T., 2007. A Practical Guide to Geostatistical Mapping of Environmental Variables. EUR 22904 EN Scientific and Technical Research series, Office for Official Publications of the European Communities, Luxemburg, 143 pp. Dear List, I've calculated a kernel density estimation (splancs function) and now I want to export the result as a kml-file to open it in the google earth viewer, but I get stuck at converting the SpatialGridDataFrame to a SpatialPixelDataFrame... here (http://www.zshare.net/download/689742375ef5fb/) you can download some testdata and my R-code so far: ### data.shp - readOGR(C:/, layer=events) prj - data.shp@ proj4string@ projargs dat - data.shp str(dat) poly.shp - readOGR(C:/, layer=hull) str(poly.shp) dat.SP - as(dat, SpatialPoints) str(dat.SP) pp_poi - as.points([EMAIL PROTECTED],1], [EMAIL PROTECTED],2]) poly - getPolygonCoordsSlot(getPolygonsPolygonsSlot(getSpPpolygonsSlot(poly.shp)[[1]])[[1]]) polymap(poly) points(pp_poi) grd - GridTopology(cellcentre.offset=c(590511, 396191), cellsize=c(2000, 2000), cells.dim=c(30,25)) kbw2000 - spkernel2d(pp_poi, poly, h0=2000, grd) spplot(SpatialGridDataFrame(grd, data=data.frame(kbw2000)), col.regions=terrain.colors(16)) test - SpatialGridDataFrame(grd, data=data.frame(kbw2000)) proj4string(test) - CRS(prj) str(test) test1 - spTransform(test, CRS(+proj=longlat +datum=WGS84)) ## here I got following error message: # validityMethod(as(object, superClass)): Geographical CRS given to non-conformant data test2 - spsample(test1, type=regular, cellsize=c(2000,2000)) # export the SpatialPixelDataFrame (Code (not testet) from the help-file) tf - tempfile() SGxx - GE_SpatialGrid(test2) png(file=paste(tf, .png, sep=), width=SGxx$width, height=SGxx$height, bg=transparent) par(mar=c(0,0,0,0), xaxs=i, yaxs=i) plot(x, xlim=SGxx$xlim, ylim=SGxx$ylim, setParUsrBB=TRUE) dev.off() kmlOverlay(SGxx, paste(tf, .kml,
Re: [R-sig-Geo] R/gstat UK with indicators
David, Thanks for you note. It is a very clear example and I have experienced this problem many times. If I can be of any assistance, here are few discussion points that might put some light on the whole thing. Gstat does not support any stat models different than linear. This is mainly because gstat implements the so-called Kriging with External Drift algorithm (covariances of residuals extended with auxiliary predictors). This algorithm is rather elegant and gives then equivalent results as if you would first estimate linear regression by GLS and then interpolate and sum the residuals. But the drawback of KED is that it only works with linear models (multiple linear regression). As far as I know, KED has been implemented in all software that can run universal kriging interpolation with auxiliary predictors (maps). The alternative is to run the so-called Regression-kriging approach where you separately model the deterministic and stochastic part (as you have also demonstrated; see also sec 2.1 of my lecture notes). However, I think that there is a complete theory missing (maybe you should report on this?)! For example, if you apply binomial link function, GAMs or regression trees, I have a feeling that you should also consider the spatio-correlation during the estimation of the model parameters. In the case of multiple regression, covariance matrix is used to account for spreading (clustering) of the points in the space. In your example, you fit a GLM that completely ignores location of the points, so I have a feeling that it is not statistically optimal. We have also reported recently on interpolation of soil categories (http://dx.doi.org/10.1016/j.geoderma.2007.04.022). We discovered that the easiest thing to do is to actually have memberships instead of indicators (it is also easier to fit variograms if you work with memberships). Then the system is equivalent to regression-kriging (after logit transformation), so then you do not have to worry about the link function (see also section 4.3.3 and Fig. 4.15 of my lecture notes). All the best, Tom Hengl http://spatial-analyst.net Hengl, T., 2007. A Practical Guide to Geostatistical Mapping of Environmental Variables. EUR 22904 EN Scientific and Technical Research series, Office for Official Publications of the European Communities, Luxemburg, 143 pp. http://bookshop.europa.eu/uri?target=EUB:NOTICE:LBNA22904:EN:HTML I am sure Edzer can answer this, but I post it to the list because it may be of general interest. If I compute an indicator: data(meuse) meuse$i.Pb.med - meuse$lead median(meuse$lead) I can of course compute an indicator variogram and krige at new locations: vi - variogram(i.Pb.med ~1, meuse) plot(vi, pl=T) (vim - fit.variogram(vi, vgm(0.17, Sph, 600, 0.1))) model psillrange 1 Nug 0.08772438 0. 2 Sph 0.18157082 691.2488 plot(vi, pl=T, model=vim) ki - krige(i.Pb.med ~1, meuse, newdata=meuse.grid, model=vm) [using ordinary kriging] summary(ki$var1.pred) Min. 1st Qu. MedianMean 3rd Qu.Max. 0.03628 0.33830 0.64020 0.58080 0.81610 0.95040 I can also fit a logistic regression model to the indicator with a GLM: (gm - glm(i.Pb.med ~ dist, family=binomial, [EMAIL PROTECTED])) (Intercept) dist -2.562 11.662 Degrees of Freedom: 154 Total (i.e. Null); 153 Residual Null Deviance: 214.9 Residual Deviance: 131.2AIC: 135.2 Here the ~ in the formula is interpreted according to the link function. I can compute residuals from the GLM, get the residual variogram, and model it: meuse$glm.r - residuals(gm) v - variogram(glm.r ~ 1, meuse, cutoff=1000) plot(v, pl=T) vm - vgm(0.45, Exp, 450, 0.6) # eyeball fit But, what happens when I use the same formula in gstat? First, just to compute the residual variogram: vr - variogram(i.Pb.med ~ dist, meuse, cutoff=1000) The two residual variograms are very different: max(v$gamma) [1] 1.001015 max(vr$gamma) [1] 0.1741679 So clearly the formula (i.Pb.med ~ dist) is being interpreted as a linear model, not a logit model. Second, an attempt at kriging with external drift: ked - krige(i.Pb.med ~ dist, loc=meuse, newdata=meuse.grid, model=vm) [using universal kriging] Since I could not specify a link function, it seems that a linear model is being fit (not a logit model). And we see very strange results for an indicator: summary(ked$var1.pred) Min. 1st Qu. MedianMean 3rd Qu.Max. -0.1235 0.2719 0.6550 0.6056 0.9139 1.4990 So (after all that buildup), I just want to know if my interpretation is correct, that gstat is not interpreting the formula in krige() or in variogram() as a GLM, rather as a linear model, which doesn't make sense with indicators. Thank you, D G Rossiter Senior University Lecturer Department of Earth Systems Analysis (DESA) International Institute for Geo-Information Science and Earth Observation (ITC) Enschede, The Netherlands
Re: [R-sig-Geo] RGB Image Plotting!
If you import the bands/images as integers [0-255], then you can visualize them using e.g.: library(colorspace) library(rgdal) vismaps = readGDAL(band_R.tif) vismaps$red = vismaps$band1 vismaps$green = readGDAL(band_G.tif)$band1 vismaps$blue = readGDAL(band_B.tif)$band1 # Display as a RGB image: RGBimg = SGDF2PCT(vismaps[c(red, green, blue)], ncolors=256, adjust.bands=FALSE) vismaps$idx - RGBimg$idx image(vismaps, idx, col=RGBimg$ct) see also http://spatial-analyst.net/VMM/R_whitening.zip On Mon, 25 Feb 2008, PUJAN RAJ REGMI wrote: Dear Madam/ Sir I am New user of R. I am looking for the way to plot true color image from matrix data of of order (m,n,3) as in MAT LAB where 3 represent the three column matrix where red, green and blue values of pixels are stored. I hope you understand my question. Looking for your help for the same. Thanking you in advance Yours, Please look at the help page and example for RGB2PCT() in the rgdal package. That should get you most of the way, using other functions in the package if need be to convert your input matrix to a three-band GDALReadOnlyDataset. Hope this helps, Roger Pujan Raj Regmi.Abroad Address: Studentenwijk Arenberg 22/408, 3001 Heverlee, Belgium. Cell Number: 0032-0486-530-340 Home Address: Prabachan Marg-182/15, Baneshwor, Kathmandu, Nepal. Home Telephone Number: +977-1-4470900Alternative e-mail: [EMAIL PROTECTED] [EMAIL PROTECTED] _ Helping your favorite cause is as easy as instant messaging. You IM, we give. [[alternative HTML version deleted]] -- Roger Bivand Economic Geography Section, Department of Economics, Norwegian School of Economics and Business Administration, Helleveien 30, N-5045 Bergen, Norway. voice: +47 55 95 93 55; fax +47 55 95 95 43 e-mail: [EMAIL PROTECTED] ___ 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] Applied Spatial Data Analysis with R from Springer
Dear Roger, Edzer and Rubio, I would like to join the congratulations (hopefully we will not overload the mailing list :))) ). I have been already recommending it to many of my students (I saw the draft at Edzer's place few months ago). BTW, the book's website is not available. Is the URL correct? http://www.asdar-book.org all the best, PS: Congratulations on the affordable price also! Tom Hengl http://spatial-analyst.net Thanks for such good news. I'm really looking forward to reading it. Thanks a lot Edzer Pebesma, Virgilio Gómez-Rubio and Roger Bivand for your contribution to the universal knowledge. -- hzambran Linux user #454569 -- Ubuntu user #17469 ___ 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] question about regression kriging
I completely agree with Thierry.
Take a look at this also:
https://stat.ethz.ch/pipermail/r-sig-geo/2008-February/003176.html
The instructions on how to run RK with binary variables in R you can find in
sec 4.3.3 (Fig. 4.15)
of my lecture notes.
Hengl, T., 2007. A Practical Guide to Geostatistical Mapping of Environmental
Variables. EUR 22904
EN Scientific and Technical Research series, Office for Official Publications
of the European
Communities, Luxemburg, 143 pp.
http://bookshop.europa.eu/uri?target=EUB:NOTICE:LBNA22904:EN:HTML
Tom Hengl
http://spatial-analyst.net
-Original Message-
From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] On Behalf Of
ONKELINX, Thierry
Sent: woensdag 2 april 2008 11:56
To: Vanessa Stelzenmuller (Cefas); r-sig-geo@stat.math.ethz.ch
Subject: Re: [R-sig-Geo] question about regression kriging
Dear Vanessa,
What residuals did you use? The ones in the original scale or in the logit
scale? Interpolate the
residuals in the logit scale and add these to the model predictions in the
logit scale. And the
transform those values back to the original scale. This will prevent values
outside the 0-1 range.
Maybe you should have a loot at the geoRglm package.
HTH,
Thierry
ir. Thierry Onkelinx
Instituut voor natuur- en bosonderzoek / Research Institute for Nature and
Forest
Cel biometrie, methodologie en kwaliteitszorg / Section biometrics, methodology
and quality
assurance
Gaverstraat 4
9500 Geraardsbergen
Belgium
tel. + 32 54/436 185
[EMAIL PROTECTED]
www.inbo.be
To call in the statistician after the experiment is done may be no more than
asking him to perform a
post-mortem examination: he may be able to say what the experiment died of.
~ Sir Ronald Aylmer Fisher
The plural of anecdote is not data.
~ Roger Brinner
The combination of some data and an aching desire for an answer does not ensure
that a reasonable
answer can be extracted from a given body of data.
~ John Tukey
-Oorspronkelijk bericht-
Van: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] Namens Vanessa
Stelzenmuller (Cefas)
Verzonden: woensdag 2 april 2008 11:32
Aan: r-sig-geo@stat.math.ethz.ch
Onderwerp: [R-sig-Geo] question about regression kriging
Hello,
We work on the application of regression kriging to presence / absence data in
the context of
species distribution modelling. In R in a first step we fit the trend surfaces
with logistic
regression models. Then we fit a variogram to the regression residuals and
interpolate the residuals
with OK. Now we face the situation that when combining trend surfaces with
residual surfaces for
some locations our occurrence probability is 0 or 1. Thus taking into account
the spatial
structure of the data (residuals) has the potential to convert a predicted
high occurrence
probability into a low occurrence probability or vice versa. Are there some
restriction for
presence/ absence data for this approach? How to deal with these estimations
(0 and 1)?
Many thanks
Vanessa
Dr. Vanessa Stelzenmüller
Marine Scientist (GIS), CEFAS
Pakefield Road, Lowestoft, NR33 0HT, UK
Tel.: +44 (0)1502 527779
www.cefas.co.uk
***
This email and any attachments are intended for the =\...{{dropped:8}}
___
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] question about regression kriging
The two components of the regression-kriging model are not independent, hence you are doing a wrong thing if you are just summing them. You should use instead the universal kriging variance that is derived in gstat. The complete derivation of the Universal kriging variance is available in Cressie (1993; p.154), or even better Papritz and Stein (1999; p.94). See also pages 7-8 of our technical note: Hengl T., Heuvelink G.B.M. and Stein A., 2003. Comparison of kriging with external drift and regression-kriging. Technical report, International Institute for Geo-information Science and Earth Observation (ITC), Enschede, pp. 18. http://www.itc.nl/library/Papers_2003/misca/hengl_comparison.pdf Edzer is right, you can not back-transform prediction variance of the transformed variable (logits). However, you can standardize/normalize the UK variance by diving it with global variance (see e.g. http://dx.doi.org/10.1016/j.geoderma.2003.08.018), so that you can evaluate the success of prediction in relative terms (see also http://spatial-analyst.net/visualization.php). Tom Hengl http://spatial-analyst.net -Original Message- From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] On Behalf Of Edzer Pebesma Sent: dinsdag 8 april 2008 20:50 To: David Maxwell (Cefas) Cc: r-sig-geo@stat.math.ethz.ch Subject: Re: [R-sig-Geo] question about regression kriging David Maxwell (Cefas) wrote: Hi, Tom and Thierry, Thank you for your advice, the lecture notes are very useful. We will try geoRglm but for now regression kriging using the working residuals gives sensible answers even though there are some issues with using working residuals, i.e. not Normally distributed, occasional very large values and inv.logit(prediction type=link + working residual) doesn't quite give the observed values. Our final question about this is how to estimate standard errors for the regression kriging predictions of the binary variable? On the logit scale we are using rk prediction (s0) = glm prediction (s0) + kriged residual prediction (s0) for location s0 Is assuming independence of the two components adequate? var rk(s0) ~= var glm prediction (s0) + var kriged residual prediction (s0) In principle, no. The extreme case is prediction at observation locations, where the correlation is -1 so that the final prediction variance becomes zero. I never got to looking how large the correlation is otherwise, but that shouldn't be hard to do in the linear case, as you can get the first and second separately, and also the combined using universal kriging. Another question: how do you transform this variance back to the observation scale? -- Edzer ___ 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] question about regression kriging
Sorry, I meant UK variance of the interpolated logits. I guess the solution is to either: (A) 1. Run 100 SG simulations using logits; 2. Back-transform the values to 0-1 scale; 3. Determine the variance per grid node (now in 0-1 scale). Running simulations could also be rather time-consuming (I often give up running SG simulations in gstat with multiple raster maps as predictors -- memory limit problems, too time-consuming). (B) 1. Determine the upper and lower 1-s confidence intervals -- 2 maps (in logit scale); 2. Back-transform the maps and derive the difference in 0-1 scale (estimation error); But I guess that a method to directly back-transform the UK variance from logit scale to 0-1 scale does not exist. Tom Hengl -Original Message- From: Edzer Pebesma [mailto:[EMAIL PROTECTED] Sent: woensdag 9 april 2008 14:48 To: Tomislav Hengl Cc: r-sig-geo@stat.math.ethz.ch Subject: Re: [R-sig-Geo] question about regression kriging Tomislav Hengl wrote: PS: How do you back-transform the GLM prediction variance to original scale (I was not aware of this, apologies)? A reference would do. 1. The white book, Statistical Models in S 2. The source of predict.glm (not for the faint of heart) 3. McCullogh Nelder (I guess it's in there, but my copy is still in a box!) -- Edzer ___ 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] Problem with universal kriging using gstat
Auxiliary variables that are used to explain the trend-part of variation need to be available also at all new prediction locations. see ?krige: newdata - data frame or Spatial object with prediction/simulation locations; should contain attribute columns with the independent variables (if present) and (if locations is a formula) the coordinates with names as defined in locations Gstat obviously can not find X and Y coordinates of your new locations (data.grid). Try instead: data.grid - expand.grid(x=seq(xLL,xUL,xN), y=seq(xLL,xUL,xN)) names(data.grid) = c(X,Y) gridded(data.grid) - ~X+Y cheers, Tom Hengl http://spatial-analyst.net/ -Original Message- From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] On Behalf Of Yong Li Sent: maandag 14 april 2008 6:45 To: r-sig-geo@stat.math.ethz.ch Subject: [R-sig-Geo] Problem with universal kriging using gstat Hi ALL, Can any expert to see my mistake in the following R script when I try the gstat universal kriging? The error occurs at the last step. It sounds I missed the setup for the names of coordinates of the SpatialGrid data.grid. Regards Yong memory.size() [1] 47930616 memory.size(TRUE) [1] 85377024 memory.limit(size=2048) NULL round(memory.limit()/1048576.0, 2) [1] 2048 options(scipen=3) options(digits=15) rm(list=ls()) require(sp) [1] TRUE require(gstat) [1] TRUE require(maptools) [1] TRUE require(foreign) [1] TRUE cellsize - 10 # 1 2 3 4 VarName=c(logit(SOM),logit(TN),logit(OLSENP),logit(EXTK)) out=paste(Variable,No,Model,C0,C,Spatial,Range,R2,p,SSE,sep=,) setwd(E:\\Rcode\\EJSS\\Hongtong\\village) data - read.dbf(libaocun_gs.dbf) names(data) [1] ID XCOORYCOORNO CDBH NO3N NH4N MINERALN SOM TOTALN OLSENP EXTK [13] TN XY class(data) [1] data.frame names(data)[14] - X names(data)[15] - Y data$X - data$X + 1750 data$Y - data$Y + 650 coordinates(data) - ~X+Y proj4string(data) - CRS(+proj=tmerc +lat_0=0 +lon_0=111 +k=1.0 +x_0=50 +y_0=0 +ellps=krass +units=m +no_defs +towgs84=22,-118,30.5,0,0,0,0) #read in a shape file of the boundary aSHAPE - readShapePoly(libao_bd.shp, IDvar=BSM, proj4string=CRS(proj4string(data))) #generate the grid coordinates (in meters) xLL - round(slot(aSHAPE,bbox)[1]) - 100 xUR - round(slot(aSHAPE,bbox)[3]) + 100 yLL - round(slot(aSHAPE,bbox)[2]) - 100 yUR - round(slot(aSHAPE,bbox)[4]) + 100 xN - round((xUR-xLL)/cellsize) yN - round((yUR-yLL)/cellsize) data.grid = SpatialGrid(GridTopology(c(xLL,yLL),c(cellsize,cellsize),c(xN,yN))) proj4string(data.grid) - CRS(proj4string(aSHAPE)) #data logit transformation SOM_MAX - max(data$SOM)*1.1 SOM_MIN - min(data$SOM)*0.9 SOM_DIFF - SOM_MAX - SOM_MIN data$SOMt- log(((data$SOM-SOM_MIN)/SOM_DIFF)/(1-(data$SOM-SOM_MIN)/SOM_DIFF)) TN_MAX - max(data$TN)*1.1 TN_MIN - min(data$TN)*0.9 TN_DIFF - TN_MAX - TN_MIN data$TNt- log(((data$TN-TN_MIN)/TN_DIFF)/(1-(data$TN-TN_MIN)/TN_DIFF)) OLSENP_MAX - max(data$OLSENP)*1.1 OLSENP_MIN - min(data$OLSENP)*0.9 OLSENP_DIFF - OLSENP_MAX - OLSENP_MIN data$OLSENPt- log(((data$OLSENP-OLSENP_MIN)/OLSENP_DIFF)/(1-(data$OLSENP-OLSENP_MIN)/OLSENP_DIFF)) EXTK_MAX - max(data$EXTK)*1.1 EXTK_MIN - min(data$EXTK)*0.9 EXTK_DIFF - EXTK_MAX - EXTK_MIN data$EXTKt- log(((data$EXTK-EXTK_MIN)/EXTK_DIFF)/(1-(data$EXTK-EXTK_MIN)/EXTK_DIFF)) #SOM # plot X11(width=8, height=6.5) par(mfrow=c(1,1)) spplot(aSHAPE[BSM], scales=list(draw=TRUE, cex=0.7), xlab=Easting (m), ylab=Northing (m), + sp.layout=list(sp.points, pch=+, col=black, data), col.regions=FALSE, colorkey=FALSE) X11(width=8, height=6.5) par(mfrow=c(1,2)) hist(data$SOMt, main=, xlab=VarName[1], n=25) boxplot(data$SOMt, xlab=VarName[1]) #estimate semivariogram model soil.v - variogram(SOMt~X+Y+X*Y+X^2+Y^2, data=data, cutoff=1050, width=90) vmodel - vgm(0.5, Sph, 1000, 0.1) v.sph - fit.variogram(object=soil.v, model=vmodel, fit.sills=TRUE, fit.ranges=TRUE, fit.method=7) v.sph model psillrange 1 Nug 0.269063860535796 0. 2 Sph 0.217370623520010 815.994438576304 X11(width=8, height=6.5) par(mfrow=c(1,1)) plot(soil.v, model=v.sph, col=red, cex=1, lwd=2.0, main=VarName[1], xlab=Separation distance (m)) #Univeral kriging data.gstat - gstat(id=SOMt, formula=SOMt~X+Y+X*Y+X^2+Y^2, data=data, nmax=15, model=v.sph) #Predicting surface data.grid - predict(data.gstat, data.grid) Error in eval(expr, envir, enclos) : object X not found [[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
Re: [R-sig-Geo] Problem with universal kriging using gstat
sorry, a type error: cellsize - 10 # see also http://spatial-analyst.net/pixel.php xLL - round(slot(aSHAPE,bbox)[1])+cellsize/2 xUR - round(slot(aSHAPE,bbox)[3])-cellsize/2 yLL - round(slot(aSHAPE,bbox)[2])+cellsize/2 yUR - round(slot(aSHAPE,bbox)[4])-cellsize/2 xN - round((xUR-xLL)/cellsize) yN - round((yUR-yLL)/cellsize) data.grid - expand.grid(x=seq(xLL,xUR,xN), y=seq(yLL,yUR,yN)) names(data.grid) = c(X,Y) gridded(data.grid) - ~X+Y -Original Message- From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] On Behalf Of Tomislav Hengl Sent: maandag 14 april 2008 10:37 To: 'Yong Li'; r-sig-geo@stat.math.ethz.ch Subject: Re: [R-sig-Geo] Problem with universal kriging using gstat Auxiliary variables that are used to explain the trend-part of variation need to be available also at all new prediction locations. see ?krige: newdata - data frame or Spatial object with prediction/simulation locations; should contain attribute columns with the independent variables (if present) and (if locations is a formula) the coordinates with names as defined in locations Gstat obviously can not find X and Y coordinates of your new locations (data.grid). Try instead: data.grid - expand.grid(x=seq(xLL,xUL,xN), y=seq(xLL,xUL,xN)) names(data.grid) = c(X,Y) gridded(data.grid) - ~X+Y cheers, Tom Hengl http://spatial-analyst.net/ -Original Message- From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] On Behalf Of Yong Li Sent: maandag 14 april 2008 6:45 To: r-sig-geo@stat.math.ethz.ch Subject: [R-sig-Geo] Problem with universal kriging using gstat Hi ALL, Can any expert to see my mistake in the following R script when I try the gstat universal kriging? The error occurs at the last step. It sounds I missed the setup for the names of coordinates of the SpatialGrid data.grid. Regards Yong memory.size() [1] 47930616 memory.size(TRUE) [1] 85377024 memory.limit(size=2048) NULL round(memory.limit()/1048576.0, 2) [1] 2048 options(scipen=3) options(digits=15) rm(list=ls()) require(sp) [1] TRUE require(gstat) [1] TRUE require(maptools) [1] TRUE require(foreign) [1] TRUE cellsize - 10 # 1 2 3 4 VarName=c(logit(SOM),logit(TN),logit(OLSENP),logit(EXTK)) out=paste(Variable,No,Model,C0,C,Spatial,Range,R2,p,SSE,sep=,) setwd(E:\\Rcode\\EJSS\\Hongtong\\village) data - read.dbf(libaocun_gs.dbf) names(data) [1] ID XCOORYCOORNO CDBH NO3N NH4N MINERALN SOM TOTALN OLSENP EXTK [13] TN XY class(data) [1] data.frame names(data)[14] - X names(data)[15] - Y data$X - data$X + 1750 data$Y - data$Y + 650 coordinates(data) - ~X+Y proj4string(data) - CRS(+proj=tmerc +lat_0=0 +lon_0=111 +k=1.0 +x_0=50 +y_0=0 +ellps=krass +units=m +no_defs +towgs84=22,-118,30.5,0,0,0,0) #read in a shape file of the boundary aSHAPE - readShapePoly(libao_bd.shp, IDvar=BSM, proj4string=CRS(proj4string(data))) #generate the grid coordinates (in meters) xLL - round(slot(aSHAPE,bbox)[1]) - 100 xUR - round(slot(aSHAPE,bbox)[3]) + 100 yLL - round(slot(aSHAPE,bbox)[2]) - 100 yUR - round(slot(aSHAPE,bbox)[4]) + 100 xN - round((xUR-xLL)/cellsize) yN - round((yUR-yLL)/cellsize) data.grid = SpatialGrid(GridTopology(c(xLL,yLL),c(cellsize,cellsize),c(xN,yN))) proj4string(data.grid) - CRS(proj4string(aSHAPE)) #data logit transformation SOM_MAX - max(data$SOM)*1.1 SOM_MIN - min(data$SOM)*0.9 SOM_DIFF - SOM_MAX - SOM_MIN data$SOMt- log(((data$SOM-SOM_MIN)/SOM_DIFF)/(1-(data$SOM-SOM_MIN)/SOM_DIFF)) TN_MAX - max(data$TN)*1.1 TN_MIN - min(data$TN)*0.9 TN_DIFF - TN_MAX - TN_MIN data$TNt- log(((data$TN-TN_MIN)/TN_DIFF)/(1-(data$TN-TN_MIN)/TN_DIFF)) OLSENP_MAX - max(data$OLSENP)*1.1 OLSENP_MIN - min(data$OLSENP)*0.9 OLSENP_DIFF - OLSENP_MAX - OLSENP_MIN data$OLSENPt- log(((data$OLSENP-OLSENP_MIN)/OLSENP_DIFF)/(1-(data$OLSENP-OLSENP_MIN)/OLSENP_DIFF)) EXTK_MAX - max(data$EXTK)*1.1 EXTK_MIN - min(data$EXTK)*0.9 EXTK_DIFF - EXTK_MAX - EXTK_MIN data$EXTKt- log(((data$EXTK-EXTK_MIN)/EXTK_DIFF)/(1-(data$EXTK-EXTK_MIN)/EXTK_DIFF)) #SOM # plot X11(width=8, height=6.5) par(mfrow=c(1,1)) spplot(aSHAPE[BSM], scales=list(draw=TRUE, cex=0.7), xlab=Easting (m), ylab=Northing (m), + sp.layout=list(sp.points, pch=+, col=black, data), col.regions=FALSE, colorkey=FALSE) X11(width=8, height=6.5) par(mfrow=c(1,2)) hist(data$SOMt, main=, xlab=VarName[1], n=25) boxplot(data$SOMt, xlab=VarName[1]) #estimate semivariogram model soil.v - variogram(SOMt~X+Y+X*Y+X^2+Y^2, data=data, cutoff=1050, width=90) vmodel - vgm(0.5, Sph, 1000, 0.1) v.sph - fit.variogram(object=soil.v, model=vmodel, fit.sills=TRUE, fit.ranges=TRUE, fit.method=7) v.sph model psillrange 1 Nug 0.269063860535796 0. 2 Sph
[R-sig-Geo] Memory limit problems in R / import of maps
Dear list, I know that much has already been said about the memory limit problems. If there is any progress about this problem, we would be interested to hear. In our project, we are importing 24 maps/bands, each consists of 1,450,000 pixels. We further would like to glue all maps into a single data frame (e.g. 'ksc' class in adehabitat package; or 'SpatialGridDataFrame' in sp package), but this seems to be impossible. We tried to run this under windows (after following http://cran.r-project.org/bin/windows/base/rw-FAQ.html#There-seems-to-be-a-limit-on-the-memory-it-us es_0021 and setting the --max-mem-size) and under Linux Ubuntu, but still get the same error message (seems that there is no difference in use of memory under the two OS): Error: Cannot allocate vector of size 11.1 Mb The R workspace with 24 loaded grids is also really small (18 MB), but any further gluing and calculation is blocked due the vector size error message. For a comparison, in a GIS such as ArcGIS or SAGA/ILWIS (open source) we have no problems of loading and processing 3-4 times more grids. Should we simply give up on running spatial analysis using large grids (10 million grids) in R? Any suggestion or comment is more than welcome, Tom Hengl http://spatial-analyst.net ___ 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] Memory limit problems in R / import of maps
Dear Edzer, Roger, Thanks for the tips! BTW, I was using Dell Latitude D630 laptop with 2Gb of RAM and 2GHz Intel chip with: (a) Windows XP professional (b) Linux Ubuntu Under both OS I get the same error message. And yes, we would like to import about 25 maps of 1.5 millions pixels each. I will try to summarize your strategies: (1) Import the data in tiles (e.g. via rgdal), e.g.: info - GDALinfo(dem50m.asc) gridmaps01 = readGDAL(dem50m.asc, region.dim = round(c(info[[rows]]/2,info[[columns]]))) gridmaps02 = readGDAL(dem50m.asc, region.dim = round(c(info[[rows]]/2,info[[columns]])), offset=round(c(info[[rows]]/2,0))) Or try reading directly from disk (any documentation on how to achieve this?). (2) Obtain a 64 bits OS/PC with 10Gb of RAM. We would like to use the package adehabitat that reads ArcInfo ASCII maps via 'import.asc'. Unfortunately, I could not find that it can read the data in tiles, but we might try importing via rgdal first and then converting to the 'ksc' class. If anybody else has experience/solution with working with large maps we would be interested to hear/learn from his/hear opinion. Just one last thing, if R is reporting an error message, that does not necessarily mean that there is a memory limit problem with the machine - shouldn't there be a way to implement memory handling in R in a more efficient way? Thanks in any case, Tom Hengl -Original Message- From: Edzer Pebesma [mailto:[EMAIL PROTECTED] Sent: dinsdag 22 april 2008 16:14 To: Tomislav Hengl Cc: r-sig-geo@stat.math.ethz.ch; 'Michalis Vardakis' Subject: Re: [R-sig-Geo] Memory limit problems in R / import of maps Hi Tom, Tomislav Hengl wrote: Should we simply give up on running spatial analysis using large grids (10 million grids) in R? Yes, and I would be very interested to hear along which other path you were then successful to finish the job. Other options I can see are: - buy a decent pc with 16 or 32 Gb memory, and use 64 bits linux (have you checked how much this would cost, and compared it to the budget of your project?). There's nothing special about it, I use it 100% of my time on my 1.2 kg laptop (with much less RAM). OR: - don't go through the grid in a single pass, but do it by tiles, e.g. use rgdal to read part of the grid and do that for 100 tiles, should reduce memory needs with a factor 100. Of course this takes a little bit more effort, in terms of administration (as Roger mentioned), OR: - rewrite the memory-hungry parts such that the bulky data is not first read into memory, but read directly from disk. Several attempts can be found in various packages. I believe you don't mean it like that, but your question (above) sounds a bit like you want us to solve your problems. That's always a dangerous attitude on lists where help comes only voluntary. You haven't even told us how much memory your computer or OS has. Best wishes, -- Edzer ___ 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] Memory limit problems in R / import of maps
Dylan, Thanks for your note. A student of mine would like to run habitat suitability analysis using the adehabitat package (http://dx.doi.org/10.1890%2F0012-9658%282002%29083%5B2027%3AENFAHT%5D2.0.CO%3B2). I encouraged him to use R, for many reasons. At the moment, he is thinking of doing the whole thing in Matlab (or using the original Biomapper software), because we would not like to give up on the original resolution (250 m). As a GIS person, I definitively do not see ~20 millions pixels as a Huge data set. cheers, Tom Hengl -Original Message- From: Dylan Beaudette [mailto:[EMAIL PROTECTED] Sent: dinsdag 22 april 2008 17:22 To: Tomislav Hengl Cc: r-sig-geo@stat.math.ethz.ch; Michalis Vardakis Subject: Re: [R-sig-Geo] Memory limit problems in R / import of maps On Tue, Apr 22, 2008 at 6:49 AM, Tomislav Hengl [EMAIL PROTECTED] wrote: Dear list, I know that much has already been said about the memory limit problems. If there is any progress about this problem, we would be interested to hear. In our project, we are importing 24 maps/bands, each consists of 1,450,000 pixels. We further would like to glue all maps into a single data frame (e.g. 'ksc' class in adehabitat package; or 'SpatialGridDataFrame' in sp package), but this seems to be impossible. We tried to run this under windows (after following http://cran.r-project.org/bin/windows/base/rw-FAQ.html#There-seems-to-be-a-limit-on-the-memory-it-us es_0021 and setting the --max-mem-size) and under Linux Ubuntu, but still get the same error message (seems that there is no difference in use of memory under the two OS): Error: Cannot allocate vector of size 11.1 Mb The R workspace with 24 loaded grids is also really small (18 MB), but any further gluing and calculation is blocked due the vector size error message. For a comparison, in a GIS such as ArcGIS or SAGA/ILWIS (open source) we have no problems of loading and processing 3-4 times more grids. Should we simply give up on running spatial analysis using large grids (10 million grids) in R? Hi, What exactly were you hoping to do with such a massive data frame once you overcame the initial memory problems associated with loading the data? Any type of multivariate, classification, or inference testing would certainly require just as much memory to perform any analysis on the stack of grids. Not knowing what the purpose of this operation is (although I would guess something related to soil property or landscape modeling of some sort), it is hard to suggest a better approach. For grid that size I would use an algorithm that operates on strips or tiles. There are several great starting points in the GRASS source code. Doing all of the pre-processing, and possibly some aggregating to larger support size, in GRASS would allow you to test any R-centric operations on a coarser version of the original dataset. Cheers, Dylan ___ 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] Memory limit problems in R / import of maps
Dear Miltinho, Please clarify what do you mean by the most important thing is the dimension (number of cols and rows). I the case of rgdal, R reads everything and puts as one variable e.g.: str(maskHG) Formal class 'SpatialGridDataFrame' [package sp] with 6 slots ..@ data :'data.frame': 180 obs. of 1 variable: .. ..$ band1: int [1:18000] 0 0 0 0 0 0 1 1 1 1 ... ..@ grid :Formal class 'GridTopology' [package sp] with 3 slots .. .. ..@ cellcentre.offset: Named num [1:2] 3946500 3247500 .. .. .. ..- attr(*, names)= chr [1:2] x y .. .. ..@ cellsize : num [1:2] 1000 1000 .. .. ..@ cells.dim: int [1:2] 1500 1200 ..@ grid.index : int(0) ..@ coords : num [1:2, 1:2] 3946500 4095500 3247500 3366500 .. ..- attr(*, dimnames)=List of 2 .. .. ..$ : NULL .. .. ..$ : chr [1:2] x y ..@ bbox : num [1:2, 1:2] 3946000 3247000 4096000 3367000 .. ..- attr(*, dimnames)=List of 2 .. .. ..$ : chr [1:2] x y .. .. ..$ : chr [1:2] min max ..@ proj4string:Formal class 'CRS' [package sp] with 1 slots .. .. ..@ projargs: chr +proj=laea +lat_0=52 +lon_0=10 +k=1.0 +x_0=4321000 +y_0=321 +ellps=GRS80 +datum=WGS84 +units=m +towgs84=0,0,0 so I do not think that col/rows makes difference. We typically work with 1200 rows x 1300 columns maps, but would like to use even bigger maps very soon. I really see this problem of not being able to load larger maps to R as the biggest bottleneck (worth investing even time to work on the R code). Tom -Original Message- From: milton ruser [mailto:[EMAIL PROTECTED] Sent: dinsdag 22 april 2008 18:21 To: Tomislav Hengl Cc: Dylan Beaudette; r-sig-geo@stat.math.ethz.ch; Michalis Vardakis Subject: Re: [R-sig-Geo] Memory limit problems in R / import of maps Hi all, In fact, sometimes I fill a little frustrated with some map handling on R. I also had problem with (almost to me) not so large maps reading. The curious is that the simple ArcView 3.2 read the files without problem (in GRD or ASC formats). In fact, I think that the most important thing is the dimension (number of cols and rows) and not the spatial resolution (10 meters, 250 meters etc). So, Tom, how large (cols and rows) are your maps? Kind redards, miltinho On 4/22/08, Tomislav Hengl [EMAIL PROTECTED] wrote: Dylan, Thanks for your note. A student of mine would like to run habitat suitability analysis using the adehabitat package (http://dx.doi.org/10.1890%2F0012-9658%282002%29083%5B2027%3AENFAHT%5D2.0.CO%3B2). I encouraged him to use R, for many reasons. At the moment, he is thinking of doing the whole thing in Matlab (or using the original Biomapper software), because we would not like to give up on the original resolution (250 m). As a GIS person, I definitively do not see ~20 millions pixels as a Huge data set. cheers, Tom Hengl -Original Message- From: Dylan Beaudette [mailto:[EMAIL PROTECTED] Sent: dinsdag 22 april 2008 17:22 To: Tomislav Hengl Cc: r-sig-geo@stat.math.ethz.ch; Michalis Vardakis Subject: Re: [R-sig-Geo] Memory limit problems in R / import of maps On Tue, Apr 22, 2008 at 6:49 AM, Tomislav Hengl [EMAIL PROTECTED] wrote: Dear list, I know that much has already been said about the memory limit problems. If there is any progress about this problem, we would be interested to hear. In our project, we are importing 24 maps/bands, each consists of 1,450,000 pixels. We further would like to glue all maps into a single data frame (e.g. 'ksc' class in adehabitat package; or 'SpatialGridDataFrame' in sp package), but this seems to be impossible. We tried to run this under windows (after following http://cran.r-project.org/bin/windows/base/rw-FAQ.html#There-seems-to-be-a-limit-on-the-memory-it-us es_0021 and setting the --max-mem-size) and under Linux Ubuntu, but still get the same error message (seems that there is no difference in use of memory under the two OS): Error: Cannot allocate vector of size 11.1 Mb The R workspace with 24 loaded grids is also really small (18 MB), but any further gluing and calculation is blocked due the vector size error message. For a comparison, in a GIS such as ArcGIS or SAGA/ILWIS (open source) we have no problems of loading and processing 3-4 times more grids. Should we simply give up on running spatial analysis using large grids (10 million grids) in R? Hi, What exactly were you hoping to do with such a massive data frame once you
[R-sig-Geo] Vacancies for 2 Scientific Programmers, Computational Geo-Ecology, Amsterdam
The Computational Geo-Ecology (CGE) research group of the University of Amsterdam is looking for two enthusiastic Scientific Programmers who will develop virtual laboratories that support data analyses and modelling. This concerns data exploration through interactive visualization, statistical analyses that demand much computer time, and complex dynamic simulation models. In cooperation with other members of the team, they will develop and implement workflows, develop new analysis methodologies and simulation models, and introduce new eScience developments and tools. More info about the CGE research group: http://www.science.uva.nl/ibed-cbpg Salary: The salary will depend on experience, and will range from €2330 to a maximum of €3678 gross per month (salary scale 9 or 10) based on a full-time appointment. For candidates with a doctorate, there is a development opportunity to Ontwikkelaar ICT 2 (salary scale 11, maximum €4284 gross per month). The collective labour agreement of Dutch universities is applicable. Application deadline: 20 May 2008. Vacancies for 2 Scientific Programmers http://www.science.uva.nl/ibed/vacancies.cfm/947EF185-1321-B0BE-A431677965C36F66 ___ R-sig-Geo mailing list R-sig-Geo@stat.math.ethz.ch https://stat.ethz.ch/mailman/listinfo/r-sig-geo
[R-sig-Geo] FYI: 1st call Summer school on spatio-temporal data analysis using R + SAGA + Google Earth, 25-30 August 2008, Amsterdam
*** From 25-30 August 2008 the research school ICG organizes a five-day course targeted at PhD students, academic staff and project teams working with spatio-temporal data. The course aims at providing a balanced combination of theoretical and hands-on-software training. The participants will receive many practical instructions on how to generate R scripts and use SAGA GIS / Google Earth to visualize outputs of analysis. All these packages are available as open-source software or as freeware. After completing the course, the course participants are expected to obtain sufficient knowledge to independently produce R scripts, perform various analyzes and generate outputs that can be directly distributed to the end-users. Flyer: http://www.science.uva.nl/ibed/icg/UvA_geostat2008_1st.pdf Moderators: T. Hengl (uva.nl), D.G. Rossiter (itc.nl), G.B.M. Heuvelink (wur.nl) and V. Olaya Ferrero (unex.es) The course fees is 300 euros and includes the cost of coffee and lunch breaks, materials, a workshop dinner and local transportation. Deadline to register: 15th June 2008 The course is limited to 25 participants. International/overseas applicants are welcome! More info/registration: http://www.science.uva.nl/ibed-agenda/events.cfm/F08AAEA4-1321-B0BE-A4707C0EADE7B462 ___ 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] displayDataset command
If it is of any help, you might want to take a look at this guide: Creating Maps for Publication using R Graphics http://www.nceas.ucsb.edu/scicomp/GISSeminar/UseCases/MapProdWithRGraphics/OneMapProdWithRGraphics.h tml The second example demonstrates how to visualize a remote sensing band using grey levels and then overlay various vector layers. Tom Hengl http://spatial-analyst.net -Original Message- From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] On Behalf Of Roger Bivand Sent: dinsdag 20 mei 2008 16:38 To: Monica Pisica Cc: r-sig-geo@stat.math.ethz.ch Subject: Re: [R-sig-Geo] displayDataset command On Tue, 20 May 2008, Monica Pisica wrote: Hi, I am trying to display a geotiff using displayDataset command from rgdal . which i can do pretty successfully but i cannot display axes for my map i tried to use axes = TRUE and i got the error that axes is matched by multiple arguments and if i try to set xaxs and yaxs as equal to i or r they are simply ignored . Suppose my geotiff is named d1 ... x - GDAL.open(d1) displayDataset(x, zlim = c(-2, 5.7), col = rainbow(20), xaxs = r, yaxs = r) And i get a nice image display but no axes what so ever ;-) If you look in the code, you'll see that axes are set to FALSE internally, and the drawing units are manipulated, so the graphics display does not know where they should go or which units they should be in. Since you can use output.dim= in readGDAL too, I'd go with image(readGDAL(d1), col=grey.colors(20), axes=TRUE) adding arguments to readGDAL() to decimate the image if there are many pixels. You don't get the color palette, though, but it could be retrieved if need be. Roger I will really appreciate any help you can provide, Thanks, Monica _ E-mail for the greater good. Join the i’m Initiative from Microsoft. ood [[alternative HTML version deleted]] -- Roger Bivand Economic Geography Section, Department of Economics, Norwegian School of Economics and Business Administration, Helleveien 30, N-5045 Bergen, Norway. voice: +47 55 95 93 55; fax +47 55 95 95 43 e-mail: [EMAIL PROTECTED] ___ 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] Defining a grid for interpolations ?
FYI I have just finishing preparing an R script that accompanies my paper in Computers in Geosciences on finding the grid cell size for various applications. The script and the data sets can be obtained from: http://spatial-analyst.net/pixel.php This gives a step-by-step guide to estimate suitable grid cell sizes given an input data with its inherent properties: scale, positional accuracy, inspection density, spatial correlation and complexity of terrain. The bottom line is that there are objective ways to calculate a suitable cell size (or at least a range of suitable cell sizes), and that the selection of grid cell size can be automated. For example, using the akima package we can automate estimation of the grid: lbrary(maptools) elevations - readShapePoints(elevations.shp, proj4string=CRS(as.character(NA))) library(akima) elevations.interp - interp([EMAIL PROTECTED],1],[EMAIL PROTECTED],2],z=elevations$VALUE, extrap=T, linear=F) image(elevations.interp, col=bpy.colors(), asp=1) library(spatstat) dem.interp - as.SpatialGridDataFrame.im(as.im(elevations.interp)) [EMAIL PROTECTED] But I would rather first fit a variogram and then select the cell size based on the number of point pairs and given the effective range of spatial autocorrelation. For more info see also: Hengl T., 2006. Finding the right pixel size. Computers and Geosciences, 32(9): 1283-1298. http://dx.doi.org/10.1016/j.cageo.2005.11.008 Tom Hengl http://spatial-analyst.net -Original Message- From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] On Behalf Of Mauricio Zambrano Sent: dinsdag 27 mei 2008 15:27 To: Edzer Pebesma; Paul Hiemstra Cc: r-sig-geo@stat.math.ethz.ch Subject: Re: [R-sig-Geo] Defining a grid for interpolations ? Dear Paul and Edzer, Thanks a lot for your answers. Before reading the solution proposed by Paul, I had tried with: # reading the boundary of the catchment library(maptools) catchment - readShapePoly(only_catchment.shp) catchment.grid - spsample(catchment, cellsize=100, offset = c(0.5, 0.5)) # reading the boundary of the catchment and then I did an IDW interpolation with: pp.idw - krige(PP_DAILY_MEAN_MM~1, meteo, catchment.grid) and it works. However, after reading your solution, I realised that I didn't use the command: gridded(catchment.grid) but I got an interpolation anyway. Where is performed the interpolation when I don't use a gridded grid ( gridded(catchment.grid) ) ? Best regards Mauricio 2008/5/26 Edzer Pebesma [EMAIL PROTECTED]: Please note that this regular grid is still random, as the first point is sampled randomly from the first cell. For a fixed, reproducable grid, add the argument offset = c(0.5, 0.5) -- Edzer Paul Hiemstra wrote: Hi Mauricio, To get a grid based on a shapefile you can use the command spsample. library(sp) library(maptools) source_poly = readShapePoly(/path/to/poly) # cellsize is in map units (e.g. km), also see ?spsample grd = spsample(source_poly, cellsize = c(10e3,10e3), type = regular) gridded(grd) = TRUE # Make it a grid summary(grd) # Visualize the grid spplot(source_poly, sp.layout = list(sp.points, grd)) grd can now be used for interpolations. hth and cheers, Paul Mauricio Zambrano wrote: Dear List, My question is about how to define the grid to be used for the interpolations, using R 2.7.0 and gstat 0.9-47 I'm working in a catchment of ~3000 km2, with daily rainfall data of several stations (more than 40), and I would like to interpolate daily values within all the catchment using ordinary Kriging. For defining the grid in which the interpolation will be carried out, at the beginning, I tried with #setting a grid each 100m vertical and horizontal dx - seq(674400,730700,by=100) dy - seq(4615100,4744400,by=100) catch.grid - expand.grid(dx,dy) #setting the names of the columns of the grid names(catch.grid) - c(x,y) #setting the coordinates of the grid coordinates(catch.grid) - ~x+y #interpolating with the inverse distance pp.idw - krige(PP_DAILY_MEAN_MM~1, meteo_catch_nNA, catch.grid) #plotting the interpolated values spplot(pp.idw[var1.pred], main=Daily Mean PP, [mm]) but looking at the figure I realized that the interpolations were carried out considering all the cells within the squared grid, and not only within the catchment. After, I tried reading a raster file (each 100m) and using it as the grid, DEMM100m - read.asciigrid(Catch_DEM_c100m) but the results that I got seems to be wrong, because I got high values where low values were expected and viceversa. (I really appreciate any help clarifying me what I'm doing wrong ) Is there any way to define a grid, starting from a shapefile of the catchment boundaries ?. For example, I would like to define something similar to the meuse.grid dataset. Thanks in advance and best regards -- hzambran Linux user #454569 -- Ubuntu user #17469 ___
Re: [R-sig-Geo] alter the lon lat lines and coastline according to the grid coordinates
For European projects, we commonly use the official European Terrestrial Reference System (www.euref.eu). You can convert your data from the longlat system to the ETRS using: library(maptools) proj4string(pointmap) - CRS(+proj=longlat +datum=WGS84) pointmap.etrs - spTransform(pointmap, CRS(+init=epsg:3035)) spplot(pointmap.etrs[1]) To simply (without any intervention) generate a gridded surface showing the change of values (the pattern), see: https://stat.ethz.ch/pipermail/r-sig-geo/2008-June/003703.html Many European GI data you can obtain from: EEA (http://dataservice.eea.europa.eu); but also take a look at: 2. Land cover maps (http://www-gem.jrc.it/glc2000/) 3. Landsat images (http://image2000.jrc.it/) 4. General type maps (http://www.inspire-geoportal.eu/) Tom Hengl http://spatial-analyst.net -Original Message- From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] On Behalf Of Costas Douvis Sent: dinsdag 10 juni 2008 17:43 To: r-sig-geo@stat.math.ethz.ch Subject: [R-sig-Geo] alter the lon lat lines and coastline according to the grid coordinates Hi everyone My problem might be easy to deal with but I have already spent quite some time and effort with no result. Any help will be appreciated My data come from an RCM that I ran, Regcm3. My domain spans over the wider area of Europe. The model grid is Lambert conformal. This means that the data values that come from the same line or column of a data matrix do not correspond to the exact same latitude or longitude. I have the values of latitude and longitude of each grid point in 2 matrices with the same dimensions as my data So what I need to do is to plot those data on a map. I guess that the best way to do this is to plot them on a rectangle map, preferably using filled.contour (forgetting for a moment that the grid is actually bended), and afterwards overlay the latitude and longitude lines and the coastline as they really are on that map (i.e. the lat and lon lines should not be straight) Is that possible? Or do you have a better suggestion? -- Kostas Douvis PhD Student University of Athens - Department of Geography and Climatology Academy of Athens - Research Centre for Atmospheric Physics and Climatology email: [EMAIL PROTECTED] tel: +30-210-8832048, +30-210-8847280 fax: +30-210-8842098 ___ 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] kriging
Dear Dave, I separate fitting of the deterministic (trend) and residual part of the universal kriging model all the time. Adding OK of residuals to the trend is fine, as long as the regression model is estimated using GLS (but many do it even if they use only OLS; the difference is often minor). Both KED and RK give the same results, but in order to run KED, you need to have the residuals, so you always have to model the trend-part anyway. See some code at: https://stat.ethz.ch/pipermail/r-sig-geo/2008-April/003433.html In future, please keep the whole history of correspondence. Tom Hengl http://spatial-analyst.net -Original Message- From: Dave Depew [mailto:[EMAIL PROTECTED] Sent: dinsdag 17 juni 2008 14:57 To: [EMAIL PROTECTED] Cc: r-sig-geo@stat.math.ethz.ch Subject: re:kriging Thanks Tom, I've been able to fit a polynomial function to the data quite well. The residuals are behaving (i.e normal distribution and no skewness of variance). I'm assuming this means that I could krige the residuals (Ordinary K?) and then add the trend back to the predicted residual grid? I realize that I won't be able to place confidence limits on the predictions, but the data is primarily to show that we might be able to use GPS hydroacoustic signals to show macrophyte cover and estimate the standing crop. I'm still a newbie when it comes to the theory involved in kriging, but I think I am familiar with the basics...the variogram of the residuals is a nice spherical model (i.e. it reaches a sill about 50% or so above the nugget, and there is little scatter). I am assuming (perhaps wrongly) that the residuals may be modelled with a variogram and then kriged... ___ 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] How to get Applied Spatial Data Analysis with Rbefore 9th of August
I forgot to mention - the book's website is already on-line: http://www.bias-project.org.uk/asdar-book/ The datasets and R scripts are still missing, but I imagine it should all soon be there? (for me personally, R scripts/datasets/examaples are equally important as the book itself) cheers, Tom -Original Message- From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] On Behalf Of [EMAIL PROTECTED] Sent: donderdag 10 juli 2008 14:27 To: r-sig-geo@stat.math.ethz.ch Subject: [R-sig-Geo] How to get Applied Spatial Data Analysis with Rbefore 9th of August Hello, I would like to get my hands on Applied Spatial Data Analysis with R before I leave for summer vacation in the 2nd week of August (it promises to be a great beach reading). It is announced for the 25th of July, so I would like to know if this is still publication date and what is the best vendor (e.g. at Springer or Amazon or maybe the author directly?) to order? I have made the experience that, in particular with new books, that the vendor promises delivery within a few days but in fact needs weeks before actually having it in stock. To complicate things, I live in the Netherlands and not in the US. Any thoughst on that vacation saving question? Best, Stefan ___ 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] How to get Applied Spatial Data Analysis withRbefore 9th of August
Dear Virgilio, Thanks for the info. I still do not mind if you would inform us, via the r-sig-geo mailing list, once the website is 90% finished (the official launch). Again, I feel that for many the website (R scripts, datasets, exercises) will be equally important as the book itself. It would be even better if the website would support some interactivity (e.g. to post comments, see the most popular scripts etc.; see e.g. http://geomorphometry.org/view_scripts.asp?id=6). Looking forward (and congratulations on your book), Tom Hengl http://spatial-analyst.net -Original Message- From: Gomez Rubio, Virgilio [mailto:[EMAIL PROTECTED] Sent: donderdag 10 juli 2008 16:40 To: Tomislav Hengl; [EMAIL PROTECTED] Cc: r-sig-geo@stat.math.ethz.ch Subject: RE: [R-sig-Geo] How to get Applied Spatial Data Analysis withRbefore 9th of August Hi, It is good to see that the book has rised such interest... :) We are not sure about the exact release date yet, but I believe it should be in production now. We are working on the web site at the moment and I am afraid that it is not ready yet. The R scripts and figures are already there and the data sets will come soon. We plan to have everything ready at http://www.asdar-book.org in two or three weeks. What you can see at http://www.bias-project.org.uk/asdar-book is a test that I am doing to produce a web site from the book source files. The final site may look different... hopefully, it will look better. :) But the R scripts should run smoothly (there is no guarantee at the moment though!). Best regards, Virgilio P.S: If you have any other question about the book I would personally prefer to be contacted off list for these issues (self-advertising is not pleasant...). -Original Message- From: [EMAIL PROTECTED] on behalf of Tomislav Hengl Sent: Thu 7/10/2008 2:30 PM To: [EMAIL PROTECTED] Cc: r-sig-geo@stat.math.ethz.ch Subject: Re: [R-sig-Geo] How to get Applied Spatial Data Analysis withRbefore 9th of August I forgot to mention - the book's website is already on-line: http://www.bias-project.org.uk/asdar-book/ The datasets and R scripts are still missing, but I imagine it should all soon be there? (for me personally, R scripts/datasets/examaples are equally important as the book itself) cheers, Tom -Original Message- From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] On Behalf Of [EMAIL PROTECTED] Sent: donderdag 10 juli 2008 14:27 To: r-sig-geo@stat.math.ethz.ch Subject: [R-sig-Geo] How to get Applied Spatial Data Analysis with Rbefore 9th of August Hello, I would like to get my hands on Applied Spatial Data Analysis with R before I leave for summer vacation in the 2nd week of August (it promises to be a great beach reading). It is announced for the 25th of July, so I would like to know if this is still publication date and what is the best vendor (e.g. at Springer or Amazon or maybe the author directly?) to order? I have made the experience that, in particular with new books, that the vendor promises delivery within a few days but in fact needs weeks before actually having it in stock. To complicate things, I live in the Netherlands and not in the US. Any thoughst on that vacation saving question? Best, Stefan ___ 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 mailing list R-sig-Geo@stat.math.ethz.ch https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] generating distance map for Large map
Yes, buffer distances can be smoothly calculated using SAGA (I do not if it will work for your HUGE image, but you can try). Here is an example: . # Import the contour data (from the 1:5k and 1:50k topo maps): library(maptools) getinfo.shape(contours.shp) contours.50 - readShapeLines(contours.shp, proj4string=CRS(as.character(NA))) contours.5 - readShapeLines(contours2.shp, proj4string=CRS(as.character(NA))) # estimate the initial cell size based on the spacing between the contours: dem.area - ([EMAIL PROTECTED],[EMAIL PROTECTED],1])*([EMAIL PROTECTED],[EMAIL PROTECTED],1]) # for beginning, we take a cell size that corresponds to the effective scale: dem.pixelsize - 25 # estimate the buffer distance between the contour lines library(RSAGA) # Download SAGA GIS from www.saga-gis.org and unzip the binaries to: rsaga.env(path=C:/Program Files/saga_vc) # first, convert the contour map to a raster map: rsaga.geoprocessor(lib=grid_gridding, module=3, param=list(GRID=contour50_buff.sgrd, INPUT=contours.shp, FIELD=1, LINE_TYPE=0, USER_CELL_SIZE=dem.pixelsize, [EMAIL PROTECTED],1], [EMAIL PROTECTED],2], [EMAIL PROTECTED],1], [EMAIL PROTECTED],2])) # now extract a buffer distance map (for 'contours') and load it back to R: # the parameters DIST and IVAL are estimated based on the grid properties: rsaga.geoprocessor(lib=grid_tools, module=10, param=list(SOURCE=contour50_buff.sgrd, DISTANCE=contours50_dist.sgrd, ALLOC=contours50_alloc.sgrd, BUFFER=contours50_bdist.sgrd, DIST=sqrt(dem.area)/3, IVAL=dem.pixelsize)) rsaga.sgrd.to.esri(in.sgrds=contours50_dist.sgrd, out.grids=contours50_dist.asc, out.path=getwd(), prec=1) contours50.dist - readGDAL(contours50_dist.asc) . and the complete sample R script: http://spatial-analyst.net/PIXEL/pixelR.zip cheers, Tom Hengl http://spatial-analyst.net -Original Message- From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] On Behalf Of milton ruser Sent: maandag 21 juli 2008 21:11 To: r-sig-geo@stat.math.ethz.ch Subject: [R-sig-Geo] generating distance map for Large map Dear all, I have a very large BINARY image (60,000x50,000 pixels) and need to generate the distance map from the 1 values. When the images are small than my image, I can generate the distance map using ArcGis (my image is in .IMG format), but ArcGis return error for this large one. I am running R 2.7.1 under windows/xp. Is there a way of do this job using R ou RSaga? (unfortunately I don´t have GRASS installed), and my experience with RSAGA is near zero. My computer is 4GB ram. I can load the images on ArcGis and Erdas without problem, but I only don´t know how do this task. Sorry for not send a reproducible code. Any help are welcome. Miltinho Astronauta Brazil [[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] Spatial Resampling
You can also try RSAGA. It is very fast and handy to up-scale large grids (it is perfect to automate processing), and you only need to define the new grid cell size and the up/down-scaling method: library(RSAGA) # Download SAGA GIS from www.saga-gis.org and unzip the binaries to: rsaga.env(path=C:/Program Files/saga_vc) # first, convert to the SAGA format: rsaga.esri.to.sgrd(in.grids=image12m.asc, out.sgrds=image12m.sgrd, in.path=getwd()) # now resample to 18 m grid resolution using the bilinear resampling: rsaga.get.usage(lib=grid_tools, module=0) rsaga.geoprocessor(lib=grid_tools, module=0, param=list(INPUT=image12m.sgrd, GRID=image18m.sgrd, METHOD=1, DIMENSIONS_CELLSIZE=18, SCALE_UP_METHOD=1)) # convert back to ESRI format: rsaga.sgrd.to.esri(in.sgrds=image18m.sgrd, out.grids=image18m.asc, out.path=getwd(), prec=1) The resampling methods implemented in SAGA (see GRID Tools Construction Resampling) are available using the rsaga.get.usage(lib=grid_tools, module=0): -SCALE_UP_METHOD:numInterpolation Method Choice Available Choices: [0] Nearest Neighbor [1] Bilinear Interpolation [2] Inverse Distance Interpolation [3] Bicubic Spline Interpolation [4] B-Spline Interpolation [5] Mean Value -SCALE_DOWN_METHOD:num Interpolation Method Choice Available Choices: [0] Nearest Neighbor [1] Bilinear Interpolation [2] Inverse Distance Interpolation [3] Bicubic Spline Interpolation [4] B-Spline Interpolation Tom Hengl http://spatial-analyst.net -Original Message- From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] On Behalf Of [EMAIL PROTECTED] Sent: vrijdag 25 juli 2008 4:00 To: r-sig-geo@stat.math.ethz.ch Subject: [R-sig-Geo] Spatial Resampling Pujan Raj Regmi [EMAIL PROTECTED] writes: I have a satellite image of 12m spatial resolution (4 bands). I want to resample to coarser level. i.e. to 18m, 24m,30m,36m and so on. The result I got from ENVI are not realistic so I want to do it in R using weighted average from neighboring pixels. In package 'spatstat' the function 'blur' will perform spatial averaging of a pixel image, using Gaussian weights. The function 'as.im' can be used to change the dimensions of a pixel image. So for example if Z is a pixel image (object of class im) then old.pixelsize - Z$xstep blurZ - blur(Z, 1.5 * old.pixelsize) old.dim - Z$dim new.dim - floor(old.dim/2) Y - as.im(blurZ, dimyx=new.dim) computes a Gaussian kernel-weighted average of pixel values, then subsamples the image. Adrian Baddeley ___ 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] [Gstat-info] LOCAL universal kriging with GLOBAL reg. coeff.estimation
Dear Radim, Edzer, I was thinking about the same problem few years ago (I assume that you work with auxiliary maps and not only coordinates). I think (have a feeling) that local and global Universal kriging should be treated as two things (especially if you put a very narrow search radius). This is because, in the case of KED algorithm, both regression and kriging part of predictions are solved at the same time. Hence, if we limit the search window, but keep a constant variogram model, we could obtain very different predictions then if we use a global (regression-kriging) model. Only if the variogram of residuals if absolutely stationary (constant), then we can limit the search window to fit the KED weights. In fact, I am confident that software packages should not allow UK/KED predictions with a limited search radius, because this conflicts the model assumptions (a bit deeper discussion about this problem can be found in my lecture notes section 2.2 Localized versus local models). In the case of regression-kriging, this is less problematic because it is computationally cheap to fit the regression part, and then you are allowed to limit the search radius to interpolate the residuals (we still cheat but not so much). But then the problem is that we do not have an estimate of the RK kriging variance (you could estimate the OLS prediction variance and then OK-residuals prediction variance and then add them together - both can be done in gstat - but then you completely ignore correlation of the two terms; although this is often insignificant). If gstat already does something like this (global GLS simulations of trend with local kriging), then this would be a legitime solution. But I still can not run much geostatistics with big datasets (1000 points, 10^6 grids) in R+gstat (on standard PCs), not to mention geostatistical simulations (that are computationally much more demanding than predictions)... if you would have to do it 100 times, this would really take a lot of processing. One alternative is to run UK/KED in SAGA software (I have tested it - it runs about 5-10 times faster than R+gstat and has no memory limit problems), which is possible through R also, e.g.: rsaga.get.usage(geostatistics_kriging, 3) rsaga.geoprocessor(lib=geostatistics_kriging, module=3, param=list(GRIDS=DEM.sgrd;SLOPE.sgrd;PLANC.sgrd;TWI.sgrd;SINS.sgrd;SMU1.sgrd;SMU3.sgrd;SMU4.sgrd;SM U9.sgrd, GRID=SOLUM_rk.sgrd, SHAPES=baranja.shp, FIELD=0, MODEL=1, NUGGET=0, SILL=200, RANGE=500, INTERPOL=0)) see also http://geomorphometry.org/view_scripts.asp?id=24 But the last version of SAGA had some bug in this module, so I remember that I did not get correct results (I did not test the most recent version of SAGA). I personally think that we should simply think about implementing local RK (regressions/variograms in moving window; maps of regression coefficients, variogram parameters, R-square etc.), and this would then solve this problem (+give us much more insight into the data). Tom Hengl http://spatial-analyst.net -Original Message- From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] On Behalf Of Edzer Pebesma Sent: donderdag 21 augustus 2008 13:07 To: Vasat, Radim Cc: [EMAIL PROTECTED]; sig-geo Subject: Re: [Gstat-info] LOCAL universal kriging with GLOBAL reg. coeff.estimation Good question, I'll include r-sig-geo as well. actually the prediction equations you end up with are kind of funny, and I've never seen them written out. Two different covariance matrices being inversed. Package gstat can do the two-step approach: global BLUE, then simple kriging of residual, but that ignores the correlation of both terms, when added. I'm quite sure that if you do Gaussian simulation however with gstat, the trend is fitted globally (and simulated from the corresponding GLS mean covariance), whereas the kriging is done locally, so that should result in realisations that have the variance you want. You just have to compute the average and variance, to get estimates of the kriging mean and variance you're looking for. Also, it might be computationally demanding, both in time and space (memory). -- Edzer Vasat, Radim wrote: Dear all, My question is how to performe LOCAL universal kriging with estimation of regression coefficients from the whole area (GLOBAL estimation) at once. For large data (I have more than 8000 records) is computationaly too demanding to handle with all the data for UK prediction. Hence, I prefer to do it LOCALLY (namx=100), but the reg. coeff. would be estimated only with limited data in this case (not all data are included). To estimate the reg. coeff. first (with BLUE=TRUE) and then perform the simple kriging with beta equals to the reg. coeff. and nmax=100 might be the solution. But the reg. coeff. estimation error is not included into the final kriging variance in this case. Anybody knows how to join these two requirements (LOCAL UK and GLOBAL
Re: [R-sig-Geo] interp.old and write.ascii.grid
try this: *** library(akima) library(spatstat) library(maptools) ew - as.vector(c(1394180, 1398190, 1399690, 1395780, 1397820, 1395290, 1394995, 1399475, 1410825)) ns - as.vector(c(6200295, 6201020, 6200420, 6205105, 6206590, 6208305, 6210370, 6213715, 6214870)) cures - as.vector(c(48, 52, 59, 84, 50, 46, 114, 92, 48)) sn8 - data.frame(ew, ns, cures) coordinates(sn8) - ~ew+ns sn8.interpold - interp(ew, ns, cures, xo= seq(min(ew), max(ew), length = 400), yo=seq(min(ns), max(ns), length = 300), linear=T) sn8.interpold.grid - as.SpatialGridDataFrame.im(as.im(sn8.interpold)) spplot(sn8.interpold.grid, col.regions=bpy.colors(), sp.layout=list(sp.points, pch=+, sn8)) slot(slot(sn8.interpold.grid, grid), cellsize) - rep(mean(slot(slot(sn8.interpold.grid, grid), cellsize)), 2) write.asciigrid(sn8.interpold.grid[1], asciisn8, na.value = -) *** cheers, Tom Hengl http://spatial-analyst.net -Original Message- From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] On Behalf Of Mikael Carlsson Sent: vrijdag 10 oktober 2008 9:15 To: r-sig-geo@stat.math.ethz.ch Subject: [R-sig-Geo] interp.old and write.ascii.grid Hello! I have a text file that looks like ew;ns;cures 1394180;6200295;48 1398190;6201020;52 1399690;6200420;59 1395780;6205105;84 1397820;6206590;50 1395290;6208305;46 1394995;6210370;114 1399475;6213715;92 1410825;6214870;48 and so on I manage to interpolate with sn8.interpold - interp.old(ew, ns, cures, xo= seq(min(ew), max(ew), length = 400), yo=seq(min(ns), max(ns), length = 300), ncp = 4, extrap=FALSE, duplicate = error, dupfun = NULL) then I want to crate a asciigrid... write.asciigrid(sn8.interpold,asciisn8, na.value = -) and here I got problems... Any tips? Thanks Mikael [[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
Re: [R-sig-Geo] Points density in R
Kernel density can be derived in spatstat package, but also in adehabitat and
many other packages.
You need to loop your operations (output as a list). Note that I use the same
file names to save
space - you only need the output of course!
***
library(maptools)
library(rgdal)
library(spatstat)
txt.list - dir(pattern=glob2rx(*.txt))
txt.list
vegetation.list - as.list(paste(1:length(txt.list),m,sep=))
density.list - as.list(paste(1:length(txt.list),m,sep=))
for(i in seq(along=txt.list)) {
vegetation - read.delim(paste(txt.list[[i]]), sep= )
coordinates(vegetation) - ~X+Y
proj4string(vegetation) - CRS(+init=epsg:26911)
vegetation.ppp - as(vegetation[1], ppp) # here you will probably be
better of if you define
the grid topology yourself (otherwise you get 100x100 pixels by default)
density.list[[i]] - density.ppp(vegetation.ppp, sigma=120)
writeOGR(vegetation[1], paste(vegetation.list[[i]], .shp, sep=),
vegetation, driver=ESRI
Shapefile)
}
image(density.list[[1]])
str(density.list, max.level=2)
***
For future reference, we can not help you if you do not provide any data (there
is a lot of data on
CRAN already). Please also read carefully
http://www.r-project.org/posting-guide.html
yours,
T. Hengl
http://spatial-analyst.net
-Original Message-
From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] On Behalf Of
Alessandro
Sent: Tuesday, October 14, 2008 9:04 PM
To: r-sig-geo@stat.math.ethz.ch
Subject: [R-sig-Geo] Points density in R
Hi all,
I have a txt file with X,Y,Z value. I wish to create a density points Map in
R but sorry I didn't find a right code in R.
Part of my code is this. I have a another question. With many data files
to import, Is It possible don't write 10 lines like my not beautiful code?
Thanks Ale
library(maptools)
library(rgdal)
rsaga.env(path=C:/Progra~1/saga_vc)
library(RSAGA)
library(gstat)
vegetation_2m - read.delim(vegetation_2m.txt, sep= )
vegetation_3m - read.delim(vegetation_3m.txt, sep= )
vegetation_4m - read.delim(vegetation_4m.txt, sep= )
vegetation_5m - read.delim(vegetation_5m.txt, sep= )
vegetation_6m - read.delim(vegetation_6m.txt, sep= )
vegetation_7m - read.delim(vegetation_7m.txt, sep= )
vegetation_8m - read.delim(vegetation_8m.txt, sep= )
vegetation_9m - read.delim(vegetation_9m.txt, sep= )
vegetation_10m - read.delim(vegetation_10m.txt, sep= )
coordinates(vegetation_2m) =~X+Y
coordinates(vegetation_3m) =~X+Y
coordinates(vegetation_4m) =~X+Y
coordinates(vegetation_5m) =~X+Y
coordinates(vegetation_6m) =~X+Y
coordinates(vegetation_7m) =~X+Y
coordinates(vegetation_8m) =~X+Y
coordinates(vegetation_9m) =~X+Y
proj4string(vegetation_2m) - CRS(+init=epsg:26911)
proj4string(vegetation_3m) - CRS(+init=epsg:26911)
proj4string(vegetation_4m) - CRS(+init=epsg:26911)
proj4string(vegetation_5m) - CRS(+init=epsg:26911)
proj4string(vegetation_6m) - CRS(+init=epsg:26911)
proj4string(vegetation_7m) - CRS(+init=epsg:26911)
proj4string(vegetation_8m) - CRS(+init=epsg:26911)
proj4string(vegetation_9m) - CRS(+init=epsg:26911)
#convert in .shp to use RSAGA
writeOGR(vegetation_2m[Z], vegetation_2m.shp, vegetation_2m,
driver=ESRI Shapefile)
writeOGR(vegetation_3m[Z], vegetation_3m.shp, vegetation_3m,
driver=ESRI Shapefile)
writeOGR(vegetation_4m[Z], vegetation_4m.shp, vegetation_4m,
driver=ESRI Shapefile)
writeOGR(vegetation_5m[Z], vegetation_5m.shp, vegetation_5m,
driver=ESRI Shapefile)
writeOGR(vegetation_6m[Z], vegetation_6m.shp, vegetation_6m,
driver=ESRI Shapefile)
writeOGR(vegetation_7m[Z], vegetation_7m.shp, vegetation_7m,
driver=ESRI Shapefile)
writeOGR(vegetation_8m[Z], vegetation_8m.shp, vegetation_8m,
driver=ESRI Shapefile)
writeOGR(vegetation_9m[Z], vegetation_9m.shp, vegetation_9m,
driver=ESRI Shapefile)
[[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] Problems with obtaining a copy of your ASDAR-book
Dear Roger, Edzer and Virgilio, Somehow, I was not lucky enough to receive a copy of your ASDAR book (I only had a chance to browse it yesterday at Edzer's place). First, I followed the R course Edzer gave at GEOENV2008, and our copy did not arrive on time (the organizers made a wrong shipment order; then, they had to send it back after the conference). Second, I just order your book 10 days ago, and it is already out of print. I would like to congratulate you on this excellent book. It helps me systematically adjust many of my data analysis strategies and get a quick review of what and how things should be done (and what should we think about in the near future). The accompanying website is also excellent. The fact that the book has been sold-out so quickly proves my evaluation (I sincerely hope that Springer did not print only a small number of copies?!). An on-line version of each chapter is available at: http://www.springerlink.com/content/978-0-387-78170-9 Unfortunately, it is not a part of standard package that my University buys from Springer, hence I would need to pay 25 USD per chapter (I would rather buy a complete copy for 45 EUR). Hopefully Springer will push this issue and order new prints very soon, so that your book gets in the hands of its biggest fans. (otherwise, you might also want to check some open access publishers in the near future: e.g. http://publications.copernicus.org/). all the best, Tom Hengl http://spatial-analyst.net ___ 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] Problems with obtaining a copy of your ASDAR-book
Dear Virgilio, Thanks for the info. I will try to order it via USA. You are right! The whole book is in fact available on-line for reading via Amazon: http://www.amazon.com/gp/reader/0387781706/ I hope that nobody will be naive and waste 25 USD on Springer. cheers, Tom Hengl http://spatial-analyst.net -Original Message- From: Virgilio Gomez Rubio [mailto:[EMAIL PROTECTED] Sent: Wednesday, October 22, 2008 12:34 PM To: Tomislav Hengl Cc: r-sig-geo@stat.math.ethz.ch Subject: Re: [R-sig-Geo] Problems with obtaining a copy of your ASDAR-book Dear Tom, Many thanks for your nice comments on the book. We realised a few days ago that the book was out of stock in Europe, but it seems that you can still get it from the USA (via Springer or Amazon.com). If you order via through the Springer website (remember the 20% discount) you may need to adjust the zone to North America to see that the book in still in stock. Hopefully they will ship more copies to Europe. As we have mentioned before, any comments on how to improve the book are welcome (off-list, please). Best wishes, Virgilio P.S: By the way, Amazon also allows you to search inside the book and read a few pages at a time. This is useful even if you have the book to look for specific keywords! ___ 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] export .sgrd in tif with R+SAGA
Yes, you can export sgrd to geoTIFFs via the io_grid_gdal module: library(RSAGA) rsaga.get.modules(io_grid_gdal) $io_grid_gdal code name interactive 10Import Raster via GDAL FALSE 21Export Raster via GDAL FALSE 32 Export Raster to GeoTIFF via GDAL FALSE rsaga.get.usage(io_grid_gdal, 2) SAGA CMD 2.0.3 library path: C:/Progra~1/saga_vc/modules library name: io_grid_gdal module name : Export Raster to GeoTIFF via GDAL Usage: 2 -GRIDS str [-FILE str] -GRIDS:str Grid(s) Grid list (input) -FILE:str File File path But it seems that you can not control much about that geoTIFFs, hence I would instead recommend converting the sgrd files to ascii format, and then directly converting them to geoTIFF using the GDAL binaries (DOS prompt; http://www.gdal.org/gdal_translate.html) e.g.: gdal_translate '-of' GTiff 'landsat1.asc' 'landsat1.tif' Unfortunately, SAGA GRID format is not registered at GDAL (maybe you should push this issue? try to write to Olaf Conrad), so that you will always need to first convert to some other GDAL-recognized format before you can do any conversion. Another alternative would be to read the SAGA grid format directly to R, e.g. using the readBin method (https://stat.ethz.ch/pipermail/r-sig-geo/2008-January/003079.html). But I do not really know the structure of the SAGA sgrd format, so I do not know how to set the parameters. Maybe Alex has more info? cheers, Tom Hengl http://spatial-analyst.net -Original Message- From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] On Behalf Of Alessandro Sent: Friday, October 24, 2008 1:59 AM To: r-sig-geo@stat.math.ethz.ch Subject: [R-sig-Geo] export .sgrd in tif with R+SAGA Hi All, Is It possible export a .sgrd grid in .tif image with R+SAGA? Thanks Ale [[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
Re: [R-sig-Geo] subsetting
If the problem is memory handling in R (still a big headache to do much of GIS analysis in R), you should instead consider reading pieces of grids that you need e.g.: info - GDALinfo(dem50m.asc) info gridmaps01 = readGDAL(dem50m.asc, region.dim = round(c(info[[rows]]/2,info[[columns]]))) gridmaps02 = readGDAL(dem50m.asc, region.dim = round(c(info[[rows]]/2,info[[columns]])), offset=round(c(info[[rows]]/2,0))) will split the original map into two tiles. You could practically read any 'tile'/subset of a grid using the 'region.dim' argument (rows and column coordinates). You could first prepare an empty grid topology, then overlay it over your polygons just to get the region dimension, and then read each piece of the grid that you need. In GIS terms, getting a raster from one grid topology to (any) other grid topology (also subgrids) is referred to 'resampling' (I do not think that this is yet implemented in any R package; a combination of akima and spTransform will do the trick but it is not really compact). A generic function to combine R+ILWIS to put a map to a different grid (e.g. to LatLon coordinates) would be: # write to ILWIS: writeGDAL(dem25m[1], dem25m.mpr, ILWIS) # create a new grid: ... gridparameters(geoarc) # resample the map (Bilinear) to the new geographic grid: shell(cmd=paste(ILWIS, crgrf geoarc.grf ,[EMAIL PROTECTED]@cells.dim[[2]], ,[EMAIL PROTECTED]@cells.dim[[1]], -crdsys=LatlonWGS84 -lowleft=(,[EMAIL PROTECTED]@cellcentre.offset[[1]],,,[EMAIL PROTECTED]@cellcentre.offset[[2]],) -pixsize=,[EMAIL PROTECTED]@cellsize[[1]],sep=), wait=F) shell(cmd=paste(ILWIS, dem25m_ll_c.mpr = MapResample(dem25m_c.mpr, geoarc, BiLinear)), wait=F) shell(cmd=paste(ILWIS, open dem25m_ll_c.mpr -noask), wait=F) see http://geomorphometry.org/R.asp hth, Tom Hengl http://spatial-analyst.net -Original Message- From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] On Behalf Of Pieter Beck Sent: Thursday, October 23, 2008 9:15 PM To: r-sig-geo@stat.math.ethz.ch Subject: [R-sig-Geo] subsetting Dear all, I am a new user to the world of spatial applications of R and am wondering about this trivial question: What is the easiest way to subset a SpatialGridDataFrame, based on 4 corner coordinates? The reason I ask is that I want to do an overlay of a polygon file on a large raster dataset, to extract summary statistics from it. Both files are extremely large, however, and just applying overlay causes memory trouble. I was hoping to loop through the polygons, subset the SpatialGridDataFrame using the bounding box surrounding a polygon and then create the summary statistics one polygon at-a-time. To do this, I'd need to subset the SpatialGridDataFrame based on the bounding box corner coordinates. Thanks in advance for the help, Kind Regards, Pieter Beck [[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
Re: [R-sig-Geo] subsetting
to get statistics of grids over polygons, you might also consider using the Grid Statistics for Polygons method in SAGA GIS: rsaga.get.modules(shapes_grid) $shapes_grid code name interactive 1 0Add Grid Values to Points FALSE 2 1 Get Grid Data for Shapes FALSE 3 2 Grid Statistics for Polygons FALSE 4 3Grid Values to Points FALSE 5 4 Grid Values to Points (randomly) FALSE 6 5 Contour Lines from Grid FALSE 7 6 Vectorising Grid Classes FALSE 8 7 Clip Grid with Polygon FALSE 9 8 Gradient from Grid FALSE rsaga.get.usage(shapes_grid, 2) *Note that SAGA can handle even large grids/polygons (up to 2GB size). hth Tom Hengl http://spatial-analyst.net -Original Message- From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] On Behalf Of Pieter Beck Sent: Thursday, October 23, 2008 9:15 PM To: r-sig-geo@stat.math.ethz.ch Subject: [R-sig-Geo] subsetting Dear all, I am a new user to the world of spatial applications of R and am wondering about this trivial question: What is the easiest way to subset a SpatialGridDataFrame, based on 4 corner coordinates? The reason I ask is that I want to do an overlay of a polygon file on a large raster dataset, to extract summary statistics from it. Both files are extremely large, however, and just applying overlay causes memory trouble. I was hoping to loop through the polygons, subset the SpatialGridDataFrame using the bounding box surrounding a polygon and then create the summary statistics one polygon at-a-time. To do this, I'd need to subset the SpatialGridDataFrame based on the bounding box corner coordinates. Thanks in advance for the help, Kind Regards, Pieter Beck [[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
Re: [R-sig-Geo] bad position of png klm file in google earth
Dear Marta,
I am not sure what is the projection system you use for the Algarve case study?
If it is different
from CRS(+proj=longlat +datum=WGS84), then you do need to first reproject
your maps.
You should really look at the examples from the book by Bivand et al.:
http://www.asdar-book.org/code.php?chapter=3figure=1
http://www.asdar-book.org/code.php?chapter=3figure=2
Here you will find small examples of how to export any raster map to GE (it
reprojects a map using
the MapResample method of ILWIS GIS):
http://geomorphometry.org/R.asp
I personally do not use much the maptools to generate GE KMLs, but instead
write directly a KML file
from R e.g. to export a time-series of interpolated temperature images to KML
in a loop:
filename = croclim12.kml
write(?xml version=\1.0\ encoding=\UTF-8\?, filename)
write(kml xmlns=\http://earth.google.com/kml/2.2\;, filename, append =
TRUE)
write(Folder, filename, append = TRUE)
write(\tnameMean Daily TEMP/name, filename, append = TRUE)
write(\topen1/open, filename, append = TRUE)
for (i in slices) {
slice = floor(unclass(as.POSIXct(i))/86400)[[1]]
write(\t\tGroundOverlay, filename, append = TRUE)
write(paste(\t\tname, i, /name, sep = ),
filename, append = TRUE)
write(\t\tTimeSpan, filename, append = TRUE)
write(paste(\t\t\tbegin, slice, /begin, sep = ),
filename, append = TRUE)
write(paste(\t\t\tend, slice + 1, /end, sep = ),
filename, append = TRUE)
write(\t\t/TimeSpan, filename, append = TRUE)
write(\t\tcolor99ff/color, filename, append = TRUE)
write(\t\tIcon, filename, append = TRUE)
write(paste(\t\t\thref, getwd(), /MDTEMP, slice, .png/href, sep =
), filename, append =
TRUE)
write(\t\t\tviewBoundScale0.75/viewBoundScale,
filename, append = TRUE)
write(\t\t/Icon, filename, append = TRUE)
write(\t\taltitude50/altitude, filename, append = TRUE)
write(\t\taltitudeModerelativeToGround/altitudeMode,
filename, append = TRUE)
write(\t\tLatLonBox, filename, append = TRUE)
write(paste(\t\t\tnorth, grids.kml$ylim[[2]],
/north, sep = ), filename, append = TRUE)
write(paste(\t\t\tsouth, grids.kml$ylim[[1]],
/south, sep = ), filename, append = TRUE)
write(paste(\t\t\teast, grids.kml$xlim[[2]],
/east, sep = ), filename, append = TRUE)
write(paste(\t\t\twest, grids.kml$xlim[[1]],
/west, sep = ), filename, append = TRUE)
write(\t\t/LatLonBox, filename, append = TRUE)
write(\t/GroundOverlay, filename, append = TRUE)
}
write(ScreenOverlay, filename, append = TRUE)
write(paste( namevalues from, MDTEMPxlim[[1]], to , MDTEMPxlim[[2]],
/name, sep = ),
filename, append = TRUE)
write( Icon, filename, append = TRUE)
write(paste(\t\t\thref, getwd(), /legend.png/href, sep = ),
filename, append = TRUE)
write( /Icon, filename, append = TRUE)
write( overlayXY x=\0\ y=\1\ xunits=\fraction\
yunits=\fraction\/, filename, append =
TRUE)
write( screenXY x=\0\ y=\1\ xunits=\fraction\ yunits=\fraction\/,
filename, append =
TRUE)
write( rotationXY x=\0\ y=\0\ xunits=\fraction\
yunits=\fraction\/, filename, append =
TRUE)
write( size x=\0\ y=\0\ xunits=\fraction\ yunits=\fraction\/,
filename, append = TRUE)
write( /ScreenOverlay, filename, append = TRUE)
write(/Folder, filename, append = TRUE)
write(/kml, filename, append = TRUE)
I guess you need to read the KML tutorial
(http://code.google.com/apis/kml/documentation/kml_tut.html) to find out what
exactly you need to
use.
hth
Tom Hengl
http://spatial-analyst.net
-Original Message-
From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] On Behalf Of
Marta Rufino
Sent: Friday, October 24, 2008 11:50 AM
To: r-sig-geo@stat.math.ethz.ch
Subject: [R-sig-Geo] bad position of png klm file in google earth
Dear list,
I would like to start by acknowledging the wonderful work of this
support list :-)
I have a question, which should be prety dummy,... but I was looking
around in the net, in reports, etc. and could not sort it out.
I am trying to export a image plot to google earth. I manage to do it
with points data... ok. but when I do it with the images, they simply
go to Fiji (which I am sure it is a nice place... but a bit far from the
Algarve, where it should be).
Even the example given in the help file, the image goes to the same
place. What am I doing wrong?
require(maptools); require(rgdal)
# This is the example, in the file
opt_exask - options(example.ask=FALSE)
qk - SpatialPointsDataFrame(quakes[, c(2:1)], quakes)
proj4string(qk) - CRS(+proj=longlat)
tf - tempfile()
SGqk - GE_SpatialGrid(qk)
png(file=paste(tf, .png, sep=), width=SGqk$width,
height=SGqk$height, bg=transparent)
par(mar=c(0,0,0,0), xaxs=i, yaxs=i)
plot(qk, xlim=SGqk$xlim, ylim=SGqk$ylim, setParUsrBB=TRUE)
dev.off()
kmlOverlay(SGqk, paste(tf, .kml, sep=), paste(tf, .png, sep=))
# than the files produced are in
tf()
I go to the file, pick it the two of them and drag it in google earth...
and there they are, in Fiji!!! I am sure it is silly, but it is my first
steps in googleearth--R
Re: [R-sig-Geo] problem with merge grids in R+SAGA
Alessandro, You need to separate the put the list of input grids as a character vector e.g.: GRIDS=c(DSM1.sgrd,DSM2.sgrd,DSM3.sgrd) See also: http://spatial-analyst.net/wiki/index.php?title=Software#SAGA_GIS Maybe you should also study this posting guide: http://spatial-analyst.net/wiki/index.php?title=Scripting_in_R -Original Message- From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] On Behalf Of [EMAIL PROTECTED] Sent: Thursday, November 13, 2008 9:20 AM To: r-sig-geo@stat.math.ethz.ch Subject: [R-sig-Geo] problem with merge grids in R+SAGA Importance: High Hi All when I use this function with R and SAGA: rsaga.get.usage(grid_tools, 3) rsaga.geoprocessor(lib=grid_tools, module=3, param=list(GRIDS=DSM1.sgrd,DSM2.sgrd,DSM3.sgrd,GRID_TARGET=DSM.sgrd,TYPE=7, INTERPOL=2, OVERLAP=0)) I have this problem SAGA CMD 2.0.3 library path: C:\Windows/modules library name: grid_tools module name : Merging author : (c) 2003 by O.Conrad Load grid: DSM1.sgrd,DSM2.sgrd,DSM3.sgrd... failed error: Grid file could not be opened. error: input file [DSM1.sgrd,DSM2.sgrd,DSM3.sgrd] error: empty input list [GRIDS] error: executing module [Merging] I don't understand very well where is the problem of the code. thanks for help Ale ___ 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] Bug report in Merged module di R+SAGA
Hi Alessandro, Please do submit such and similar bugs to Olaf Conrad (the main developer of SAGA). For your information, the C++ code of SAGA (SAGA 2.0.2 API - Python Interface) is available from http://sourceforge.net/projects/saga-gis/ which means that you could obtain the code and fix the bug yourself. You could also try using the ILWIS MapGlue command e.g. ILWIS - C:\\Progra~1\\N52\\SetupIlwis\\IlwisClient.exe -C shell(cmd=paste(ILWIS, GRID_AB.mpr = MapGlue(GRID_A, GRID_B, replace)), wait=F) shell(cmd=paste(ILWIS, open GRID_AB.mpr -noask), wait=F) This is very flexible as it allows you to take unlimited (?) number of maps and resample them to a new common grid (something like the mosaicing tool in Erdas Imagine). See also: http://spatial-analyst.net/ILWIS/htm/ilwisapp/glue_raster_maps_functionality.htm (I guess that merging of maps is also possible in the GRASS GIS?) In any case, welcome to world of the open source software! Tom Hengl http://spatial-analyst.net/wiki/index.php?title=Software -Original Message- From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] On Behalf Of Alessandro Sent: Friday, November 14, 2008 8:01 PM To: r-sig-geo@stat.math.ethz.ch Subject: [R-sig-Geo] Bug report in Merged module di R+SAGA Hi All, There is a bug in the SAGA merged module in R. I tried it with SAGA 2.0.3 from within R, with different grids, and It's ran forever. This must be a bug in the Merged module, specifically in the command line version of SAGA (saga_cmd). It's been a problem in some other problems too that the saga_cmd binding crashed but the saga_gui version of the same modules worked. Today I report this bug in Sourcefourge project web site for SAGA. Ale [[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
Re: [R-sig-Geo] ramping four color in xy space
Dear Holger, This is a small script that uses [colorspace] library to run calculations with colors (some methods are limited to 256 colors only): http://spatial-analyst.net/scripts/whitening.R If you assign R,G,B values spatially and then interpolate them (they have to be in the 0-1 range), then it is relatively easy to produce an RGB image. cheers, Tom Hengl http://spatial-analyst.net/wiki/index.php?title=Uncertainty_visualization PS: maybe you should also consider color mixing: Hengl T., Walvoort D.J.J., Brown A. 2004. A double continuous approach to visualisation and analysis of categorical maps. Int. Jou. of Geographical Information Science, 18(2): 183-202. http://dx.doi.org/10.1080/13658810310001620924 -Original Message- From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] On Behalf Of Holger Kreft Sent: Tuesday, November 18, 2008 9:38 AM To: R-sig-Geo@stat.math.ethz.ch Subject: [R-sig-Geo] ramping four color in xy space Hi, I have run a metaMDS(vegan) on ecological community data and have a nice plot in xy space. What I want to do next is to assign a unique color to each point in the plot and then look at how these communities map out in geographic space. So, each corner of the plot has a color (let's say green, blue, red and yellow). I guess I have to ramp the colors for the space inbetween. The right plot here illustrate how the palette should look like: http://www.flickr.com/photos/erniebooth/183647844/ Does anyone have experience in ramping four colours in 2-D? I looked at a lot of packages, but they are all not doing what I want. Any hint / help is greatly appreciated! Thanks, Holger ___ 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] Geonames elevation help
You could also loop e.g. (this way you have a bit more control):
long - c(-81.66,-82.66)
lat - c(36.21,37.21)
df - data.frame(lat,long)
df$srtm - rep(NA, length(df$lat))
for(i in 1:length(df$srtm)){
+ df$srtm[i] - GNsrtm3(df$lat[i],df$long[i])$srtm3
+ }
df
lat long srtm
1 36.21 -81.66 990
2 37.21 -82.66 500
But Barry's solution is more straight forward.
Tom Hengl
-Original Message-
From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] On Behalf Of
Barry Rowlingson
Sent: Tuesday, November 25, 2008 5:40 PM
To: r-sig-geo@stat.math.ethz.ch
Subject: Re: [R-sig-Geo] Geonames elevation help
2008/11/25 Michael Denslow [EMAIL PROTECTED]:
I am hoping that someone can help me with the GNsrtm3 function in the
Geonames package. I am
trying to get elevation values for more than one lat,long stored in a data
frame. I am not sure if
this function can process more than one request since it is a web query. But
since R seems to be
able to do anything, I figured I was missing something.
I am an R novice so any help would be most appreciated. Here is a shortened
example of what I am
trying to do.
library(geonames)
long - c(-81.66,-82.66)
lat - c(36.21,37.21)
df - data.frame(lat,long)
GNsrtm3(df$lat,df$long)
# the result below is only for the first set of coordinates
srtm3lng lat
1 990 -81.66 36.21
I also tried using apply, but it also only seemed to work for the first set
of coordinates.
Thanks in advance,
Michael
Maybe you were using apply wrong.
With your df, you can do:
df$srtm3=apply(df,1,function(l){GNsrtm3(l[1],l[2])$srtm3})
df
lat long srtm3
1 36.21 -81.66 990
2 37.21 -82.66 500
apply() passes the row of the dataframe as a vector to the function,
so I just get the first and last element and call GNsrtm3, and strip
off just the result.
Oh, make sure I've got the lat-long the right way round.
Barry
___
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] Regression Kriging
Dear Pritam, Please forward us the R code, and, if possible, upload your data to some ftp server. The coordinate system string of the maps imported to R can be set using: data(meuse) coordinates(meuse) - ~x+y proj4string(meuse) - CRS(+init=epsg:28992) A list of epsg code is available via http://epsg-registry.org FYI: ILWIS GIS can be also controlled from R (see http://spatial-analyst.net/wiki/index.php?title=Software#ILWIS_GIS), regression-kriging/universal kriging you can run from a variety of software (http://spatial-analyst.net/wiki/index.php?title=Regression-kriging_guide). Personally, I would go for a combination of R+SAGA (e.g. http://geomorphometry.org/view_scripts.asp?id=24) to generate maps of large areas. First try all these combinations and then get back to us once you have a concrete problem with the code (note that this is an R mailing list, and not an ILWIS mailing list). I am sure we will find a solution to your problem. Tom Hengl http://spatial-analyst.net -Original Message- From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] On Behalf Of Pritam Sharma Sent: Thursday, November 27, 2008 9:51 AM To: r-sig-geo@stat.math.ethz.ch Subject: [R-sig-Geo] Regression Kriging Hello Sir/Madam This is Pritam Chand, working as Technical Associate in Forest Survey of India, Dehrdaun. I have some doubts regarding regression kriging. Actually I am trying to make final regression kriging predicted map and using R and ILWIS software. I am able to reach steps upto variogram analysis of residuals but in last command for regression kriging it is giving error that the predictor variables are in different coordinate systems it might be because of using ILWIS.. as it is its basic problem I am using 8 predictor variable maps..and with the help of these predicting soil organic carbon..All the steps are running except the final one command of regression kriging.. Please help me regarding this query or kindly send me some link or email from where I get help regrading this query. If possible then plz send me sample code for this problem so I can get idea and go forward with my objectives. With regards, Pritam Chand Sharma, FSI, Dheradun, India -- The Earth was not given to us by our parents. It was loaned to us by our children. [[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] FW: 1st Call for Papers: Geomorphometry 2009, 29 August - 2nd September, Zurich, Switzerland
Apologies for cross postings! FIRST CALL FOR PAPERS Geomorphometry 2009 29 August - 2 September 2009 Zurich, Switzerland http://2009.GEOMORPHOMETRY.ORG e-mail: [EMAIL PROTECTED] PROGRAM CHAIRS Ross Purves University of Zurich Stephan Gruber University of Zurich Tomislav Hengl University of Amsterdam KEY DATES Workshop proposals due 14 January 2009 Extended abstracts due 1 March 2009 Notification of acceptance 1 April 2009 Final camera-ready digital manuscripts due 1 May 2009 Author registration deadline 15 May 2009 Early registration deadline 15 May 2009 Geomorphometry 2009 Workshops 29 August 30 August 2009 Geomorphometry 2009 31 August - 2 September 2009 AIMS AND SCOPE The aim of Geomorphometry 2009 is to bring together researchers to present and discuss developments in the field of quantitative modelling and analysis of elevation data. Geomorphometry is the science of quantitative land-surface analysis and description at diverse spatial scales. It draws upon mathematical, statistical and image-processing techniques and interfaces with many disciplines including hydrology, geology, computational geometry, geomorphology, remote sensing, geographic information science and geography. The conference aims to attract leading researchers in geomorphometry presenting methodological advances in the field and to provide young researchers with an opportunity to present new results. The Geomorphometry 2009 conference will continue a series initiated by the Terrain Analysis and Digital Terrain Modelling conference hosted by Nanjing Normal University in November 2006. Possible topics include, but are not limited to: - Extraction of land-surface parameters from DEMs - Implications of novel data sources - Identification and classification of land-surface objects - Uncertainty in geomorphometry - Semantics of land-surface description - Visualisation in geomorphometry - Implications of scale and resolution - Flow and hydrological modelling using DEMs - Efficient methods for application to large data sets - Novel applications of geomorphometry - Planetary geomorphometry We specifically aim at papers with new methodological insights and thus papers which simply describe the application of GIS are discouraged. CONFERENCE PROGRAMME The conference programme will be based around a single track of papers, all of which will be subject to review in the form of extended abstracts by members of the programme committee. Criteria for paper acceptance will include relevance to the conference, novelty, scientific significance, relation to previous work in the domain and the quality of presentation. The proceedings will be made available both digitally and as printed working materials to attendees at the time of the conference and archived online. A special issue of a journal, to which authors will be invited to submit full papers after the conference is also planned. WORKSHOPS Geomorphometry will host up to three workshops, each with 15-30 attendees on the 29th and 30th August. We invite applications to host a workshop on a theme related to the main conference. Workshops should primarily take the form of either tutorials in a particular method or technique, or provide the opportunity for detailed discussion of upcoming topics. They should not simply be mini-conferences. If you are interested in organising a workshop, please mail [EMAIL PROTECTED] with WORKSHOP as the subject line, as well as a 1-page description of the workshop with the following headings: intended aims and scope, intended audience, outline workshop programme and technical requirements. SUBMISSIONS Prospective authors will be invited to submit extended abstracts of up to 2000 words by the above deadline through the EasyChair system. Formatting instructions and detailed information on using the submission system will be available in due course. Extended abstracts must be original works by the authors, not be currently under review in the same form by another outlet and not submitted elsewhere prior to the notification date. SCIENTIFIC COMMITTEE Alexander Brenning University of Waterloo, Canada Ian Evans Durham University, UK Peter Fisher University of Leicester, UK John Gallant CSIRO, Australia Paul Gessler University of Idaho, USA Stephan Gruber University of Zurich, Switzerland Tomislav Hengl University of Amsterdam, Netherlands Oliver Korup WSL, Switzerland Helena Mitasova North Carolina State University, USA Scott Peckham Rivix, USA Hannes Reuter Joint Research Centre, Italy Robert Weibel University of Zurich, Switzerland John Wilson University of Southern California, USA Jo Wood
Re: [R-sig-Geo] resampling grids
A quick way to rescale a set of maps is to use the grid tools of SAGA GIS (see
rsaga.get.usage(grid_tools, 0)). First, download all ASCII file maps to a
working directory, then
obtain the list of maps in the folder, and then automate:
ascmaps - list.files(getwd(), pattern=\\.asc$, full=F)
levels(ascmaps) - list.files(getwd(), pattern=\\.asc$, full=F)
ascmaps500 - ascmaps
levels(ascmaps500) - paste(500m,levels(ascmaps),sep=)
library(RSAGA)
rsaga.esri.to.sgrd(in.grids=levels(ascmaps),
out.sgrds=set.file.extension(levels(ascmaps),.sgrd), in.path=getwd())
rsaga.get.usage(grid_tools, 0)
for (i in 1:length(ascmaps)) {
rsaga.geoprocessor(lib=grid_tools, module=0,
param=list(INPUT=set.file.extension(levels(ascmaps)[i],.sgrd),
GRID=set.file.extension(levels(ascmaps500)[i],.sgrd), METHOD=0,
DIMENSIONS_CELLSIZE=500,
SCALE_UP_METHOD=1))
}
rsaga.sgrd.to.esri(in.sgrds=set.file.extension(levels(ascmaps500),.sgrd),
out.grids=set.file.extension(levels(ascmaps500),.asc), out.path=getwd(),
prec=1)
unlink(*.hgrd) # delete the temporary files (SAGA)
unlink(*.sgrd)
unlink(*.sdat)
Note that you only need to set the resolution of the output maps and/or the
rescaling method:
-SCALE_UP_METHOD:numInterpolation Method
Choice
Available Choices:
[0] Nearest Neighbor
[1] Bilinear Interpolation
[2] Inverse Distance Interpolation
[3] Bicubic Spline Interpolation
[4] B-Spline Interpolation
[5] Mean Value
-SCALE_DOWN_METHOD:num Interpolation Method
Choice
Available Choices:
[0] Nearest Neighbor
[1] Bilinear Interpolation
[2] Inverse Distance Interpolation
[3] Bicubic Spline Interpolation
[4] B-Spline Interpolation
Mechanical down-scaling from 20 to 10 m is certainly doable but you will be
'cheating' right?
More stuff on this topic (downscaling) in e.g.:
Marc F. P. Bierkens, Peter A. Finke, P. de Willigen, 2000. Upscaling and
Downscaling Methods for
Environmental Research. Springer, ISBN 0792363396, 9780792363392, 190 pages.
http://www.amazon.com/gp/reader/0792363396/
Tom Hengl
http://spatial-analyst.net/wiki/index.php?title=Software
-Original Message-
From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] On Behalf Of
Jason Gasper
Sent: Wednesday, December 10, 2008 7:23 PM
To: Dave Depew
Cc: r-sig-geo@stat.math.ethz.ch
Subject: Re: [R-sig-Geo] resampling grids
Hello David,
I am not sure I completly understand what you want to do, but Hawths
tools (a free add in for arcMAP) has a resampling function that may be
of use: http://www.spatialecology.com/htools/rndselss.php. Spatial
analyst also has some resampling functions. If you have access to ESRI
you could read the data into R using the maptools. package.
Dave Depew wrote:
Hi all,
I apologize in advance, I haven't the time today to comb the help list
for tips and hints, but does anyone know how to resample a 20 by 20m
grid so the operational spacing is down to a 10 by 10m grid?
Any assistance is MUCH appreciated!
___
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] Converting kml to shape (readOGR)
Dear Henk,
My experiences is that the easiest way to read KML files and convert them to
sp-type of objects is
to use the XML package.
Here is a small example:
library(XML)
meuse_lead.kml - xmlTreeParse(meuse_lead.kml)
lengthp - length(meuse_lead.kml$doc[[1]][[1]][[1]])-1
lead_sp - data.frame(Longitude=rep(0,lengthp), Latitude=rep(0,lengthp),
Var=rep(0,lengthp))
for(j in 1:lengthp) {
LatLon -
unclass(meuse_lead.kml$doc[[1]][[1]][[1]][j+1][[1]][2][[1]][[1]][[1]])$value
Var -
unclass(meuse_lead.kml$doc[[1]][[1]][[1]][j+1][[1]][1][[1]][[1]])$value
lead_sp$Longitude[[j]] - as.numeric(matrix(unlist(strsplit(LatLon,
split=,)), ncol=2)[1])
lead_sp$Latitude[[j]] - as.numeric(matrix(unlist(strsplit(LatLon,
split=,)), ncol=2)[2])
lead_sp$Var[[j]] - as.numeric(matrix(unlist(strsplit(strsplit(Var,
split=i)[[1]][2],
split=/i)), ncol=2)[1])
}
coordinates(lead_sp) -~Longitude+Latitude
proj4string(lead_sp) - CRS(+proj=longlat +ellps=WGS84)
bubble(lead_sp, Var)
every KML can have a completely different structure, so you will need to modify
position of
coordinates / variables where needed.
see also:
http://spatial-analyst.net/wiki/index.php?title=Export_maps_to_GE
HTH
-Original Message-
From: r-sig-geo-boun...@stat.math.ethz.ch
[mailto:r-sig-geo-boun...@stat.math.ethz.ch] On Behalf Of
Edzer Pebesma
Sent: Tuesday, January 06, 2009 12:16 PM
To: Henk Sierdsema
Cc: r-sig-geo@stat.math.ethz.ch
Subject: Re: [R-sig-Geo] Converting kml to shape (readOGR)
Henk Sierdsema wrote:
Hi,
I've been trying to convert a vector kml-file to an ESRI shape with
readOGR(), but I don't manage
to read the kml-file. Can anyone help me with this?
http://www.gdal.org/ogr/ogr_formats.html mentions at KML that read
support needs libexpat. Further it says KML reading is only available
if GDAL/OGR is built with the Expat XML Parser, otherwise only KML
writing will be supported.
Maybe this was missing while your gdal/rgdal was built?
--
Edzer
Thnx,
Henk Sierdsema
SOVON Vogelonderzoek Nederland / SOVON Dutch Centre for Field Ornithology
Rijksstraatweg 178
6573 DG Beek-Ubbergen
The Netherlands
tel: +31 (0)24 6848145
fax: +31 (0)24 6848122
___
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.springer.com/978-0-387-78170-9 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 mailing list
R-sig-Geo@stat.math.ethz.ch
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] combining two image files into one database
Dear Alok, Maybe you should consider converting the im data to sp class. Then, you have as many grids in a single SpatialGridDataFrame e.g.: library(maptools) library(rgdal) grids - as(bei.extra[[1]], SpatialGridDataFrame) names(grids)[1] - elev grids$grad - as(bei.extra[[2]], SpatialGridDataFrame)$v spplot(grids[elev]) ... grids$achan - readGDAL(achan.asc)$band1 save(grids, grids.image) see also: http://spatial-analyst.net/wiki/index.php?title=Species_Distribution_Modelling HTH, T. Hengl -Original Message- From: r-sig-geo-boun...@stat.math.ethz.ch [mailto:r-sig-geo-boun...@stat.math.ethz.ch] On Behalf Of Alok K. Bohara, PhD Sent: Wednesday, January 07, 2009 4:21 PM To: Greg Lee Cc: r-sig-geo@stat.math.ethz.ch Subject: Re: [R-sig-Geo] combining two image files into one database Thanks Greg for your response: I want to do what the examples in Spatstat do with regard to the point and pixel/image data bei and bei.extra.. (bei.extra contains a list of two image files --elevation and gradient. They plot these two items (and also take log, I believe) In my case: X_image Y_image (I converted both from text matrix to pixel/image using im(...) command) I want to combine (or list) X_image and Y.image as a part of a single file -- XY_image, and save and do the following load(XY.image) plot(XY_image)# plot both images side by side plot(XY_image$X_image) # plot just the X_image etc... Thank you for your help. # I tried doing the list command: XY_image -list(X_image,Y_image) but it returns NULL when I do XY$X_image Sincerely, Alok Bohara Greg Lee wrote: Hello Alok, We really need more details about what you are trying to acheive in order to help you. c(X, Y) concatenates the objects X and Y. The result will be a vector. If X and Y happen to be m x n matrices, the result will be a vector of length 2 * m * n: probably not what you want. You can see this for yourself. Eg: X - matrix(1, 3, 3) Y - matrix(2, 3, 3) c(X, Y) If you were hoping to combine the information in corresponding pixels of the matrices X and Y, you will need to have a clear idea of *how* you would like the information to be combined. ie: Are you looking for the mean value of the pixels, or is some other operation required? Eg 2: X - matrix(1:9, 3, 3) Y - matrix(2, 3, 3) X + Y / 2 Hope that helps. Greg 2009/1/7 Alok K. Bohara, PhD boh...@unm.edu mailto:boh...@unm.edu Hi: Let me start by saying that I am new to R. I converted two pixel ascii matrices into image files: X_im - im(X) # X is an ascii pixel matrix Y_im - im(Y) # Y is an ascii pixel matrix How can I combine these two images into 1) one datafile XY_im, 2) save it and 3) read it after loading (e.g., x - XY_im$X_im. and plot it.) I tried doing this: XY_im - c(X_im,Y_im) save(XY _im, file = XM_im.Rdata) load(XY_im.Rdata) x- XY_im$X_im# It does not like $ and says it is atomic something... plot(x) # of course it does not work either. What am I doing wrong? Any tips? Thanks you. Alok Bohara -- Alok K. Bohara, Ph.D. Professor Department of Economics University of New Mexico Albuquerque, NM 87131,USA Ph: 505-277-5903/5304(w) Fax:505-277-9445 email: boh...@unm.edu mailto:boh...@unm.edu http://www.unm.edu/~econ/faculty/professors.html http://www.unm.edu/%7Eecon/faculty/professors.html Nepal Study Center: http://nepalstudycenter.unm.edu http://nepalstudycenter.unm.edu/ ___ R-sig-Geo mailing list R-sig-Geo@stat.math.ethz.ch mailto:R-sig-Geo@stat.math.ethz.ch https://stat.ethz.ch/mailman/listinfo/r-sig-geo -- Alok K. Bohara, Ph.D. Professor Department of Economics University of New Mexico Albuquerque, NM 87131,USA Ph: 505-277-5903/5304(w) Fax:505-277-9445 email: boh...@unm.edu http://www.unm.edu/~econ/faculty/professors.html Nepal Study Center: http://nepalstudycenter.unm.edu ___ 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] Spatial data analysis in MATLAB / Comparison of MATLAB vs R
Dear Barry, Thanks for your reply. Funnily enough, I got the same reply from the module coordinator - you will enjoy learning MATLAB. I have accepted to do everything in MATLAB (I do not have much choice really). Then, I might also try Scilab afterwards and controlling Scilab from R. I will keep you informed about what I discovered about MATLAB (at least considering the spatial analysis capabilities). The real problem is that students want to get training in ESRI and Mathworks products, because they think that it gives them better chances to find work in large/government companies (and this is still largely true). So I think that we really need to 'work' on government agencies. At least when I look at the Google trend graphs, I can see that some of the commercial players are irreversibly flowing into a decay function :) Tom Hengl http://spatial-analyst.net -Original Message- From: b.rowling...@googlemail.com [mailto:b.rowling...@googlemail.com] On Behalf Of Barry Rowlingson Sent: Tuesday, January 20, 2009 10:17 AM To: Tomislav Hengl Cc: r-sig-geo@stat.math.ethz.ch Subject: Re: [R-sig-Geo] Spatial data analysis in MATLAB / Comparison of MATLAB vs R 2009/1/20 Tomislav Hengl t.he...@uva.nl: I am preparing to teach spatial data analysis at BSc level (environmental and Earth sciences). The study programme is completely based on MATLAB, which means that I will also need to adjust (I do not have much experience with MATLAB). If you convert your entire study programme to R you'll annoy your colleagues who have to learn and rewrite their courses in R. If you just convert your spatial data analysis module to R you'll annoy your students who'll have to learn R for this and MATLAB for everything else. You choose :) Ideally all your colleagues will take to R and feel joy at the prospect of learning new software and rewriting course notes - but how often does that happen? As a first step you could try introducing Scilab or Octave to the students so they can use a MATLAB-ish package at zero cost. For our undergrad maths course they use Scilab for linear algebra-type stuff, R for stats, and Maxima for computer algebra. You may need to really rework your programme if you want to use just R throughout. It could be worth it. Don't forget the big argument. R is free and open source - it can, like the research we publish, be freely shared. For me that's ultimately compelling - if it's not open, it's not science; and if it's not science, it's not a BSc. Barry ___ R-sig-Geo mailing list R-sig-Geo@stat.math.ethz.ch https://stat.ethz.ch/mailman/listinfo/r-sig-geo
[R-sig-Geo] The geonames.org server changed / geonames package
Hi Barry, Edzer, Just to let you know that geonames.org server changed from ws.geonames.org to ws5.geonames.org See the info at: http://www.geonames.org/export/ Maybe you should add to your R package geonames an option to specify the server name manually, e.g.: setGNserver - ws5.geonames.org PS: Does anybody know more about how to utilize the Google maps services in R (http://code.google.com/apis/maps/documentation/services.html)? (i.e. to get the same type of output as in the case of geonames.org). These need to be combined with some Java script or similar right? My idea is to automate the way to save an image shown in Google maps (e.g. as TIF), then obtain the bounding coordinates, import to a GIS and use the image (high res satellite imagery) as a GIS layer. Tom Hengl http://spatial-analyst.net ___ R-sig-Geo mailing list R-sig-Geo@stat.math.ethz.ch https://stat.ethz.ch/mailman/listinfo/r-sig-geo
[R-sig-Geo] A small repository of publicly available Global maps
Hi all, I have just finished preparing a small repository of publicly available Global maps on environment. For more info see: http://spatial-analyst.net/wiki/index.php?title=Global_datasets The maps can be downloaded directly from: http://spatial-analyst.net/worldmaps/ I would be interested to find out if any of you are aware of a similar dataset that is not mentioned in this article? (of course, I am only interested in the publicly available datasets) If somebody knows a way how to make these maps even more accessible to the R-spatial community, please let me know (e.g. to include them in the maps package). Tom Hengl http://home.medewerker.uva.nl/t.hengl/ ___ 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] calculate raster values based on vector regions
Dear Herry, If I understand what you problem, one solution is to use R+SAGA. You should first convert the polygon map to the same grid, and then you can load it to R and do any type of aggregation: library(maptools) library(rgdal) library(RSAGA) # load the gridded map: rastermap - readGDAL(rastermap.asc) # load the shapefile: rsaga.esri.to.sgrd(in.grids= rastermap.asc, out.sgrds= rastermap.sgrd, in.path=getwd()) # convert the polygon map to a raster map: cellsize - raster...@grid@cellsize[1] rsaga.geoprocessor(lib=grid_gridding, module=3, param=list(GRID=polygons.sgrd, INPUT=polygons.shp, FIELD=1, LINE_TYPE=0, USER_CELL_SIZE=cellsize, user_x_extent_min=raster...@bbox[1,1]+cellsize, user_x_extent_max=raster...@bbox[1,2]-cellsize, user_y_extent_min=raster...@bbox[2,1]+cellsize, user_y_extent_max=raster...@bbox[2,2]-cellsize)) rsaga.sgrd.to.esri(in.sgrds=polygons.sgrd, out.grids=polygons.asc, out.path=getwd(), prec=0) rastermap$polygons - as.factor(readGDAL(polygons.asc)) # summary statistics per polygon class: raster.polygons - boxplot(band1 ~ polygons, rastermap @data, col=rainbow(length(levels(negotingrid$SOIL str(raster.polygons) You will of course need to adjust the FIELD position in the attribute table and the precision prec. To obtain and use SAGA, check this article: http://spatial-analyst.net/wiki/index.php?title=Software cheers, How is the weather now in Canberra? In Amsterdam is freezing brrr. T. Hengl -Original Message- From: r-sig-geo-boun...@stat.math.ethz.ch [mailto:r-sig-geo-boun...@stat.math.ethz.ch] On Behalf Of alexander.h...@csiro.au Sent: Wednesday, January 28, 2009 2:10 AM To: r-sig-geo@stat.math.ethz.ch Subject: [R-sig-Geo] calculate raster values based on vector regions Hi List, How would I go about assigning values to a raster based on polygon regions in R. Examples would be most appreciated. I have (vector) regions assigned to a specific value. The raster has NAs and some pixels where these values are likely to occur on ground. I need to assign these values to the raster and can do this in ArcInfo through vectors and rasterizing but only to a limited raster size. And R is much more preferable anyway... Any help appreciated. Cheers Herry ___ 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] calculate raster values based on vector regions
My appologies. The previous script had few (5) typos: library(rgdal) library(RSAGA) # load the gridded map: rastermap - readGDAL(rastermap.asc) rsaga.esri.to.sgrd(in.grids= rastermap.asc, out.sgrds= rastermap.sgrd, in.path=getwd()) # convert the polygon map to a raster map: cellsize - raster...@grid@cellsize[1] rsaga.get.usage(lib=grid_gridding, module=3) rsaga.geoprocessor(lib=grid_gridding, module=3, param=list(GRID=polygons.sgrd, INPUT=polygons.shp, FIELD=1, LINE_TYPE=0, USER_CELL_SIZE=cellsize, user_x_extent_min=raster...@bbox[1,1]+cellsize, user_x_extent_max=raster...@bbox[1,2]-cellsize, user_y_extent_min=raster...@bbox[2,1]+cellsize, user_y_extent_max=raster...@bbox[2,2]-cellsize)) rsaga.sgrd.to.esri(in.sgrds=polygons.sgrd, out.grids=polygons.asc, out.path=getwd(), prec=0) rastermap$polygons - as.factor(readGDAL(polygons.asc)$band1) # summary statistics per polygon class: raster.polygons - boxplot(band1 ~ polygons, raster...@data, col=rainbow(length(levels(rastermap $polygons str(raster.polygons) HTH Tom Hengl http://spatial-analyst.net -Original Message- From: r-sig-geo-boun...@stat.math.ethz.ch [mailto:r-sig-geo-boun...@stat.math.ethz.ch] On Behalf Of Tomislav Hengl Sent: Wednesday, January 28, 2009 9:56 AM To: r-sig-geo@stat.math.ethz.ch Cc: alexander.h...@csiro.au Subject: Re: [R-sig-Geo] calculate raster values based on vector regions Dear Herry, If I understand what you problem, one solution is to use R+SAGA. You should first convert the polygon map to the same grid, and then you can load it to R and do any type of aggregation: library(maptools) library(rgdal) library(RSAGA) # load the gridded map: rastermap - readGDAL(rastermap.asc) # load the shapefile: rsaga.esri.to.sgrd(in.grids= rastermap.asc, out.sgrds= rastermap.sgrd, in.path=getwd()) # convert the polygon map to a raster map: cellsize - raster...@grid@cellsize[1] rsaga.geoprocessor(lib=grid_gridding, module=3, param=list(GRID=polygons.sgrd, INPUT=polygons.shp, FIELD=1, LINE_TYPE=0, USER_CELL_SIZE=cellsize, user_x_extent_min=raster...@bbox[1,1]+cellsize, user_x_extent_max=raster...@bbox[1,2]-cellsize, user_y_extent_min=raster...@bbox[2,1]+cellsize, user_y_extent_max=raster...@bbox[2,2]-cellsize)) rsaga.sgrd.to.esri(in.sgrds=polygons.sgrd, out.grids=polygons.asc, out.path=getwd(), prec=0) rastermap$polygons - as.factor(readGDAL(polygons.asc)) # summary statistics per polygon class: raster.polygons - boxplot(band1 ~ polygons, rastermap @data, col=rainbow(length(levels(negotingrid$SOIL str(raster.polygons) You will of course need to adjust the FIELD position in the attribute table and the precision prec. To obtain and use SAGA, check this article: http://spatial-analyst.net/wiki/index.php?title=Software cheers, How is the weather now in Canberra? In Amsterdam is freezing brrr. T. Hengl -Original Message- From: r-sig-geo-boun...@stat.math.ethz.ch [mailto:r-sig-geo-boun...@stat.math.ethz.ch] On Behalf Of alexander.h...@csiro.au Sent: Wednesday, January 28, 2009 2:10 AM To: r-sig-geo@stat.math.ethz.ch Subject: [R-sig-Geo] calculate raster values based on vector regions Hi List, How would I go about assigning values to a raster based on polygon regions in R. Examples would be most appreciated. I have (vector) regions assigned to a specific value. The raster has NAs and some pixels where these values are likely to occur on ground. I need to assign these values to the raster and can do this in ArcInfo through vectors and rasterizing but only to a limited raster size. And R is much more preferable anyway... Any help appreciated. Cheers Herry ___ 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 mailing list R-sig-Geo@stat.math.ethz.ch https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] FW: Interpolcation option: IDW or OK?
Dear Yong Li, I hope you will not mind me joining this interesting discussion. If there is no evident spatial auto-correlation structure (pure nugget effect), IDW/OK are as good as randomly drawing a value from the global (normal) distribution. You can even test this using cross-validation! In principle, there is no justification to use distance-based interpolators if there is no evident spatial auto-correlation structure (maybe only the moving-window kriging, as implemented in e.g. Vesper, or stratified kriging techniques could discover some local spatial dependence). In addition, IDW should be considered an outdated technique, applicable only for situations where the variogram is close to linear (e.g. elevation data and similar smooth surfaces). What you should really consider using are the globaly available free maps/images (e.g. MODIS EVI, SRTM DEM parameters etc.), and then see if you can explain some of the variability in your target variable. But there will always be situations (especially in DSM applications) where you simply can not explain much of the target variability, neither with auxiliary maps nor with spatial auto-correlation. What to do then? I guess you simply have to collect more samples / more auxiliary maps and then try again. HTH T. Hengl See also: Compendium of Global datasets: http://spatial-analyst.net/wiki/index.php?title=Global_datasets Regression-kriging: http://spatial-analyst.net/wiki/index.php?title=Regression-kriging Pebesma, E., 2006. The Role of External Variables and GIS Databases in Geostatistical Analysis. Transactions in GIS, 10(4): 615-632. http://dx.doi.org/10./j.1467-9671.2006.01015.x -Original Message- From: r-sig-geo-boun...@stat.math.ethz.ch [mailto:r-sig-geo-boun...@stat.math.ethz.ch] On Behalf Of Edzer Pebesma Sent: Monday, February 09, 2009 9:08 AM To: Yong Li Cc: r-sig-geo@stat.math.ethz.ch Subject: Re: [R-sig-Geo] FW: Interpolcation option: IDW or OK? Yong Li wrote: Hi Edzer, I would say the spatial structure is regarded not significant when c0/c0+c1 is very much greater than 75%. In my case I used even distance intervals and calculated c0/c0+c1 for log(OLSENP) greater than 85%. I knew this index sometimes is very fragile, very much depending on how we fit the model. However when I zoomed in by using variable distance intervals (boundaries=c(100,200,300,400,600,900,1000,1500,2000))and maxdist=2000 meters, I found a pretty good model-fitted experimental variogram. But the local OK interpolation using such a fitted model did not make sense when compared the predictions to the observations as in most areas values of OLSENP were severely underestimated. You may have seen my code with which I have tried the nested models, but unfortunately no luck either. I maybe think the parameters for local ordinary kriging are not optimized, but I have tried lots of sets of nmin, nmax and maxdist and did see the hopeful end. The journal editor insists in OK being better than IDW. I need to collect my evidence to defend my IDW choice. That is my intention raised such a question in our forum here. I cannot find evidence in your data for such a claim; the cross validation statistics (rmse) seem to favour OK with your nested model. In your first email, you stated the following: Normally if we do not find a significant spatial structure for a soil variable, we may choose IDW or other methods. What is the argumentation behind this? Who claimed this? -- 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.springer.com/978-0-387-78170-9 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 mailing list R-sig-Geo@stat.math.ethz.ch https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] Calculating descriptive stats from many maps
Dear Rainer,
This is of course possible in R, and can be done in several ways:
1) for example, you can derive the average value using the rowSums function:
maps$Nsum - rowSums(m...@data, na.rm=T, dims=1)
maps$avg - maps$Nsum/(length(names(meuse.g...@data))-1)
You could also loop the sd, mean and quantile function over a range of cells:
for(i in length(names(m...@data))) {
m...@data$sd[i] - sd(m...@data[i,])
m...@data$mean[i] - mean(m...@data[i,])
...
}
This could take a lot of time!
2) if your maps are rather large, try also using the SAGA function:
rsaga.get.usage(lib = geostatistics_grid, module=5)
SAGA CMD 2.0.3
library path: C:/Progra~1/saga_vc/modules
library name: geostatistics_grid
module name : Statistics for Grids
This is probably the fastest method you can use.
HTH
T. Hengl
-Original Message-
From: r-sig-geo-boun...@stat.math.ethz.ch
[mailto:r-sig-geo-boun...@stat.math.ethz.ch] On Behalf
Of Rainer M Krug
Sent: Monday, February 09, 2009 4:58 PM
To: R-sig-Geo@stat.math.ethz.ch
Subject: [R-sig-Geo] Calculating descriptive stats from many maps
Hi
I have 25000 maps, generated by simulation predictions, covering the
same area, and would like to calculate some descriptive stats, like
mean, standard deviation, median, quartiles of all cells, to create a
variability map.
Is there an easy way of doing this in R?
Thanks,
Rainer
--
Rainer M. Krug, PhD (Conservation Ecology, SUN), MSc (Conservation
Biology, UCT), Dipl. Phys. (Germany)
Centre of Excellence for Invasion Biology
Faculty of Science
Natural Sciences Building
Private Bag X1
University of Stellenbosch
Matieland 7602
South Africa
___
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] Wiki?
Jonathan, I think this is an excellent idea. Ideally we should have a system that automatically scans e-mails, distinguishes between a question, answer/suggestion, and R code, and then generates a wiki article (such articles could then be classified into some hierachical structure using various similarity measures; e.g. based on the e-mail thread, users, libraries used etc.). Otherwise, nothing stops you from sorting the R-sig-geo traffic manually, extracting R codes and putting them into wiki articles. :) Tom Hengl http://spatial-analyst.net -Original Message- From: r-sig-geo-boun...@stat.math.ethz.ch [mailto:r-sig-geo-boun...@stat.math.ethz.ch] On Behalf Of Jonathan Greenberg Sent: Friday, February 20, 2009 5:39 AM To: r-sig-geo@stat.math.ethz.ch Subject: [R-sig-Geo] Wiki? R-sig-geo'ers -- do we have a wiki/script repository? If not, can we start a r-sig-geo wiki subpage on one of the main R wiki sites (http://wiki.r-project.org/rwiki/doku.php)? It would be nice to have a centralized location to paste help and/or scripts that we develop. --j -- Jonathan A. Greenberg, PhD Postdoctoral Scholar Center for Spatial Technologies and Remote Sensing (CSTARS) University of California, Davis One Shields Avenue The Barn, Room 250N Davis, CA 95616 Cell: 415-794-5043 AIM: jgrn307, MSN: jgrn...@hotmail.com, Gchat: jgrn307 ___ 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] Feature Extraction - object based image analysis
Hi Barry,
SAGA has few simple but useful modules for extraction of features from grids.
See:
library(RSAGA)
rsaga.get.modules(grid_discretisation)
$grid_discretisation
code name interactive
10 Supervised Classification FALSE
21 Cluster Analysis for Grids FALSE
32 Grid Segmentation FALSE
43 Grid Segmentation (b) FALSE
54 Grid Skeletonization FALSE
6 NA NA FALSE
7 NA NA FALSE
For multiple images, you can use either cluster analysis or supervised
classification:
rsaga.get.usage(grid_discretisation, 0)
SAGA CMD 2.0.3
library path: C:/Progra~1/saga_vc/modules
library name: grid_discretisation
module name : Supervised Classification
Usage: 0 -GRIDS str -POLYGONS str [-FIELD num] -CLASSES str -RESULT
str [-ML_PROB str]
[-METHOD num] [-NORMALISE] [-ML_THRESHOLD str]
-GRIDS:str Grids
Grid list (input)
-POLYGONS:str Training Areas
Shapes (input)
-FIELD:num Class Identifier
Table field
-CLASSES:strClass Information
Table (output)
-RESULT:str Classification
Grid (output)
-ML_PROB:strDistance/Probability
Grid (optional output)
-METHOD:num Method
Choice
Available Choices:
[0] Minimum Distance
[1] Maximum Likelihood
-NORMALISENormalise
Boolean
-ML_THRESHOLD:str Probability Threshold (Percent)
Floating point
Value Range: 0.00 - 100.00
You can also consider getting the statistics from grids for polygons:
rsaga.get.modules(shapes_grid)
$shapes_grid
code name interactive
1 0Add Grid Values to Points FALSE
2 1 Get Grid Data for Shapes FALSE
3 2 Grid Statistics for Polygons FALSE
4 3Grid Values to Points FALSE
5 4 Grid Values to Points (randomly) FALSE
6 5 Contour Lines from Grid FALSE
7 6 Vectorising Grid Classes FALSE
8 7 Clip Grid with Polygon FALSE
9 8 Gradient from Grid FALSE
10 NA NA FALSE
11 NA NA FALSE
You will soon discover that there is not much literature that explains how to
run processing and
what does each parameter in SAGA means :(
Unfortunately, Olaf does not have time to write up a user's manual for SAGA,
but you can always zoom
into the original code. I am sure that all algorithms come from some
literature. E.g.:
Kothe, R., Bock, M., 2006. Development and use in practice of SAGA modules for
quality analysis of
geodata. In Böhner, J., McCloy, K.R., Strobl, J. [Eds.]: SAGA – Analysis and
Modelling Applications.
Göttinger Geographische Abhandlungen, Vol.115.
http://surfnet.dl.sourceforge.net/sourceforge/saga-gis/gga115_08.pdf
HTH
Tom Hengl
http://spatial-analyst.net
-Original Message-
From: r-sig-geo-boun...@stat.math.ethz.ch
[mailto:r-sig-geo-boun...@stat.math.ethz.ch] On Behalf Of
Rawlins, Barry G
Sent: Friday, February 27, 2009 3:06 PM
To: r-sig-geo@stat.math.ethz.ch
Subject: [R-sig-Geo] Feature Extraction - object based image analysis
I have a colleague who is not familiar with R and we were discussing feature
extraction from
multiple images (air photos, satellite based etc.). He uses a commercial tool
for feature
extraction. Can anyone point me to a package in R that does this? I have done a
search but not found
such a function. It is not a procedure I am familiar with.
Many thanks, Barry
Dr Barry Rawlins
Sustainable Soils Team Leader
British Geological Survey
Keyworth
Nottingham
NG12 5GG
Direct tel: 0115 9363140
Mobile: 0788 4235473
homepage:
http://barryrawlins.googlepages.comhttp://barryrawlins.googlepages.com/
--
This message (and any attachments) is for the recipient ...{{dropped:8}}
___
R-sig-Geo mailing list
R-sig-Geo@stat.math.ethz.ch
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
[R-sig-Geo] GEOSTAT 2009 Summer School, 3-10 May 2009 MEDILS, Split, Croatia 2nd call
Dear R-sig-geo, There are still some open places for our GEOSTAT 2009 summer school (the registration closes on 15h of March). We will finally be three lecturers: (1) Roger Bivand; (2) Olaf Conrad (the main creator of SAGA) and (3) Tomislav Hengl. The fifth day of the summer school, we will run a R-sig-geo LIVE discussion - via a live broadcast (under preparation). Current registrations: uni-bonn.de; zg.htnet.hr; gmx.de; izor.hr; agr.hr; ugent.be; geo.uzh.ch; geodata.hr; uclouvain.be; iastate.edu; ceh.ac.uk; F.A.Q. 1. Where can I find more info about the course programme? - The programme of the school will be published after we collect the registrations. We keep some flexibility considering the programme because we try to adjust to the topics that applicants select as much as possible. 2. Is there a possibility to get a discount or funding from the hosting institute? - In principle no. We have tried our best to reduce the price as much as possible (all the lecturers are volunteers), while still being able to use excellent facilities (not to mention that MEDILS has its private beach). Possibly the only way to save some money is to try to organize the accommodation yourself. The registration costs for external accommodation will be 300-350 € (but you would soon find out that accommodation in Split is as expensive as on the institute). 3. Can two or more people from the same institute apply? - Of course! The selection criteria are only: (1) scientific excellence, (2) relevance, (3) distance from the course (more distant applications have advantage), (4) time of application. my apologies for cross posting! 2nd call: GEOSTAT 2009 Summer School Spatio-temporal data analysis with R + SAGA + Google Earth Period: 3 – 10 May 2009 Location: Mediterranean Institute For Life Sciences (MEDILS), Split, Croatia The objective of the GEOSTAT Summer School is to promote use of open source tools for the analysis of spatio-temporal data. The course aims at providing a balanced combination of theoretical and hands-on-software (R + open source GIS) training, and includes a one-day workshop where each participant can present his/her case study and receive feedback. Course moderators: * Roger S. Bivand (Norwegian School of Economics and Business Administration) * Tomislav Hengl (University of Amsterdam) * Olaf Conrad (University of Hamburg) The MEDILS institute features fully equipped lecture and seminar rooms. Accommodation is available in-house (hotel complex on the upper floors) - comfortable two-bed rooms with private bathrooms. All rooms are fully air-conditioned and provide telephone and high-speed Internet access. Three meals per day (breakfast, lunch, dinner) and intermediate coffee/snacks are provided at the Institute. The institute also features sport facilities, and access to sea (beach). Registration fees: Single rooms*: 840 € Shared rooms (2 beds): 680 € *A limited number of single rooms is available. In the case of high interest for single rooms, the organizers will try to arrange an equivalent accommodation in the neighboring hotel. The registration fees are ALL INCLUSIVE: they include workshop fees, accommodation at the institute (six nights), meals during the whole duration of the summer school (excluding Friday 8.05 dinner and Sunday 10.05 lunch), and costs of the Saturday 9.05 excursion to island Brac (Pustinja Blaca - old Monastery). Registration dead-line: 15 March 2009 The course is limited to 25 participants. The candidates will be selected based on: (a) academic excellent; (b) thematic relevance to the topic of the course; (c) geographic location (more distant applications have advantage) and (d) time of application. To register for this summer school, please fill out and forward the registration form available at: http://www.medils.org/index.php/archives/category/events/schools-workshops/ For all other information: Tomislav Hengl, course coordinator E-mail: t.he...@uva.nl URL: http://spatial-analyst.net ___ 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] residuals by subtracting point values from grid values
Dear Else, I've read your message twice, but I am still not sure what exactly you want to run: (1) which prediction technique, (2) for which purpose (validation only?). Using the krige.cv method, you can also run the cross validation for inverse distance interpolation (note: you only need to omit the variogram): idw.zn.cv - krige.cv(zinc~1, meuse, nfold=2) [inverse distance weighted interpolation] ... This is all explained in Edzer's paper on gstat: Pebesma, E.J. 2004. Multivariable geostatistics in S: The gstat package. Computers and Geosciences 30: 683–91. http://dx.doi.org/10.1016/j.cageo.2004.03.012 see also pages 222-226 in Bivand et al. (2008) code available at: http://www.asdar-book.org/book/geos_mod.R here are some similar simple examples: http://spatial-analyst.net/wiki/index.php?title=Spatial_prediction_of_soil_moisture HTH Tom Hengl -Original Message- From: r-sig-geo-boun...@stat.math.ethz.ch [mailto:r-sig-geo-boun...@stat.math.ethz.ch] On Behalf Of Els Verfaillie Sent: Tuesday, March 24, 2009 3:00 PM To: r-sig-geo@stat.math.ethz.ch Subject: [R-sig-Geo] residuals by subtracting point values from grid values Hi list, I've splitted my dataset (e.g. the meuse dataset) into random 2 parts (2/3 for prediction and 1/3 for validation). I've done an interpolation (e.g. IDW, but I also want to do it with different kriging techniques) based on the prediction dataset. Now I want to subtract the zinc values from the validation dataset (i.e. a SpatialPointsDataFrame with 52 points) from the var1.pred values from the interpolation (i.e. a SpatialPointsDataFrame with 3103 points) to obtain the residuals. This of course only on the location of the 52 points from the validation dataset. Then I want to calculate validation indices. This is my code: #Validation and prediction dataset permid - sample(155) idportion - 1:round(2*155/3) meuse.pred - meuse[permid[idportion], ] meuse.val - meuse[permid[-idportion], ] #idw based on meuse.pred idw.pred.zn - idw(zinc ~ 1, meuse.pred, meuse.grid) I know that I could use the krige.cv to do an n-fold cross-validation, but I would like to do it based on the above method. Any suggestions? Thanks a lot! Els __ Dr. Els Verfaillie Carto-GIS cluster Ghent University (UGent) - Department of Geography Krijgslaan 281 - S8 B- 9000 Gent Belgium __ [[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] Running analysis on DEMs using GRASS (from R?)
Dear list, I am trying to test some DEM analysis functions in GRASS GIS (under MS Windows machines). This is my first contact with GRASS scripting (I have successfully installed GRASS 6.4 for windows from http://grass.itc.it/grass63/binary/mswindows/native/). For example, with SAGA GIS, I can completely control the analysis from R and use SAGA as external application to run the processing, e.g. to extract the drainage networks from a DEM: library(rgdal) library(RSAGA) # obtain the data: download.file(http://www.geomorphometry.org/data/DEM25m.zip;, destfile=paste(getwd(),DEM25m.zip,sep=/)) fname - zip.file.extract(file=DEM25m.asc, zipname=DEM25m.zip) file.copy(fname, ./DEM25m.asc, overwrite=TRUE) # load to SAGA format: rsaga.esri.to.sgrd(in.grids=dem25m.asc, out.sgrds=dem25m.sgrd, in.path=getwd()) # Fill the spurious sinks: rsaga.geoprocessor(lib=ta_preprocessor, module=2, param=list(DEM=dem25m.sgrd, RESULT=dem25m_f.sgrd, MINSLOPE=0.05)) # extract the channel network: rsaga.geoprocessor(lib=ta_channels, module=0, param=list(ELEVATION=dem25m_f.sgrd, CHNLNTWRK=channel_ntwrk.sgrd, CHNLROUTE=channel_route.sgrd, SHAPES=channels.shp, INIT_GRID=dem25m_f.sgrd, DIV_CELLS=3, MINLEN=40)) # read back to R: channels - readOGR(channels.shp, channels) spplot(channels[LENGTH], col.regions=bpy.colors()) How could I run the same processing using GRASS? Can you control GRASS from R at all? (I would be a bit disappointed with GRASS if I am not able to run the analysis from R). Many thanks for your time and suggestions in advance! Tom Hengl http://home.medewerker.uva.nl/t.hengl/ ___ 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] Running analysis on DEMs using GRASS (from R?)
Dear Augustin,
Thank you for your message and examples shown (flavored with a bits of
Portuguese). I am aware of
the spgrass6 package. I also know that there are plenty of examples for the R
GRASS connection e.g.
http://casoilresource.lawr.ucdavis.edu/drupal/node/438, but I would still like
to see how an R
script would look like for the case study shown below
(http:/www.geomorphometry.org/data/DEM25m.zip; extraction of channel
network).
I am still confused whether GRASS can be run from R (or is it only the other
way around?)? how to
import the data into GRASS (e.g. ArcInfo ASCII)? how to set-up the parameters
of the GRASS
environment (working directory, projection system etc)? how to control the
processing from R (e.g.
using system(...)), and then read the results of analysis back to R?
Please excuse my ignorance (and the fact that I am a MS Windows user). I am
just a beginner with
GRASS, so it could be that things are much simpler than they seem (it could
also be that they are
also much more complicated).
Many thanks,
Tom Hengl
http://home.medewerker.uva.nl/t.hengl/
-Original Message-
From: Agustin Lobo [mailto:alobolis...@gmail.com]
Sent: Wednesday, March 25, 2009 10:46 AM
To: Tomislav Hengl
Cc: r-sig-geo@stat.math.ethz.ch; carlos.grohm...@gmail.com
Subject: Re: [R-sig-Geo] Running analysis on DEMs using GRASS (from R?)
Tomislav,
This would be part of an R script using grass:
brachspread - function(itmax=1,a=-1)
{
#Do the next once per session
#attach(/media/mifat32/Rutils/utiles.rda)
#Execute grass command: import tiff file into grass
#Data:
#manchas (importado de )
#noocup (importado de no_ocupable1.tif, 3==no ocupable)
#(senderos importado de senderos_raster10m.tif)
#trailbuf importado de senderos
#r.buffer -z --o in=senderos out=trailbuf distance=10 units=meters
#camino:1; contiguo:2; resto, NULL
require(spgrass6) #Package connecting R-grass
G - gmeta6() #Read grass environment
system(g.region rast=manchas)
system(r.mapcalc 'habqual = if(trailbuf==1,0,trailbuf)')
system(r.mapcalc 'habqual = if(habqual==2,0.8,habqual)')
system(r.mapcalc 'habqual = if(isnull(habqual),0.5,habqual)')
#Execute grass command: set region to raster
system(g.region rast=manchas)
system(r.mapcalc 'manchasit = manchas')
#inicio del bucle de itersciones tiempo
it -0
system(r.mapcalc 'prob = 0.000')
for(it in 1:itmax) {
system(r.clump --o in=manchasit out=manchasit.clmp)
#Execute grass command: create stats file
system(r.stats -cl in=manchasit.clmp out=delme.txt)
#Read stats file into R:
mistat - read.table(delme.txt,header=F,comment.char=*)
system(rm delme.txt) #delete garbage
#Write reclass file for grass
write.table(mistat, file=delme.txt, row.names=F, col.names=F ,sep==)
system(r.reclass --o in=manchasit.clmp out=manchas.sze rules=delme.txt)
system(rm delme.txt)
etc
Tomislav Hengl wrote:
Dear list,
I am trying to test some DEM analysis functions in GRASS GIS (under MS
Windows machines). This is
my
first contact with GRASS scripting (I have successfully installed GRASS 6.4
for windows from
http://grass.itc.it/grass63/binary/mswindows/native/).
For example, with SAGA GIS, I can completely control the analysis from R and
use SAGA as external
application to run the processing, e.g. to extract the drainage networks from
a DEM:
library(rgdal)
library(RSAGA)
# obtain the data:
download.file(http://www.geomorphometry.org/data/DEM25m.zip;,
destfile=paste(getwd(),DEM25m.zip,sep=/))
fname - zip.file.extract(file=DEM25m.asc, zipname=DEM25m.zip)
file.copy(fname, ./DEM25m.asc, overwrite=TRUE)
# load to SAGA format:
rsaga.esri.to.sgrd(in.grids=dem25m.asc, out.sgrds=dem25m.sgrd,
in.path=getwd())
# Fill the spurious sinks:
rsaga.geoprocessor(lib=ta_preprocessor, module=2,
param=list(DEM=dem25m.sgrd,
RESULT=dem25m_f.sgrd, MINSLOPE=0.05))
# extract the channel network:
rsaga.geoprocessor(lib=ta_channels, module=0,
param=list(ELEVATION=dem25m_f.sgrd,
CHNLNTWRK=channel_ntwrk.sgrd, CHNLROUTE=channel_route.sgrd,
SHAPES=channels.shp,
INIT_GRID=dem25m_f.sgrd, DIV_CELLS=3, MINLEN=40))
# read back to R:
channels - readOGR(channels.shp, channels)
spplot(channels[LENGTH], col.regions=bpy.colors())
How could I run the same processing using GRASS? Can you control GRASS from R
at all? (I would be
a
bit disappointed with GRASS if I am not able to run the analysis from R).
Many thanks for your time and suggestions in advance!
Tom Hengl
http://home.medewerker.uva.nl/t.hengl/
___
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] Running analysis on DEMs using GRASS (from R?)
Thanks to Roger (and with many useful comments from Augustin) GRASS just got another big user - possible even many more :) Roger has prepared an update of his package spgrass6 that can now be used to control the GRASS processing from R (functions: initGRASS, execGRASS). This way one can skip some manual settings of GRASS (environmental parameters) and just continue with the analysis (this is now really easy). These are the installation steps (MS Windows machines): 1. Download and install GRASS 6 (e.g. from http://grass.itc.it/grass63/binary/mswindows/native/) 2. Install the updated version of the spgrass6 package (from http://spatial.nhh.no/R/Devel/spgrass6_0.6-1.zip) 3. Start R and then run the following code: # obtain the data: download.file(http://geomorphometry.org/data/dem25m.zip;, destfile=paste(getwd(),dem25m.zip,sep=/)) fname - zip.file.extract(file=DEM25m.asc, zipname=dem25m.zip) file.copy(fname, ./DEM25m.asc, overwrite=TRUE) library(spgrass6) # Location of your GRASS installation: loc - initGRASS(C:/GRASS, home=tempdir()) loc # Import the ArcInfo ASCII file to GRASS: execGRASS(r.in.gdal, flags=o, parameters=list(input=DEM25m.asc, output=DEM)) execGRASS(g.region, parameters=list(rast=DEM)) gmeta6() # extract the drainage network: execGRASS(r.watershed, flags=overwrite, parameters=list(elevation=DEM, stream=stream, threshold=as.integer(50))) # thin the raster map so it can be converted to vectors: execGRASS(r.thin, parameters=list(input=stream, output=streamt)) # convert to vectors: execGRASS(r.to.vect, parameters=list(input=streamt, output=streamt, feature=line)) streamt - readVECT6(streamt) plot(streamt) Here are some screen shots: http://geomorphometry.org/R.asp Tom Hengl http://spatial-analyst.net -Original Message- From: Roger Bivand [mailto:roger.biv...@nhh.no] Sent: Wednesday, March 25, 2009 1:27 PM To: Agustin Lobo Cc: Tomislav Hengl; carlos.grohm...@gmail.com; r-sig-geo@stat.math.ethz.ch Subject: Re: [R-sig-Geo] Running analysis on DEMs using GRASS (from R?) On Wed, 25 Mar 2009, Agustin Lobo wrote: Tomislav, You set up your grass location and mapset, import your data etc. outside R. Then, from the grass console, you start R and run an scrip like the one I sent you or tits individual commands. Once you've done require(spgrass6) #Package connecting R-grass G - gmeta6() #Read grass environment you can run grass commands using system(), pass results from R to grass and from grass to R (It's done in the script I sent you), import raster files (which I normally don't do: raster files are too large for R in general, I prefer to use grass for processing and then pass data as tables to R for the stats; but have no experience with the new raster package made by Roger The raster package on R-forge is by Robert Hijmans and a team of developers - very interesting, and often the most active project on the whole (very active) site! As you write, spgrass6 itself uses system() within its functions. A year ago I did some work on an automatic interface generator for GRASS commands in R using the interface descriptions provided by GRASS modules in XML. Would there be interest in using these to try to automate the writing of GRASS commands from R? Would anyone like to help? Should spgrass6 be moved from sourceforge CVS to R-forge SVN to make it easier for others to join in? Roger Actually, I often use R as an scripting language for grass. I used to do scripts in the Bourne shell itself (actually, it was csh), but using R (perhaps python also) is much easier (and powerful). I do not see that many differences between the way you are using saga and the way you would use grass, except that for grass you must have your location, mapset and region set. This is actually very powerful, as (re)defining regions with g.region lets you restrict your analysis to an area of interest, which you can change by repeated calls to g.region that could depend upon a dynamic process. The same can be done by setting a MASK. The script that I sent you is part of a process simulating the spread of an invasive plant. Note that all control is within R. You can do the same with R alone, but if your raster layers are large, it could become very slow (note that there is for() loop in the script, and for() loops have a very low performance in R; as most of the commands run within the loop are grass commands run through system(), there is minimum memory built up). I've used the grass that comes within QGIS in windows sometimes, but I feel really not comfortable doing this type of work outside linux though. You seem experienced with saga. Actually, would like to talk to you on a possible use of sage from within qgis thanks to the fact that both saga and qgis make use of python. Agus Tomislav Hengl wrote: Dear Augustin, Thank you for your message
Re: [R-sig-Geo] Spatial Interpolation of Regularly Gridded Data
Hi Greg,
Yes, there are many possibilities for downscaling grids in R, so you are at the
right place. :)
1. If you only wish to downscale climatic grids (e.g. using splines), then
probably the most
efficient (fastest; can handle large grids) way is to use the downscaling
method in SAGA:
library(RSAGA)
library(RCurl)
out.url -
getURL(https://stat.ethz.ch/pipermail/r-sig-geo/attachments/20090406/6b460881/attachment.obj;,
ssl.verifypeer=F)
write(out.url, out.csv)
out - read.csv(out.csv, header=T, sep=,)
str(out)
coordinates(out) - ~lon+lat
proj4string(out) - CRS(+proj=longlat +datum=WGS84)
writeOGR(out, ., out, ESRI Shapefile)
# convert to a grid:
rsaga.geoprocessor(lib=grid_gridding, module=3,
param=list(GRID=out5deg.sgrd, INPUT=out.shp,
FIELD=2, LINE_TYPE=0, USER_CELL_SIZE=0.5, USER_FIT_EXTENT=T))
# downscale the grid using Bicubic spline interpolation:
rsaga.geoprocessor(lib=grid_tools, module=0,
param=list(INPUT=out5deg.sgrd,
GRID=out25deg.sgrd, METHOD=0, DIMENSIONS_CELLSIZE=0.25, SCALE_DOWN_METHOD=3))
2. Otherwise, I would really suggest that you try downscaling grids using
auxiliary information /
geostatistical modelling, which is explained in e.g. the following two papers:
Hengl, T., Bajat, B., Reuter, H.I., Blagojevic, D., 2008. Geostatistical
modelling of topography
using auxiliary maps. Computers Geosciences, 34: 1886-1899.
http://geomorphometry.org/view_scripts.asp?id=6
Grohmann, C.H., Steiner, S.S., 2008. SRTM resample with Short Distance-Low
Nugget Kriging.
International Journal of Geographical Information Science, 22 (8):895-906.
http://dx.doi.org/10.1080/13658810701730152
HTH,
T. Hengl
-Original Message-
From: r-sig-geo-boun...@stat.math.ethz.ch
[mailto:r-sig-geo-boun...@stat.math.ethz.ch] On Behalf
Of Greg King
Sent: Tuesday, April 07, 2009 12:01 AM
To: r-sig-geo@stat.math.ethz.ch
Subject: [R-sig-Geo] Spatial Interpolation of Regularly Gridded Data
Hi,
I'm fairly new to R, but finding the experience a good one. However I am a
little overwhelmed by
the number of packages and not sure I am necessarily using the most
appropriate ones.
Here is the background to what I am trying to achieve: I have a CSV file
which contains weather
forecast data for latitude and longitude points (see attached out.csv, the
data I understand is on
WGS84 http://www.ready.noaa.gov/faq/geodatums.html). The sample points are
at half degree
intervals. My objective is to work out what the forecast data is at any
specific given
latitude/longitude by interpolating data from the 0.5x0.5 degree grid. I am
doing this for a
number of different time points using the following functions:
library(akima)
library(rgdal)
gks - function(inLat,inLon,dframe,variab) {
wind.grid - expand.grid(lat=inLat, lon=inLon)
coordinates(wind.grid) - ~ lat+lon
proj4string(wind.grid) - CRS(+init=epsg:4326)
pnt-interpp(coordinates(dframe)[,1], coordinates(dframe)[,2],
z=as.data.frame(dframe)[,1],
coordinates(wind.grid)[,1],coordinates(wind.grid)[,2],linear=TRUE)
return(pnt$z)
}
interp_gk - function(lat, lon) {
wind-read.csv(file=/Users/greg/Out.csv,head=TRUE,sep=,,
colClasses=c(POSIXct,numeric,numeric,numeric,numeric))
coordinates(wind)=~lat+lon
proj4string(wind) - CRS(+init=epsg:4326)
times-unique(wind$vt)
columns-names(wind)[2:length(names(wind))]
dOut-data.frame(dateTime=times[1])
for (i in 1:length(times)) {
dOut[i,dateTime]-times[i]
for (j in 1:length(columns)) {
sset-subset(wind, wind$vt==times[i], select=c(columns[j]))
dOut[i,columns[j]]-gks(lat,lon,sset,columns[j])
}
}
dOut-cbind(dOut, mag=sqrt(dOut$ugrd_0^2+dOut$vgrd_0^2))
return(dOut)
}
However, I have the following concerns:
* Should I really be using akima? The documentation states it is not for
use on regularly
spaced grids - what are my alternatives?
* The interp funcrtion will not work for cubic spline linear=FALSE
interpolation (is it
because my data is regularly gridded?). How can I achieve cubic spline
interpolation?
* Is my function really using the Cordinate reference system specified?
When I comment out
the CRS lines, my functions return the same values?
Lots of questions I appreciate, but I am curious! It seems R can achieve
what I am trying to
do... but I may just be missing some vital bits of information.
Thanks,
Greg
___
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] problem to module MERGE GRID in RSAGA
-Original Message-
From: r-sig-geo-boun...@stat.math.ethz.ch
[mailto:r-sig-geo-boun...@stat.math.ethz.ch] On Behalf
Of gianni lavaredo
Sent: Tuesday, April 07, 2009 4:44 PM
To: r-sig-geo@stat.math.ethz.ch
Subject: [R-sig-Geo] problem to module MERGE GRID in RSAGA
Dear R Users
I have a problem to understand the loop merge in RSAGA
I have this file (not order name)
DEM2644147
DEM2644148
DEM2644149
DEM2744150
DEM2744155
DEM2744168
I did a file list in CSV:
path- (C:/mydata/)
temp - paste(path,input.csv, sep=)
input - read.csv(temp, header = TRUE)
attach(input)
names(input)
temp - file
file.list - as.list(temp)
# USE POWER TO DEM
power - seq(1.0, 10.0, by=0.1)
# power of DCM
for(i in 1:length(file.list)){
for (method in 1:length(power)){
rsaga.grid.calculus(in.grids = c(paste(C:/mydata/DEM,
file.list[[i]],.sgrd, sep=)),
out.grid =
c(paste(C:/mydata/Power_DEM,file.list[[i]],_,method+9,.sgrd,sep=)),
formula = paste(a^,power[method],sep=))
}
}
#TRY with two DEM with power^1
rsaga.geoprocessor(lib=grid_tools, module=3,
param=list(GRIDS=C:/mydata/Power_DEM2644147_10.sgrd;C:/mydata/Power_DEM2644148_10.sgrd,
MERGED=C:/mydata/merge_power10.sgrd,TYPE=6,INTERPOL=0,OVERLAP=0))
Now I am trying to merge in a LOOP every DEM to have:
all_merge_power10
all_merge_power11
all_merge_power12
all_merge_power13
etc
for(i in 1:length(file.list)){
for (method in 1:length(power)){
rsaga.geoprocessor(lib=grid_tools, module=3,
param=list(GRIDS=c(paste(C:/mydata/Power_DEM,file.list[[i]],method+9,.sgrd,sep=)),
MERGED=C:/mydata/all_merge_power10.sgrd,
TYPE=6,
INTERPOL=0,
OVERLAP=0))
}
}
why do you loop?
simly do:
GRIDS=paste(set.file.extension(file.list, .sgrd), collapse=;)
There is also no need to write the working directory every time. Just set it up
at the beginning:
setwd(C:/mydata/)
Some similar examples here:
http://spatial-analyst.net/wiki/index.php?title=Software#SAGA_GIS
T. Hengl
I have diffucult to load all DEM power 10, ect etc and every time I have a
Error :-(.
Thanks a lot
Gianni
[[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] Converting downloaded PNG images using RgoogleMaps to SpatialGridDataFrame class (proj4string?)
Hi all, Thanks to Markus Loecher and colleagues we can now easily obtain background maps from Google Earth and use it for plotting/interpretation of spatial data. This runs pretty smoothly (e.g. a map of the Netherlands): library(RgoogleMaps) # obtain the API key and save into the home folder; MyMap - GetMap.bbox(center=c(52.1551723,5.3872035), zoom=7, destfile=netherlands.png, maptype =hybrid) Read 1 item [1] http://maps.google.com/staticmap?center=52.1551723,5.3872035zoom=7size=640x640 + maptype=hybridformat=png32key=sensor=true trying URL 'http://maps.google.com/staticmap?center=52.1551723,5.3872035zoom=7size=640x640 + maptype=hybridformat=png32key==true' Content type 'image/png' length 703541 bytes (687 Kb) opened URL downloaded 687 Kb netherlands.png has GDAL driver PNG and has 640 rows and 640 columns PlotOnStaticMap(MyMap, lat=52.1551723, lon=5.3872035) But how to convert this image into a SpatialGridDataFrame to allow overlays, export to GIS formats etc? I am not even sure what the correct proj4string for this image is? str(MyMap) # this only gives coordinates of the centre List of 4 $ : num 52.2 $ : num 5.39 $ : num 7 $ : num [1:640, 1:640] 187 60 66 1 66 202 151 1 214 220 ... ..- attr(*, COL)= chr [1:256] #7C #887440 #747838 #0C1C20 ... ..- attr(*, type)= chr rgb GDALinfo(netherlands.png) # the coordinates are not embedded! Thanks, Tom Hengl http://home.medewerker.uva.nl/t.hengl/ ___ R-sig-Geo mailing list R-sig-Geo@stat.math.ethz.ch https://stat.ethz.ch/mailman/listinfo/r-sig-geo
[R-sig-Geo] Info: SAGA vs GRASS workshop, University of Zurich 29 August 30 August 2009
Title: Automated analysis of elevation data in R+SAGA/GRASS Venue: University of Zurich, Irchel Campus, 29 August 30 August 2009 Workshop moderators: Tomislav Hengl / Carlos H. Grohmann Registration fees: 150 CHF (PhD students) For more info: http://2009.geomorphometry.org/ Summary: This workshop aims at PhD students and professionals interested to use open source software packages for processing of their elevation data. R is the open-source version of the S language for statistical computing; SAGA (System for Automated Geoscientific Analyses) and GRASS (Geographic Resources Analysis Support System) are the two most used open-source desktop GIS for automated analysis of elevation data. A combination of R+SAGA/GRASS provides a full integration of statistics and geomorphometry. The topics in this workshop will range from selection of grid cell size, choice of algorithms for DEM generation and filtering, to geostatistical simulations and error propagation. The workshop moderators will demonstrate that R+SAGA/GRASS is capable of handling such demanding tasks as DEM generation from auxiliary maps, automated classification of landforms, and sub-grid parameterization of surface models. The course will focus on understanding R and SAGA/GRASS syntax and building scripts that can be used to automate DEM-data processing. Each participant should come with a laptop PC and install all software needed prior to the workshop. Registered participants will receive an USB stick with all data sets and overheads at the beginning of the course. Participants will follow a case study that focuses on generation of DEMs, extraction of DEM parameters and landform classes, and implementation of error propagation in geomorphometry. To register, fill in and forward the registrations forms available at: http://2009.geomorphometry.org/registration.htm ___ R-sig-Geo mailing list R-sig-Geo@stat.math.ethz.ch https://stat.ethz.ch/mailman/listinfo/r-sig-geo
[R-sig-Geo] From the GEOSTAT 2009 Summer School for PhD students (report)
FYI: The official report from the GEOSTAT 2009 Summer School for PhD students: http://spatial-accuracy.org/FromGEOSTAT2009 TITLE: Spatio-temporal data analysis with R+SAGA+Google Earth LOCATION: Mediterranean Institute for Life Sciences, Split, Croatia PERIOD: 3-10 May 2009 LECTURERS: R. Bivand, T. Hengl and O. Conrad I have also started maintaining a small list of the upcoming training courses (R and spatial data analysis) and people involved: http://spatial-analyst.net/wiki/index.php?title=Training_in_R Let me know if you think that I have missed some important event/institution/person. Tom Hengl http://spatial-analyst.net (Maybe R-sig-geo should consider starting a newsletter?) ___ 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] Problems with loading the mgcv package
Hi Edzer, I am sorry about that. I have been using exclusively Tinn-R to write R scripts now over a year (once you install Tinn-R over R, you do need to set-up the Rprofile; Tinn-R works only with the MDI mode of Rgui). I have just installed a somewhat older version of R (2.7) and tried loading the package spatstat in a clean session and it does load without a problem. I assume that the problem might be with my Rprofile and/or the Tinn-R integration. Tinn-R does need to load some packages: library(TinnR) library(svSocket) which might conflict with mgcv. I did not have problems with Tinn-R with the latest version of R/Tinn-R, so it is a new thing. I will try to contact also the Tinn-R people (http://www.sciviews.org/Tinn-R/). thanks, Tom Hengl http://home.medewerker.uva.nl/t.hengl/ -Original Message- From: Edzer Pebesma [mailto:edzer.pebe...@uni-muenster.de] Sent: Thursday, June 11, 2009 11:02 AM To: Tomislav Hengl Cc: r-sig-geo@stat.math.ethz.ch Subject: Re: [R-sig-Geo] Problems with loading the mgcv package In general, before you post questions of this kind to the list, it makes sense to see if the packages load on your platform in a clean session, e.g. started with R --vanilla. (Or, anyhow, without automatically loading a data base or Rprofile type modified local configurations). Then you know where to search! -- Edzer Tomislav Hengl wrote: Hi guys, Thanks for your reply (and apologies for a delay - it was night time in Europe). I am running Windows XP, R version 2.9.0. I have just reinstalled spatstat with dependencies: install.packages(spatstat, dependencies=TRUE) Then I tried: library(spatstat) and I get the same error: Loading required package: mgcv Error in formatDL(nm, txt, indent = max(nchar(nm, w)) + 3) : ... then I tried to re-install the mgcv package but this gives me an error: Warning: cannot remove prior installation of package 'mgcv' *this is version 1.5-5 When I look under my \library\mgcv\ I can still see a dll file which can not be removed. Maybe there is a problem that I am using a customized Rprofile? My Rprofile is available at: http://spatial-analyst.net/scripts/Rprofile.site Thanx! Tom Hengl http://home.medewerker.uva.nl/t.hengl/ -Original Message- From: Dylan Beaudette [mailto:debeaude...@ucdavis.edu] Sent: Wednesday, June 10, 2009 7:40 PM To: r-sig-geo@stat.math.ethz.ch Cc: Hengl, T. Subject: Re: [R-sig-Geo] Problems with loading the mgcv package works for me: library(spatstat) sessionInfo() Loading required package: mgcv This is mgcv 1.5-5 . For overview type `help(mgcv-package)'. Loading required package: gpclib General Polygon Clipper Library for R (version 1.4-4) Type 'class ? gpc.poly' for help Loading required package: deldir deldir 0.0-8 spatstat 1.15-3 R version 2.9.0 (2009-04-17) i686-pc-linux-gnu locale: LC_CTYPE=en_US.UTF-8;LC_NUMERIC=C;LC_TIME=en_US.UTF-8;LC_COLLATE=en_US.UTF- 8;LC_MONETARY=C;LC_MESSAGES=en_US.UTF-8;LC_PAPER=e n_US.UTF-8;LC_NAME=C;LC_ADDRESS=C;LC_TELEPHONE=C;LC_MEASUREMENT=en_US.UTF-8;LC_IDENTIFICATION=C attached base packages: [1] stats graphics grDevices utils datasets methods base other attached packages: [1] spatstat_1.15-3 deldir_0.0-8gpclib_1.4-4mgcv_1.5-5 loaded via a namespace (and not attached): [1] grid_2.9.0 lattice_0.17-25 nlme_3.1-92 On Wednesday 10 June 2009, Marcelino de la Cruz wrote: Hi Tom, I can load sptastat and mgcv without problem (see below). Could you please send us information about your environment? Cheers, Marcelino library(spatstat) Loading required package: mgcv This is mgcv 1.5-5 . For overview type `help(mgcv-package)'. Loading required package: gpclib General Polygon Clipper Library for R (version 1.4-3) Type 'class ? gpc.poly' for help Loading required package: deldir deldir 0.0-8 spatstat 1.15-3 Type 'help(spatstat)' for information sessionInfo() R version 2.9.0 (2009-04-17) i386-pc-mingw32 locale: LC_COLLATE=Spanish_Spain.1252;LC_CTYPE=Spanish_Spain.1252;LC_MONETARY=Spani sh_Spain.1252;LC_NUMERIC=C;LC_TIME=Spanish_Spain.1252 attached base packages: [1] stats graphics grDevices utils datasets methods base other attached packages: [1] spatstat_1.15-3 deldir_0.0-8gpclib_1.4-3mgcv_1.5-5 loaded via a namespace (and not attached): [1] grid_2.9.0 lattice_0.17-22 nlme_3.1-90 tools_2.9.0 At 16:40 10/06/2009, Hengl, T. wrote: Hi, I am having problems with loading the mgcv package, which is used by the spatstat package. I did not have similar problems with the latest version of the spatstat, so I imagine that the problem could lay with the mgcv package (new parameters). library(spatstat) Loading required
Re: [R-sig-Geo] Slope and Aspect calculations in R
Yes: SAGA runs normally on Linux (I use it also), but you can not use the
command line from R (i.e.
RSAGA) to pass command to SAGA :(
But you can use R to pass command lines to SAGA without RSAGA (i.e. the batch
files):
http://saga-gis.wiki.sourceforge.net/Executing+Modules+with+SAGA+CMD
Tom Hengl
http://spatial-analyst.net
-Original Message-
From: Thomas Adams [mailto:thomas.ad...@noaa.gov]
Sent: Friday, July 10, 2009 4:20 PM
To: Tomislav Hengl
Cc: r-sig-geo@stat.math.ethz.ch
Subject: Re: [R-sig-Geo] Slope and Aspect calculations in R
Tom,
Thanks! I have seen your posts and results of some of your work with
SAGA and I am very interested. But, I have been to lazy to this point to
try compiling SAGA for either MacOS X or Linux (is a Linux distribution
available?).
Regards,
Tom
Tomislav Hengl wrote:
Dear Josh, Thomas,
I typically use R+SAGA (R on top!) for such type of analysis (and I do have
access to all
operating
system/software). There must be some good reasons that I am using SAGA
(speed, stability, rich
geomorphometry module etc.)! So do not give up so easy on SAGA.
I especially like that SAGA is fairly efficient in processing big grids,
and the commands are
rather
compact. For example to reproject grids from one to another proj4string,
you only need to set
e.g.:
rsaga.geoprocessor(lib=pj_proj4, 2, param=list(SOURCE_PROJ=paste('',
proj4string(dem25m),
'',
sep=), TARGET_PROJ=\+proj=longlat +datum=WGS84\, SOURCE=TWI.sgrd,
TARGET=TWI_ll.sgrd,
TARGET_TYPE=0, INTERPOLATION=0))
*see http://www.geomorphometry.org/R.asp
I remember that somebody did manage to set-up SAGA also under Mac OS
(http://www.saga-gis.uni-goettingen.de/html/index.php?module=pnForumfunc=viewtopictopic=458),
but
I also remember that most of students on our GEOSTAT course used Windows
under Mac OS to run
SAGA. I
am sure that it should be possible to compile SAGA under Mac OS, but
somebody needs to put some
effort to bridge this gap (check also with Olaf Conrad).
Here is more discussion about the same topic from few months ago:
https://stat.ethz.ch/pipermail/r-sig-geo/2009-April/005418.html
cheers,
Tom Hengl
-Original Message-
From: r-sig-geo-boun...@stat.math.ethz.ch
[mailto:r-sig-geo-boun...@stat.math.ethz.ch] On
Behalf
Of Thomas Adams
Sent: Friday, July 10, 2009 3:12 PM
To: Josh London
Cc: r-sig-geo@stat.math.ethz.ch
Subject: Re: [R-sig-Geo] Slope and Aspect calculations in R
Josh,
I use R GRASS together on both MacOS X and Linux without problems.
It's very direct to write a shell script that calls both R GRASS
within a single script. Writing a R script that calls GRASS may be more
tedious to do what you want. Writing a shell script that calls both may
be the wat to go.
Tom
Josh London wrote:
Thanks for the suggestions so far regarding the use of GRASS or SAGA.
We have considered those options but SAGA doesn't appear available for
MacOS. My understanding is that spgrass6 is more for accessing R
code/functions from within GRASS vs accessing GRASS functions from
within R. We have all our other analysis coded up in R, and just need
these functions to do one part of a bigger process. Maybe I'm
misunderstanding how to use R and GRASS together. I'll do some further
research.
Thanks again
Josh
Josh London wrote:
Hello
We are looking to mimic, in R, the slope (max magnitude difference
between a cell and its neighbors) and aspect (direction of maximum
magnitude difference) functions found in ESRI's Spatial Analyst
package. We were just about to code up the slope function but thought
I would make sure we weren't re-inventing any wheels already out there.
Thanks
Josh
--
Thomas E Adams
National Weather Service
Ohio River Forecast Center
1901 South State Route 134
Wilmington, OH 45177
EMAIL: thomas.ad...@noaa.gov
VOICE: 937-383-0528
FAX: 937-383-0033
___
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
--
Thomas E Adams
National Weather Service
Ohio River Forecast Center
1901 South State Route 134
Wilmington, OH 45177
EMAIL:thomas.ad...@noaa.gov
VOICE:937-383-0528
FAX: 937-383-0033
___
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] more doubts with spTransform()
Hi Agustin, You should really work with some referent point for which you know both local and geographic (WGS84) coordinates. Here is an illustration: http://spatial-analyst.net/wiki/index.php?title=MGI_/_Balkans_coordinate_systems#Validation_of_CRS_p arameters Note that I use Google Earth imagery to validate that everything is fine. I agree - setting up the right CRS can be a headache, but there are many 'hard' (historic, technology connected etc) reasons for this. The good thing is that: once you got it right, you soon forget about all these problems. HTH, T. Hengl -Original Message- From: r-sig-geo-boun...@stat.math.ethz.ch [mailto:r-sig-geo-boun...@stat.math.ethz.ch] On Behalf Of Agustin Lobo Sent: Tuesday, July 14, 2009 9:43 AM To: sig-geo Subject: [R-sig-Geo] more doubts with spTransform() Inverse transform now: anyellaWGS84 - readOGR(dsn=., layer=parcellesAnyellaWGS84) OGR data source with driver: ESRI Shapefile Source: ., layer: parcellesAnyellaWGS84 with 6 rows and 7 columns Feature type: wkbPoint with 2 dimensions proj4string(anyellaWGS84) [1] +proj=utm +zone=31 +ellps=WGS84 +datum=WGS84 +units=m +no_defs +towgs84=0,0,0 anyellaED50notowgs84 - spTransform(anyellaWGS84,CRS(+proj=utm +zone=31 +ellps=intl +units=m +no_defs)) anyellaED50towgs84 - spTransform(anyellaWGS84,CRS(+proj=utm +zone=31 +ellps=intl +units=m +no_defs +towgs84=-87,-98,-121)) coordinates(anyellaWGS84)[1,] coords.x1 coords.x2 421005 4683849 coordinates(anyellaED50notowgs84)[1,] coords.x1 coords.x2 421001.4 4683932.4 coordinates(anyellaED50towgs84)[1,] coords.x1 coords.x2 421097.5 4684050.9 which is correct? anyellaED50notowgs84 or anyellaED50towgs84 ? I guess anyellaED50towgs84, but would like to confirm, Thanks Agus ___ 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] getting subfolders in a directory on FTP (MODIS)
Hey Barry,
Your code works fine (at least with the NASA ftps). Thanks!
I knew it - everything is possible in R :)
T. Hengl
http://spatial-analyst.net
-Original Message-
From: b.rowling...@googlemail.com [mailto:b.rowling...@googlemail.com] On
Behalf Of Barry
Rowlingson
Sent: Tuesday, July 14, 2009 4:53 PM
To: Tomislav Hengl
Cc: alexandre villers; r-sig-geo@stat.math.ethz.ch
Subject: Re: [R-sig-Geo] getting subfolders in a directory on FTP (MODIS)
On Tue, Jul 14, 2009 at 2:12 PM, Tomislav
Henglhe...@spatial-analyst.net wrote:
To tell you honestly, I could not figure out how to fetch names of
sub-folders from an ftp
(maybe it
is simply not possible?),
Apparently the FTP spec doesn't specify how to return a list of
directories, which, I read, is a source of headaches to FTP client
developers.
If you do:
resp = getURL(ftp://n4ftl01u.ecs.nasa.gov/SAN/MOST/MOD10CM.005/;)
items = strsplit(resp,\n)[[1]]
you get the folders (and files) but the folder names are in the form
of a unix directory listing. So you then want the last word of any
lines that start with 'd'. Lets get the lines first:
folderLines = items[substr(items,1,1)=='d']
that gets you all the lines. There's probably a neater way of getting
the last word than this:
lastBit - function(x){x[length(x)]}
dirs = unlist(lapply(strsplit(folderLines, ),lastBit))
Note this only works on the FTP server specified, other FTP servers
are at liberty to send back any dir format they want (I think I was
surprised the first time I connected to a DOS-based FTP server and saw
a DOS DIR listing!).
Barry
___
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] still struggling with FTP download on MODIS
-Original Message-
From: r-sig-geo-boun...@stat.math.ethz.ch
[mailto:r-sig-geo-boun...@stat.math.ethz.ch] On Behalf
Of Alexandre VILLERS
Sent: Wednesday, July 15, 2009 10:03 AM
To: Aide R SIG GEO
Subject: [R-sig-Geo] still struggling with FTP download on MODIS
Dear all,
I'm still not sure that this list is the right place for that. So if
those questions bother you, please give me the name of the list I should
contact.
Here a version of the script I'm using to download MODIS data from NSIDC
(weekly snow cover)
Tsnow -
ftp://anonymous:t...@n4ftl01u.ecs.nasa.gov/SAN/MOST/MOD10A2.005/;
###server of NSIDC with directories / id and password
myfolder-paste(tempdir(),\\ , sep=) #the folder to download data
#-getting the list of the folders in Tsnow (thanks to Barry Rowlingson)
resp -getURL(Tsnow)
items - strsplit(resp,\r\n)[[1]]
folderLines - items[substr(items,1,1)=='d']
lastBit - function(x){x[length(x)]}
dates - unlist(lapply(strsplit(folderLines, ),lastBit))
datesK-data.frame(dates=dates, keep=rep(NA,length(dates)))
# select months for which I want the data
months-c(01,02,03,04,05,06,10, 11, 12) #months for
which I want the data
for (i in months){ #loop to keep only the dates I'm interested in
garde-grep(datesK$dates, pattern=paste(.,i,., sep=))
datesK$keep[garde]-1
}
datesK-subset(datesK, datesK$keep==1) # keep the dates of interest
for (i in 1:length(datesK[,1])){## loop to download the tiles of
interest
getlist1 - strsplit(getURL(paste(Tsnow, datesK[i,][[1]], /,
sep=),.opts=curlOptions(ftplistonly=TRUE)), \r\n)[[1]] #get the
files in the folder i
BLOCKtemp - getlist1[grep(getlist1,
pattern=MOD10A2..h18v04.005.*.hdf)[1]] #select the
tile of interest (in my case h18v04.005)
download.file(paste(Tsnow, datesK[i,][[1]], /,
BLOCKtemp,sep=), destfile=paste(myfolder, BLOCKtemp, sep=),
mode='wb', method='wget') #download the tile
}
So here is my problem. I can't get into some specifc folders and receive
the error message.
Error in curlPerform(curl = curl, .opts = opts, .encoding = .encoding) :
server did not report OK, got 550
The ftp error code 550 means: Requested action not taken. File unavailable
(e.g., file not found,
no access)., which means that the ftp server is not allowing you access.
Did you try checking your wget tool? I remember I had headaches before I manage
to get the wget
running without problems. Make sure you set:
options(download.file.method=auto)
wget.exe you need to put in your Windows directory folder (Note: make sure you
disable your
antivirus tools such as Norton or McAfee otherwise it might block wget from
running);
For example, i=4 (folder 2000.03.13). And this is always true,
suggesting it's not a problem of connection. The error appears as soon
as getURL is called in the loop.
I succeed in downloading for some values of i (manually) but when
looping, it stops on those black folders. I've checked directly from
filezilla or my web browser: I can get into those folder and download
the tiles.
The tiles are pretty small (150 to 450 Ko). So if someone could try from
his place to see if he gets the same error message for the same value of
i (4, 28,31, etc.).
Since I have 550 tiles to download and process, I would prefer not to do
that manually (and I can't get to use the tool proposed by the NSIDC for
the moment).
Thanks for the help
Alex
P.S.: my next question will deal with some projection problem !!!
--
Alexandre Villers
PhD Student
Team Biodiversity
Centre d'Etudes Biologiques de Chizé-CNRS UPR1934
79360 Beauvoir sur Niort
Phone +33 (0)5 49 09 96 13
Fax +33 (0)5 49 09 65 26
__ Information from ESET Mail Security, version of virus signature
database 4244
(20090715) __
The message was checked by ESET Mail Security.
http://www.eset.com
___
R-sig-Geo mailing list
R-sig-Geo@stat.math.ethz.ch
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
T. Hengl
http://spatial-analyst.net
___
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] Problem with size of dataset
Dear Emmanuel, I have the same problem. I either can not run processing with large data set in R or I can not even load such data to R. Then, if I want to do any geostatistics, it takes forever. R (gstat/geoR) is simply not that efficient with large spatial data as e.g. GIS software. What you can definitively try is to subset your point data randomly by using e.g.: points.sub - points[runif(length(poi...@data[[1]]))0.1] # 10% sample! This will allow you to fit variograms etc. Then, if you really want to interpolate all of your 157k points, you might consider using SAGA. For example, you can also sub-set randomly points from a shape file: # too many points; subset to 5% (Split Shapes Layer Randomly in SAGA): rsaga.geoprocessor(lib=shapes_tools, module=17, param=list(SHAPES=complete.shp, A=part_A.shp, B=part_B.shp, PERCENT=5)) ground.sample - readOGR(part_B.shp,part_B) # Learn more about geostatistics in SAGA: rsaga.get.modules(lib=geostatistics_kriging) $geostatistics_kriging code name 1 0 Ordinary Kriging (VF) 2 1 Ordinary Kriging (VF, Global) 3 2 Universal Kriging (VF) 4 3 Universal Kriging (VF, Global) 5 4Semivariogram (Dialog)) 6 5 Ordinary Kriging 7 6 Ordinary Kriging (Global) 8 7 Universal Kriging 9 8 Universal Kriging (Global) 10 NA NA 11 NA NA # Read the mask map: gridmaps - readGDAL(maskmap.asc) cell.size - gridm...@grid@cellsize[[1]] # Ordinary kriging in SAGA: rsaga.geoprocessor(lib=geostatistics_kriging, module=5, param=list(GRID=var_OK.sgrd, SHAPES=var.shp, BVARIANCE=F, BLOCK=F, FIELD=1, BLOG=F, MODEL=1, TARGET=0, NPOINTS_MIN=10, NPOINTS_MAX=60, NUGGET=rvgm.Pb$psill[1], SILL=1.65, RANGE=1238, MAXRADIUS=5, USER_CELL_SIZE=cell.size, user_x_extent_min=gridm...@bbox[1,1]+cell.size/2, user_x_extent_max=gridm...@bbox[1,2]-cell.size/2, user_y_extent_min=gridm...@bbox[2,1]+cell.size/2, user_y_extent_max=gridm...@bbox[2,2]-cell.size/2)) # the same way you can run regression-kriging/universal kriging; You will soon find out that there is a big difference in the efficiency between SAGA and R - SAGA will interpolate your 157k points within few minutes or less. On the other hand, SAGA has very very limited geostatistical functionality (for example it can not fit variograms etc.), so what you really need is a combination of SAGA and R! Here are more examples: http://geomorphometry.org/view_scripts.asp?id=24 HTH, T. Hengl http://spatial-analyst.net -Original Message- From: r-sig-geo-boun...@stat.math.ethz.ch [mailto:r-sig-geo-boun...@stat.math.ethz.ch] On Behalf Of Poizot Emmanuel Sent: Thursday, July 16, 2009 9:55 AM To: r-sig-geo@stat.math.ethz.ch Subject: [R-sig-Geo] Problem with size of dataset Dear all, I would like to perform a geostatistical analysis using R. To do so, I'm classicaly use geoR ou GSTAT packages. In the present case, I have a dataset of around 157000 locations where I have a value (depth). I've been not able to create both geodata or gstat valid R object because apparently of the size of the dataset. DOes anybody have an idear of how to conduct such study with such dataset size ? Regards -- Cordialement Emmanuel Poizot Cnam/Intechmer B.P. 324 50103 Cherbourg Cedex Phone (Direct) : (00 33)(0)233887342 Fax : (00 33)(0)233887339 ___ 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] krige.bayes and data export
Hi Laura, Sorry for a bit delay in reply. I hope you did not give up on geoR. The truth is: transforming result of prediction using geoR to e.g. sp classes is not trivial (it looks like as creators of geoR did not really consider this option). Once you get it into sp class, then it is easy to export it using rgdal or to a database format (?). I had difficulties myself to get it running - but it is doable. There are certainly many benefits of using geoR (actually, I am starting to like it more and more). data(meuse) coordinates(meuse) - ~x+y zinc.geo - as.geodata(meuse.ov[zinc]) zinc.svar2 - variog(zinc.geo, lambda=0, max.dist=1500, messages=FALSE) zinc.vgm2 - likfit(zinc.geo, lambda=0, messages=FALSE, ini=c(var(log1p(zinc.geo$data)),500), cov.model=exponential) # prepare prediction locations: locs - pred_grid(c(pc.co...@bbox[1,1]+gridcell/2, pc.co...@bbox[1,2]-gridcell/2), c(pc.co...@bbox[2,1]+gridcell/2, pc.co...@bbox[2,2]-gridcell/2), by=gridcell) ... mask.bor - m...@polygons[[1]]@polygons[[...@coords ... zinc.ok2 - krige.conv(zinc.geo, locations=locs, krige=krige.control(obj.m=zinc.vgm2), borders=mask.bor) ... mask.ov - overlay(mask, locs.sp) mask.sel - !is.na(mask.ov$MASK.SGRD) locs.geo - data.frame(x=locs...@coords[mask.sel,1], y=locs...@coords[mask.sel,2], zinc.ok2=zinc.ok2[[1]], zinc.rkvar2=zinc.ok2[[2]]) coordinates(locs.geo) -~X+Y gridded(locs.geo) - TRUE write.asciigrid(locs.geo[1], zinc_ok2.asc, na.value=-1) Here are more examples (see section Geostatistics in geoR): http://spatial-analyst.net/scripts/meuse.R T. Hengl -Original Message- From: r-sig-geo-boun...@stat.math.ethz.ch [mailto:r-sig-geo-boun...@stat.math.ethz.ch] On Behalf Of Laura Poggio Sent: Wednesday, July 22, 2009 6:11 PM To: r-sig-geo@stat.math.ethz.ch Subject: [R-sig-Geo] krige.bayes and data export Dear list, I got some simulations with krige.bayes of geoR packages. I can plot them with the code provided in the example. However I also would like to export the results of the simulations as database. This in order to access them later for further analysis. To plot them I access the simulations as following: image(awc.res.bayes, val= simulation, number.col=1) image(awc.res.bayes, val= simulation, number.col=2) image(awc.res.bayes, val= simulation, number.col=3) image(awc.res.bayes, val= simulation, number.col=...) image(awc.res.bayes, val= simulation, number.col=30) I did not find a way to store the simulation values in a dataframe. And I did not find any example around. I hope I manage to explain my problem. Thank you very much in advance Laura Poggio [[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
Re: [R-sig-Geo] Polygon to raster?
SAGA has a module to convert from shapes to grids (GRASS can do the same). See: library(RSAGA) rsaga.get.usage(lib=grid_gridding, module=0) SAGA CMD 2.0.3 library path: C:/Progra~1/saga_vc/modules library name: grid_gridding module name : Shapes to Grid Usage: 0 [-GRID str] -INPUT str [-FIELD num] [-TARGET_TYPE num] [-LINE_TYPE num] [-USER_CELL_SIZE str] [-USER_FIT_EXTENT] [-USER_X_EXTENT_MIN str] [-USER_X_EXTENT_MAX str] [-USER_Y_EXTENT_MIN str] [-USER_Y_EXTENT_MAX str] [-USER_GRID_TYPE num] [-GET_SYSTEM_SYSTEM_NX num] [-GET_SYSTEM_SYSTEM_NY num] [-GET_SYSTEM_SYSTEM_X str] [-GET_SYSTEM_SYSTEM_Y str] [-GET_SYSTEM_SYSTEM_D str] [-GET_SYSTEM_GRID_TYPE num] [-GRID_GRID str] -GRID:str Grid Data Object (optional output) -INPUT:str Shapes Shapes (input) -FIELD:num Attribute ... Here is an example: ... # first, convert the contour map to a raster map: rsaga.geoprocessor(lib=grid_gridding, module=3, param=list(GRID=contours_buff.sgrd, + INPUT=z_contours.shp, FIELD=1, LINE_TYPE=0, USER_CELL_SIZE=dem.pixelsize, + user_x_extent_min=elevati...@bbox[1,1], user_x_extent_max=elevati...@bbox[1,2], + user_y_extent_min=elevati...@bbox[2,1], user_y_extent_max=elevati...@bbox[2,2])) # now extract a buffer distance map (for 'contours') and load it back to R: # the parameters DIST and IVAL are estimated based on the grid properties: rsaga.geoprocessor(lib=grid_tools, module=10, param=list(SOURCE=contours_buff.sgrd, + DISTANCE=contours_dist.sgrd, ALLOC=contours_alloc.sgrd, BUFFER=contours_bdist.sgrd, + DIST=sqrt(dem.area)/3, IVAL=dem.pixelsize)) rsaga.sgrd.to.esri(in.sgrds=contours_dist.sgrd, out.grids=contours_dist.asc, +out.path=getwd(), prec=1) contours.dist - readGDAL(contours_dist.asc) spplot(contours.dist[1], col.regions=bpy.colors(), scales=list(draw=TRUE), +sp.layout=list(sp.lines, col=cyan, z.contours)) hist(contours.dist$band1, col=grey) ... HTH T. Hengl http://spatial-analyst.net/wiki/index.php?title=Grid_size_calculator -Original Message- From: r-sig-geo-boun...@stat.math.ethz.ch [mailto:r-sig-geo-boun...@stat.math.ethz.ch] On Behalf Of Jonathan Greenberg Sent: Sunday, August 02, 2009 8:26 PM To: r-sig-geo@stat.math.ethz.ch Subject: [R-sig-Geo] Polygon to raster? I'm trying to hunt down a function that will take a polygon (derived from a shapefile) and rasterize it at a given pixel resolution -- any hints? --j -- Jonathan A. Greenberg, PhD Postdoctoral Scholar Center for Spatial Technologies and Remote Sensing (CSTARS) University of California, Davis One Shields Avenue The Barn, Room 250N Davis, CA 95616 Cell: 415-794-5043 AIM: jgrn307, MSN: jgrn...@hotmail.com, Gchat: jgrn307 ___ 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] A working draft of lecture notes A Practical Guide to Geostatistical Mapping
Dear R-sig-geo, I have just finished with updating the second edition of my lecture notes that I have used over years for the GEOSTAT summer schools for PhD students. A copy of the working draft can be obtained here: http://spatial-analyst.net/book/ I would like to invite you to read selected chapters (and review the attached R scripts) and then post comments via the editorial system (please read the instructions for reviewers first!). As a reward for your effort, I plan to send some 80 free copies to the most active reviewers (based on the web traffic), so please also make sure you enter a correct mailing address. FYI all datasets used in this book have been put on-line: http://spatial-analyst.net/DATA/ http://spatial-analyst.net/worldmaps/ I sincerely plan to keep this publication of open access. This was a difficult decision, because open access inevitably carries a risk of lower quality and lower confidence in what is said. Hopefully, you will also see the benefits of open access publishing and help me improve this book by sending me your comments and suggestions. This way it will reach some minimum quality standards that will make us both happy. cheers, T. Hengl ___ 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] Overlay bubble chart on map
Here is another example with overlay from the Reimann et al. (2008) book (this uses a raster map as background - this if often more impressive): http://www.statistik.tuwien.ac.at/StatDA/R-scripts/page151.html The code is rather long, but it will give you good ideas what you have to do. HTH T. Hengl On Sun, Aug 16, 2009 at 12:20 PM, Rebecca Rossrebecca.r...@plants.ox.ac.uk wrote: Dear All, I am new to this list, so apologies if this question has been asked and answered many times before! I would like to overlay a bubble chart onto a map so that the bubbles correspond to particular points on the map. I have the shape file for the map and the lat long co-ordinates for the points. Could anyone tell me how to put the two together? The suggestions I found so far are plot.Map and plot.Spatial, neither of which work. I am working to a bit of a deadline, so a response today would be really helpful. To read foo.shp: install.packages(rgdal) require(rgdal) map = readOGR(/path/to/maps/,foo) plot(map) then you can do: points(p$x,p$y) or whatever to overlay on it. If a 'bubble chart' is just points scaled by some factor, you can either use the 'symbols' function in R or use pch and cex options to plot - with something like: p = data.frame(x=1:10,y=1:10,size=1:10) plot(p$x,p$y,pch=19,cex=p$size) hope that helps! Barry ___ 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] Producing ENFA derived suitability maps with adehabitat
Dear Alexandre, I think that this should be doable. Once you run the ENFA to estimate the habitat model, copy the values of new rasters (future habitat) and replace the original values in the object. Then simply re-run the predict.enfa method e.g.: # prepare the dataset for ENFA: beidata - data2enfa(as.kasc(list(dem=import.asc(dem.asc), grad=import.asc(grad.asc), twi=import.asc(twi.asc), achan=import.asc(achan.asc))), bei.sub@coords) # run ENFA: enfa.bei - enfa(dudi.pca(beidata$tab, scannf=FALSE), beidata$pr, scannf=FALSE, nf=2) # same area, new environmental conditions: beidata.new - data2enfa(as.kasc(list(dem=import.asc(dem_new.asc), grad=import.asc(grad_new.asc), twi=import.asc(twi_new.asc), achan=import.asc(achan_new.asc))), bei.sub@coords) # needs to have the same mask! enfa.bei.new - enfa(dudi.pca(beidata.new$tab, scannf=FALSE), beidata.new$pr, scannf=FALSE, nf=2) # copy the model parameters fitted previously: enfa.bei.new$m - enfa.bei$m enfa.bei.new$s - enfa.bei$s enfa.bei.new$lw - enfa.bei$lw enfa.bei.new$li - enfa.bei$li enfa.bei.new$co - enfa.bei$co enfa.bei.new$mar - enfa.bei$mar bei.dist.new - predict(enfa.bei.new, beidata.new$index, beidata.new$attr) ...something like this (I am not really sure which parameters do you have to copy, but the $tab data.frame certainly needs to be replaced with new predictors). Let me know if it works. cheers, T. Hengl http://spatial-analyst.net/wiki/index.php?title=Species_Distribution_Modelling -Original Message- From: r-sig-geo-boun...@stat.math.ethz.ch [mailto:r-sig-geo-boun...@stat.math.ethz.ch] On Behalf Of Alexandre VILLERS Sent: Monday, September 07, 2009 9:31 AM To: Aide R SIG GEO; r-sig-ecol...@r-project.org Subject: [R-sig-Geo] Producing ENFA derived suitability maps with adehabitat Dear list(s) 'members, I have a dataset representing the position of a species over a large region for 3 different years (2000, 2004 and 2008) and seven ecogeographical variables. Given that the study takes place on an agricultural region, the landscape changes every year at fine scale and there is a long term trend in crops sowed, leading me to account for between years variability. I would like to use the niche determined by the landscape and location of birds in 2000 and then, appply this niche over the landscape in 2004 and 2008 (this would, I believe, give me a first answer on whether changes in birds' location are due to a decrease in habitat suitability or birds). I have already computed ENFA with adehabitat (using the doc provided with the package adehabitat and the www.spatial-analyst.net of T. Hengl) but I don't see how exactly using the result of an ENFA on a new landscape... Any link or help would be welcome. Alex -- Alexandre Villers PhD Student Team Biodiversity Centre d'Etudes Biologiques de Chizé-CNRS UPR1934 79360 Beauvoir sur Niort Phone +33 (0)5 49 09 96 13 Fax +33 (0)5 49 09 65 26 __ Information from ESET Mail Security, version of virus signature database 4401 (20090906) __ The message was checked by ESET Mail Security. http://www.eset.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 mailing list R-sig-Geo@stat.math.ethz.ch https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] Producing ENFA derived suitability maps with adehabitat
Dear Mathieu,
Thanks for the clarification. I have underestimated the complexity, but I
admit - habitat analysis is not really 'my cup of tee' (maybe this
question was more suited for the Environmentrics R-sig).
As least I was correct about the whole concept - it is doable ;)
I would be interested to learn from Alexandre how did the prediction of
future habitat work and does it satisfies their objectives.
T. Hengl
http://home.medewerker.uva.nl/t.hengl/
Dear all,
The answer Tom gave was actually incorrect. Let me first explain why,
and then give a solution to the problem.
The so-called predictions of the ENFA (function 'predict.enfa') are
based on the row coordinates ($li) and the vector of presence ($pr).
From the row coordinates, the function computes Mahalanobis distances
from the center of the niche to every pixel (row), given the niche
covariance structure. So that if you copy the $li of the first ENFA
object to the second, considering that the vector of presence $pr
remains the same, then you will have exactly the same predictions, which
is not exactly what we want.
Now, let's try to give a solution to this problem. The row coordinates
gives you the coordinates of every pixel projected into the new space
created by the ENFA. You can compute them by hand by multiplying the
(scaled) table of environmental variables with the columns coordinates.
For example, following the example of the function 'enfa':
## We load the data
data(lynxjura)
map - lynxjura$map
tmp - lynxjura$locs[,4]!=D
locs - lynxjura$locs[tmp, c(X,Y)]
map[,4] - sqrt(map[,4])
## We perform the ENFA
dataenfa1 - data2enfa(map, locs)
pc - dudi.pca(dataenfa1$tab, scannf = FALSE)
(enfa1 - enfa(pc, dataenfa1$pr, scannf = FALSE))
## We compute the row coordinates by hand and check that this is
## correct
bla - as.matrix(enfa1$tab) %*% as.matrix(enfa1$co)
all.equal(bla, as.matrix(enfa1$li))
This means that we can then reproject any new environmental variable
into the space created by the ENFA. Let's say that a second set of
environmental maps is called 'map2'. It can be then projected as follows:
li2 - as.matrix(scale(kasc2df(map2)$tab)) %*% as.matrix(enfa1$co)
Note that:
1) the table must be scaled, as is the $tab component of an ENFA object
(modified array).
2) for the sake of simplicity, I used here the function 'scale', which
used a N-1 denominator. In the ENFA example, a PCA ('dudi.pca') is ran
first on the data, before we use the scale table. In the 'dudi.pca'
function, the denominator is N, so that it might give slightly different
results (negligible for large N).
Then, to predict, you must use an adapted 'predict.enfa' function,
with a 'new' argument in which we feed the new maps:
predict.enfa - function (object, index, attr, nf, new, ...)
{
if (!inherits(object, enfa))
stop(should be an object of class \enfa\)
warning(the enfa is not mathematically optimal for prediction:\n
please consider the madifa instead)
if ((missing(nf)) || (nf object$nf))
nf - object$nf
Zli - object$li[, 1:(nf + 1)]
f1 - function(x) rep(x, object$pr)
Sli - apply(Zli, 2, f1)
m - apply(Sli, 2, mean)
cov - t(as.matrix(Sli)) %*% as.matrix(Sli)/nrow(Sli)
if (!missing(new))
Zli - as.matrix(dudi.pca(kasc2df(new)$tab, scannf = FALSE,
nf = 1)$tab) %*% as.matrix(object$co)
maha - mahalanobis(Zli, center = m, cov = cov)
map - getkasc(df2kasc(data.frame(toto = maha, tutu = maha),
index, attr), toto)
return(invisible(map))
}
And this should work like this:
pred2 - predict(enfa1, dataenfa1$index, dataenfa1$attr, new = map2)
Note that the two maps have here the same attributes. It still works if
it's not the case, simply by adjusting the index and attr in the predict
call (given by a kasc2df on map2).
Hope this helps,
Mathieu.
Tomislav Hengl a écrit :
Dear Alexandre,
I think that this should be doable.
Once you run the ENFA to estimate the habitat model, copy the values of
new rasters (future habitat)
and replace the original values in the object. Then simply re-run the
predict.enfa method e.g.:
# prepare the dataset for ENFA:
beidata - data2enfa(as.kasc(list(dem=import.asc(dem.asc),
grad=import.asc(grad.asc),
twi=import.asc(twi.asc), achan=import.asc(achan.asc))),
bei.sub@coords)
# run ENFA:
enfa.bei - enfa(dudi.pca(beidata$tab, scannf=FALSE), beidata$pr,
scannf=FALSE, nf=2)
# same area, new environmental conditions:
beidata.new - data2enfa(as.kasc(list(dem=import.asc(dem_new.asc),
grad=import.asc(grad_new.asc), twi=import.asc(twi_new.asc),
achan=import.asc(achan_new.asc))),
bei.sub@coords) # needs to have the same mask!
enfa.bei.new - enfa(dudi.pca(beidata.new$tab, scannf=FALSE),
beidata.new$pr, scannf=FALSE, nf=2)
# copy the model parameters fitted previously:
enfa.bei.new$m - enfa.bei$m
enfa.bei.new$s - enfa.bei$s
enfa.bei.new$lw - enfa.bei$lw
Re: [R-sig-Geo] DEMs
Hi Abe,
Did you try using the ETOPO1? This (0.1 degree) DEM can be well loaded to R:
URL - http://spatial-analyst.net/worldmaps/;
# list of maps:
map.list - c(globedem, chlo08, dcoast)
# download the zipped maps one by one:
for(i in 1:length(map.list)) {
download.file(paste(URL, map.list[i], .zip, sep=),
+ destfile=paste(getwd(), /, map.list[i], .zip, sep=))
unzip(zipfile=paste(map.list[i], .zip, sep=), exdir=getwd())
unlink(paste(map.list[i], .zip, sep=))
}
# read maps to R:
worldmaps - readGDAL(paste(map.list[i], .asc, sep=))
# fix the layer name:
names(worldmaps)[1] - map.list[1]
for(i in map.list[-1]) {
worldm...@data[map.list[i]] - readGDAL(paste(map.list[i], .asc,
sep=))$band1
}
proj4string(meuse.grid) - CRS(+proj=longlat +ellps=WGS84)
I do not think that it will be easy to generate kriging with any global
dataset that is 5 km resolution. You could try using the kriging in SAGA
(it can handle maps up to 2 GB in size; it is exp(2) times faster than R).
Here is an example:
http://spatial-analyst.net/wiki/index.php?title=Interpolation_of_ISRIC-WISE_international_soil_profile_data
T. Hengl
PS: Which variable are you interpolating?
Hello All,
My questions are regarding the best way to obtain DEMs for the entire
globe. I am trying to krige a dataset of global scope. The dataset has
points at all latitude/longitude intersects and so there are 64,800 sample
points. The problem that I am having is that the DEMs I have been using
are
too large (memory) to be used as R variables. I also have had trouble
finding a way to download DEMs for the entire globe.
1. Given the size of my sample dataset: What is the maximum resolution of
DEM that I should be using? Do I only need values at prediction points or
will higher res yield higher accuracy?
2. Is there a way to programmatically access DEMs using GDAL or R in
general? ie. internally or through integration with McIDAS, NetCDF, etc.
3. Do I need to break the DEMs and sample dataset into smaller tiles to
avoid the memory issues? If so how do I do that in R / GDAL?
Thanks,
Abe
[[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
Re: [R-sig-Geo] DEMs
-Original Message-
From: Aben Woha [mailto:abenw...@gmail.com]
Sent: Tuesday, September 15, 2009 7:52 AM
To: Tomislav Hengl
Cc: r-sig-geo@stat.math.ethz.ch
Subject: Re: [R-sig-Geo] DEMs
Hi Tomislav,
Thanks for the detailed response. I tried to use the globedem.zip from your
website, but then
when I call overlay() with my data, there are many NAs produced in the
output: it seems the DEM
covers X= -180 to 180 Y= -65 to 65, my dataset covers X= -180 to 180 Y=
-90 to 90
So when I call krige, it complains:
Error: dimensions do not match: locations 129600 and data 47029
129600 = (number of sample points) * 2
47029 = (number of sample points) - (number of NAs in overlay dem data)
Any suggestions on how to handle this? Should I subset the data? Find a
different dem with full
global coverage?
You should mask the NA values prior to variogram fitting using e.g.:
sel - !is.na(point.ov$SOLAR)
plot(variogram(SOLAR~1, point.ov[sel,]))
SRTM DEM is available only for 56°S to 60°N. The 10 minutes global DEM (GTOPO?)
at is probably your
best choice then (it is 2160 x 900 pixels):
http://biogeo.berkeley.edu/worldclim1_4/grid/cur/alt_10m_bil.zip (metadata
missing :( )
Also when you say a global dataset 5km resolution, is that as measured at
the equator?
Yes. The ETOPO1 (10 km) DEM is 3600 x 1300 pixels, which is already on the edge
of what R can
handle. There is also the GDEM (ASTER DEM, which has global coverage), but
nobody yet produced a
single layer (e.g. at resolution of 1 km).
The current variable I am working with is solar radiation. I have plans to
work on several other
variables as well.
But do you need to produce maps of solar radiation for artics/antartica, sea
etc? Besides, plenty of
global solar radiation maps already exist e.g.:
1. http://www.meteonorm.com/pages/en/downloads/maps.php
2. http://re.jrc.ec.europa.eu/pvgis/countries/europe.htm
...
On Mon, Sep 14, 2009 at 1:58 AM, Tomislav Hengl he...@spatial-analyst.net
wrote:
Hi Abe,
Did you try using the ETOPO1? This (0.1 degree) DEM can be well loaded
to R:
URL - http://spatial-analyst.net/worldmaps/;
# list of maps:
map.list - c(globedem, chlo08, dcoast)
# download the zipped maps one by one:
for(i in 1:length(map.list)) {
download.file(paste(URL, map.list[i], .zip, sep=),
+ destfile=paste(getwd(), /, map.list[i], .zip, sep=))
unzip(zipfile=paste(map.list[i], .zip, sep=), exdir=getwd())
unlink(paste(map.list[i], .zip, sep=))
}
# read maps to R:
worldmaps - readGDAL(paste(map.list[i], .asc, sep=))
# fix the layer name:
names(worldmaps)[1] - map.list[1]
for(i in map.list[-1]) {
worldm...@data[map.list[i]] - readGDAL(paste(map.list[i], .asc,
sep=))$band1
}
proj4string(meuse.grid) - CRS(+proj=longlat +ellps=WGS84)
I do not think that it will be easy to generate kriging with any global
dataset that is 5 km resolution. You could try using the kriging in
SAGA
(it can handle maps up to 2 GB in size; it is exp(2) times faster than
R).
Here is an example:
http://spatial-analyst.net/wiki/index.php?title=Interpolation_of_ISRIC-
WISE_international_soil_profile_data
T. Hengl
PS: Which variable are you interpolating?
Hello All,
My questions are regarding the best way to obtain DEMs for the entire
globe. I am trying to krige a dataset of global scope. The dataset
has
points at all latitude/longitude intersects and so there are 64,800
sample
points. The problem that I am having is that the DEMs I have been
using
are
too large (memory) to be used as R variables. I also have had trouble
finding a way to download DEMs for the entire globe.
1. Given the size of my sample dataset: What is the maximum
resolution of
DEM that I should be using? Do I only need values at prediction
points or
will higher res yield higher accuracy?
2. Is there a way to programmatically access DEMs using GDAL or R in
general? ie. internally or through integration with McIDAS, NetCDF,
etc.
3. Do I need to break the DEMs and sample dataset into smaller tiles
to
avoid the memory issues? If so how do I do that in R / GDAL?
Thanks,
Abe
[[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
Re: [R-sig-Geo] Importing ESRI Shapefile into RSAGA
Jonas, SAGA (GUI) can read and write shape files. This is the native vector format it uses. Maybe this is confusing, but RSAGA is NOT a GIS package under R, but only a link to send things from R to SAGA and back. Here are some examples: http://spatial-analyst.net/wiki/index.php?title=Export_maps_to_GE#Polygon_maps HTH, T. Hengl http://home.medewerker.uva.nl/t.hengl/ -Original Message- From: r-sig-geo-boun...@stat.math.ethz.ch [mailto:r-sig-geo-boun...@stat.math.ethz.ch] On Behalf Of von Rütte Jonas Sent: Friday, October 02, 2009 11:31 AM To: r-sig-geo@stat.math.ethz.ch Subject: [R-sig-Geo] Importing ESRI Shapefile into RSAGA Hello! I have a more general problem Is there a way to import ESRI shapfiles into RSAGA? And if yes how can I then import it to SAGA GIS since it seems to not support the ESRI format. Thanks in advance Jonas Soil and Terrestrial Environmental Physics ETH Zurich [[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
Re: [R-sig-Geo] Determining a coordinate given another coordinate, a great circle distance, and an angle
Hi Harry, The key issue is that you need to select the (Euclidean) coordinate system. For example, to represent the whole world, there are many possiblities (http://www.radicalcartography.net/?projectionref). If you are working with the great circle distance, make sure that the geographic coordinate system uses a sphere. Say if you want to use some equidistant coordinate system with center at 0,0 longlat (Lambert Azimuthal Equal-Area), you could then get the coordinates by using: P1 - data.frame(E=35, N=25) lambda - 60 # azimuth! dist12 - 30 coordinates(P1) - ~E+N proj4string(P1) - CRS(+proj=longlat +datum=WGS84) P1.xy - spTransform(P1, CRS(+proj=laea +lat_0 +lon_0 +x_0 +y_0)) P2.xy - data.frame(e=p1...@coords[1,1]+dist12*sin(lambda), n=p1...@coords[1,2]+dist12*cos(lambda)) coordinates(P2.xy) - ~E+N proj4string(P2.xy) - p1...@proj4string P2 - spTransform(P2.xy, CRS(+proj=longlat +datum=WGS84)) P12 - rbind(P1, P2) P12 SpatialPoints: E N [1,] 35.0 25.000 [2,] 33.55488 22.554 Coordinate Reference System (CRS) arguments: +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 load(file=url(http://spatial-analyst.net/DATA/worldborders.RData;)) spplot(worldborders[NAME], colorkey=FALSE, col.regions=rep(grey(0.5), length(levels(worldborders$NAME))), sp.layout=list(sp.points, P12, pch=+, col=red, cex=3)) But I would consider using the UTM or which ever coordinate system is of your interest; then you also need to convert the GCD to distance in the local system. For example, to represent the whole of USA (contiguous 48-state area) I typically use: +proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=23 +lon_0=-96 +x_0=0 +y_0=0 +ellps=GRS80 +datum=NAD83 +units=m +no_defs See also: http://spatial-analyst.net/DATA/usgrids5km.zip cheers, T. Hengl http://spatial-accuracy.org/FromGEOSTAT2009 PS: I never got any photos from your trip to the Pyramids! -Original Message- From: r-sig-geo-boun...@stat.math.ethz.ch [mailto:r-sig-geo-boun...@stat.math.ethz.ch] On Behalf Of Harry Kim Sent: Tuesday, October 06, 2009 7:23 PM To: r-sig-geo@stat.math.ethz.ch Subject: [R-sig-Geo] Determining a coordinate given another coordinate,a great circle distance, and an angle Dear R-sig-geoers, I am stuck with a problem and i was hoping if you could help me. Suppose i am given a coordinate in lat/long (say 25 lat 35 long) and a great circle distance (say 10km). Also assume that i am given an angle with respect to latitude (an angle from the x-axis say 30 degrees). What would be the best way to determine a coordinate 10kms away from 25 lat 35 long with a given angle (30 degrees)? I think in euclidean space, I can use tangent of the given angle and combine with the Pythagorean theorem to find the target coordinate. Thank you very much in advance, Harry ___ 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] Including SAGA grid as GDAL-supported format
-Original Message- From: Roger Bivand [mailto:roger.biv...@nhh.no] Sent: Friday, October 16, 2009 11:06 AM To: Tomislav Hengl Cc: r-sig-geo@stat.math.ethz.ch Subject: Re: [R-sig-Geo] Including SAGA grid as GDAL-supported format On Fri, 16 Oct 2009, Tomislav Hengl wrote: Dear r-sig-geo, I would like to initiate the processes of registering the SAGA grid format under GDAL (the SAGA 2.0.4 source code is available for download here http://sourceforge.net/projects/saga-gis/files/SAGA%20- %202.0/SAGA%202.0.4/saga_2.0.4_src.zip/download). Do I have to follow some formal procedure, or do I have to prepare the GDAL driver myself (as explained at http://www.gdal.org/gdal_drivertut.html)? Any help/suggestions are welcome (apparently it should not be too complicated). Tom, Writing a read-only driver would be the first step. Once that works, it could be extended to copy creation, and finally to direct creation if justified. If SAGA I/O has a separate I/O library, it could be shipped with the driver as in the pcraster driver, but I don't see one in the source tree. Writing a driver for SAGA rasters would permit lots of other software to interact with SAGA, not just in R, so would be more robust. You'd develop on a local copy of the GDAL source, and once done submit a file archive of the directory to Frank Warmerdam or other GDAL people. Maybe ask on the GDAL developers' list if they have seen questions about SAGA rasters before. Are SAGA rasters similar to any existing supported formats (Erdas?, IDRISI?, ILWIS?) - could an existing driver be copied and modified? Roger, I've got an reply from FW who suggested that I should used the existing raw file + ascii header drivers as a template: http://svn.osgeo.org/gdal/trunk/gdal/frmts/raw/genbindataset.cpp We could try adopting this to see it will work with SAGA grids also. Apparently Volker Wichmann started implementing the SAGA driver himself (this summer), but then never finished (he might send me his code). Thanks for your suggestions and clarifications. Tom Roger ___ 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] inquiry on suitable packages
Hakim, There are plenty of (SDM) algorithms that generate species distribution maps given occurrence (or density) records and environmental maps. For a systematic review, see e.g. Tsoar et al. 2007 [http://dx.doi.org/10./j.1472-4642.2007.00346.x]. In R, you can try using the adehabitat package that can be used to generate Habitat Suitability Index maps using methods such as madifa and/or enfa (read more about ecology-analysis packages via the Environmetrics views). Here are some simple examples: http://spatial-analyst.net/scripts/bei.R I myselft have recently been testing MaxEnt [http://www.cs.princeton.edu/~schapire/maxent/] and discovered that it has several advantages over competition (GARP, ENFA): - it accepts both categorical and continuous predictors; - it allows cross-validation; - it is fast(est); In short, MaxEnt should be high on your list because it is also the most used SDM by biologists at the moment. There is no package for R to use MaxEnt, but there is probably no need for this neither. The whole algorithm comes in a single .jar file that can be run via batch commands (i.e. from R). Here are some examples using the global environmental layers: http://spatial-analyst.net/scripts/worldmaps.R Note also that there is also a R mailing list for ecology (R-sig-ecology) where you will probably be able to get much more feedback. HTH, T. Hengl http://home.medewerker.uva.nl/t.hengl/ Quoting Hakim Abdi hakim.a...@uni-muenster.de: Dear All: I have a set of nine independent variables (extracted Landsat TM band DNs, land surface temperature, and NDVI) and a set of eight response variables (presence/absence breeding season survey of eight bird species). I'm doing a small exercise to try to find a model that best fits the relationship between the independent variables and species occurrence and then apply that model to another set of independent variables from another time period (and even from another area) and see what the distribution might be at that time (or place). I am currently working with the packages BIOMOD Gstat and spatstat. I'm wondering if there are any others I'm not aware of that could further help me? Thanks in advance for your time and assistance. Hakim -- Mit freundlichen Grüßen/Regards __ Abdulhakim M. Abdi Erasmus Mundus Master's Programme in Geospatial Technologies Institut für Geoinformatik Universität Münster 32U 5755408.16mN, 404463.35mE ifgi.uni-muenster.de Instituto Superior de Estatística e Gestão de Informação Universidade Nova de Lisboa 29S 4287095.5mN, 486099.7mE www.isegi.unl.pt www.geospatialtechnologist.com www.hakimabdi.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
Re: [R-sig-Geo] Spatial analysis question
Marcelo,
Surprisingly, I could not find any function in the spatstat package (or splancs
package) that
specifically derives cross-correlations between multiple point processes:
data(lansing)
plot(split(lansing)) # distribution of occurrence records for five+1 species;
plot(density(split(lansing)), ribbon = FALSE)
# fit stationary marked Poisson process with different intensity for each
species:
lansing.ppm - ppm(lansing, ~marks, Poisson())
summary(lansing.ppm)
...but this does not say anything about which species are most correlated (and
which are negatively
correlated). See also Mark correlation function in PART V. MARKED POINT
PATTERNS:
Baddeley, A., 2008. Analysing spatial point patterns in R. CSIRO, Canberra,
Australia.
http://www.csiro.au/files/files/pn0y.pdf
I guess that there is no reason NOT to do what you suggest:
dens.lansing - density(split(lansing))
dens.lansing.sp - as(dens.lansing[[1]], SpatialGridDataFrame)
names(dens.lansing.sp)[1] - names(dens.lansing)[1]
for(i in 2:length(dens.lansing)) {
dens.lansing...@data[names(dens.lansing)[i]] - as(dens.lansing[[i]],
SpatialGridDataFrame)$v
}
spplot(dens.lansing.sp, col.regions=grey(rev(seq(0,1,0.025
round(cor(log1p(dens.lansing...@data[names(dens.lansing)]),
use=complete.obs), 2)
blackoak hickory maple misc redoak whiteoak
blackoak 1.000.55 -0.73 -0.64 -0.51 0.23
hickory 0.551.00 -0.84 -0.63 -0.52-0.27
maple -0.73 -0.84 1.00 0.75 0.50-0.09
misc-0.64 -0.63 0.75 1.00 0.70 0.09
redoak -0.51 -0.52 0.50 0.70 1.00 0.25
whiteoak 0.23 -0.27 -0.09 0.09 0.25 1.00
# PCA:
sp.formula - as.formula(paste(~, paste(log1p(, names(dens.lansing), ),
collapse=+),
sep=))
PCA.sp - prcomp(sp.formula, scale=TRUE, dens.lansing...@data)
biplot(PCA.sp, arrow.len=0.1, xlabs=rep(., length(PCA.sp$x[,1])), main=PCA
biplot,
ylabs=names(dens.lansing))
which clearly shows that the most positively correlated species are hickory
and blackoak, while
the most 'competing' species are maple/redoak and hickory.
HTH
T. Hengl
http://home.medewerker.uva.nl/t.hengl/
-Original Message-
From: r-sig-geo-boun...@stat.math.ethz.ch
[mailto:r-sig-geo-boun...@stat.math.ethz.ch] On Behalf
Of Marcelo Tognelli
Sent: Friday, October 30, 2009 7:50 PM
To: r-sig-geo@stat.math.ethz.ch
Subject: [R-sig-Geo] Spatial analysis question
Dear List,
I have probability maps of the distribution of 4 species of venomous snakes
(raster files output from species distribution modeling software) and point
locality data with information on snake bite events (most of them without
the id of the species involved in the accident). I would like to run an
analysis to see what species correlates best with snake bite events. My idea
is to generate a kernel density raster from the point event data and then do
some kind of spatial correlation against the species distribution maps.
I would greatly appreciate any suggestions on the type of analysis that I
can perform with these data and on the software and/or R package to run it.
Thanks in advance,
Marcelo
[[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
Re: [R-sig-Geo] Multiple predictors for external drift kriging
Roger Bivand wrote: On Tue, 10 Nov 2009, ageel bushara wrote: Dear list members, I'm doing spatial interpolation of soil moisture using kriging with external drift following the example given by Tom Hengl (http://spatial-analyst.net/wiki/index.php?title=Regression-kriging_guide). When I used one predictor (in my case gradient or cos_aspect, see below) it works fine for me, but when I used multiple predictors (gradient +cos_aspect), still I'm able to have the the experimental variogram as well as the fitted theoretical variogram, however, I'm not able to have the kriged soil moisture. here is the code # import and define the input variables library(sp) library(lattice) trellis.par.set(sp.theme()) # plots the final predictions using blue-pink-yellow legend water - read.csv(water_content.txt) # table with x,y coordinates,and z is? soil moisture str(water) coordinates(water)=~x+y?? # this makes depth a SpatialPointsDataFrame str(water) cos_aspect = read.asciigrid(cos_aspect.asc)? # reads ArcInfo Ascii raster map str(cos_aspect) spplot(cos_aspect, scales=list(draw=T), sp.layout=list(sp.points, water, pch=+)) gradient= read.asciigrid(gradient.asc) str(gradient) spplot(gradient, scales=list(draw=T), sp.layout=list(sp.points, water, pch=+)) #Plot the xy graph target versus predictor: cos_aspect.ov = overlay(cos_aspect, water)? # create grid-points overlay str(cos_aspect...@data) water$cos_aspect.asc = cos_aspect.ov$cos_aspect.asc? gradient.ov = overlay(gradient, water) water$gradient.asc= gradient.ov$gradient.asc lm.depth - lm(z~ gradient.asc+cos_aspect.asc, as.data.frame(water)) summary(lm.depth) plot(z~ gradient.asc+cos_aspect.asc, as.data.frame(water)) abline(lm(z~ gradient.asc, as.data.frame(water))) abline(lm(z~ cos_aspect.asc, as.data.frame(water))) #Fit the variogram model of the residuals: library(gstat) null.vgm - vgm(var(water$z), Sph, sqrt(areaSpatialGrid(cos_aspect))/4, nugget=0) #initial parameters vgm_depth_r - fit.variogram(variogram(z~ cos_aspect.asc+ gradient.asc, water), model=null.vgm) plot(variogram(z~ gradient.asc+cos_aspect.asc,water), vgm_depth_r, main=fitted by gstat) # It works fine till here # Run uk in gstat: depth_uk - krige(z~ cos_aspect.asc+gradient.asc, locations=water, newdata=cos_aspect, model=vgm_depth_r) the problem it seems mainly in newdata, here is the error given by R: Error in eval(expr, envir, enclos) : object 'gradient.asc' not found. if i assigned the? newdata= gradient then the R error message is : Error in eval(expr, envir, enclos) : object 'cos_aspect.asc' not found You probably need to think through what you are doing. You should use a Spatial*DataFrame for newdata that includes cos_aspect.asc and gradient.asc among the values returned by the names() method. Since each of your two SpatialGridDataFrame objects have single variables, you should review what a SpatialGridDataFrame is - a data.frame with a description of the GridTopology. Assuming that your two objects share their GridTopology values, maybe cbind() them? Then check that the names() output includes the names of variables on the RHS in the formula. Don't think of objects as grids, think of them as data.frames, then you see what is going on. Try with lm() and predict() first, then with gstat() and predict() - the newdata= argument works in exactly the same way. Hope this helps, Roger Hi Ageel, Roger is right. You should read all rasters to a single (multilayer) grid e.g.: grids - readGDAL(cos_aspect.asc) names(grids) - cos_aspect.asc grids$gradient.asc - readGDAL(gradient.asc)$band1 ... Then overlay, filter the NA's, and then you can fit a variogram and run predictions e.g.: depth_uk - krige(z~ cos_aspect.asc+gradient.asc, locations=water, newdata=grids, model=vgm_depth_r) ... Here are some more examples: http://geomorphometry.org/content/some-examples-rsaga http://spatial-analyst.net/scripts/meuse.R all the best, T. Hengl http://home.medewerker.uva.nl/t.hengl/ as I understand that newdata is the grid where it does the interpolations How can I obtained the kriged soil moisture?, what is wrong here?, Thanks, Ageel [[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] optimizeNetwork function in intamapInteractive package (how to use auxiliary predictors)
Dear Edzer, Jon, Stephanie, Olivier et al., First, let me congratulate you on your new package intamap / intamapInteractive that I feel will significantly advance practice of geostatistical mapping. Really a great contribution (I knew about SSA already in 2000, but was always lacking an operational tool). My interest at the moment are the sampling optimization algorithms and their use in combination with universal kriging / regression-kriging. I had no problems to run spatially simulated annealing with your package (OK model), but could not grasp how to do the same thing using a universal kriging model (auxiliary predictors). Here is the code: library(gstat) library(RSAGA) library(maptools) library(intamapInteractive) # load data: data(meuse) coordinates(meuse) - ~x+y data(meuse.grid) coordinates(meuse.grid) - ~x+y gridded(meuse.grid) - TRUE # estimate variograms (OK/UK): vt.fit - fit.variogram(variogram(log1p(zinc)~1, meuse), vgm(1, Exp, 300, 1)) vr.fit - fit.variogram(variogram(log1p(zinc)~sqrt(dist), meuse), vgm(1, Exp, 300, 1)) # study area of interest: write.asciigrid(meuse.grid[mask], meuse_mask.asc) rsaga.esri.to.sgrd(in.grids=meuse_mask.asc, out.sgrd=meuse_mask.sgrd, in.path=getwd()) # raster to polygon conversion; rsaga.geoprocessor(lib=shapes_grid, module=6, param=list(GRID=meuse_mask.sgrd, SHAPES=meuse_mask.shp, CLASS_ALL=1)) mask - readShapePoly(meuse_mask.shp, proj4string=CRS(as.character(NA)), force_ring=T) # prepare observations / predGrid objects: observations - meuse[zinc] attr(observati...@coords, dimnames)[[2]] - c(x, y) predGrid - data.frame(meuse.grid[dist]) predGrid$dist - NULL # predGrid has to be SpatialPoints type? coordinates(predGrid) - ~x+y windows() # add 20 more points assuming OK model (SSA method): optim.ok - optimizeNetwork(observations, predGrid, candidates=mask, method=ssa, action=add, nDiff=20, model=vt.fit, criterion=MUKV, plot=FALSE, nr_iterations=50, nmax=40) This works fine. But how do I now optimize samples using universal kriging model (with auxiliary predictors e.g. dist)? I simply could not figure out from your examples. I also tried looking at the Package source R functions, but could not get it running. It should be something like this: observations - pts.ov[c(zinc, dist)] predGrid - data.frame(meuse.grid[dist]) coordinates(predGrid) - ~x+y optim.rk - optimizeNetwork(observations, predGrid, candidates=mask, method=ssa, action=add, nDiff=20, model=vr.fit, criterion=MUKV, plot=FALSE, nr_iterations=50, formula=log1p(zinc)~sqrt(dist), nmax=40) thanks! Tom Hengl http://home.medewerker.uva.nl/t.hengl/ ___ 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] Unconditional simulation
Hi Nick, Very interesting problem. At first thought, I imagined that you just want to simulate noise ;) In the geoR package (http://leg.ufpr.br/geoR/) there is a function to simulate Gaussian Random Fields (uses actually RandomFields package) using various models e.g.: library(geoR) - Analysis of geostatistical data For an Introduction to geoR go to http://www.leg.ufpr.br/geoR geoR version 1.6-27 (built on 2009-10-15) is now loaded - s - grf(100, grid=reg, cov.pars=c(2, 0.2), cov.model=mat, kappa=1.5) grf: generating grid 10 * 10 with 100 points grf: process with 1 covariance structure(s) grf: nugget effect is: tausq= 0 grf: covariance model 1 is: matern(sigmasq=2, phi=0.2, kappa = 1.5) grf: decomposition algorithm used is: cholesky grf: End of simulation procedure. Number of realizations: 1 image(s, col=gray(seq(1, 0.2, l=21))) hist(s$data) # normal distribution You can also simulate a regular point sample with the same spatial structure on top of that using either Poisson, Bernoulli or binomial models. For example, to simulate a Poisson model, you could use: # define your own model, e.g. poisson: lambda - exp(0.5 +s$data) y - rpois(length(s$data), lambda=lambda) points(y) text(s$coords, label=y, pos=3, offset=0.3) hist(y) For a uniform model, I would then use the Empirical Cumulative Distribution Function (ECDF): # uniform distribution: y.cdf - ecdf(s$data) y - y.cdf(s$data) image(s, col=gray(seq(1, 0.2, l=21))) points(y) text(s$coords, label=y, pos=3, offset=0.3) hist(y) This would then have the same spatial auto-correlation structure and 'perfectly' uniform distribution (I might also be wrong - I do not like that a simulated variable has a perfect histogram). I am sure that other mathematicians have maybe better ideas. HTH T. Hengl http://home.medewerker.uva.nl/t.hengl/ -Original Message- From: r-sig-geo-boun...@stat.math.ethz.ch [mailto:r-sig-geo-boun...@stat.math.ethz.ch] On Behalf Of Nick Hamm Sent: Tuesday, November 17, 2009 10:06 AM To: r-sig-geo@stat.math.ethz.ch; ai-geost...@jrc.it Subject: [R-sig-Geo] Unconditional simulation Dear all I want to simulate a spatially-correlated random field which follows a uniform rather than than Gaussian distribution. Does anybody know a straight-forward way to do this? Nick ___ 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] 'LDLfactor' error in 'krige' function
Hi Paul, Edzer, I understand why the singular matrix problem happens and I know that there is not really a mathematical solution around it: x - matrix(runif(100), nrow=10) x.i - solve(x) str(x.i) num [1:10, 1:10] 0.8191 -1.0293 0.0826 1.068 -0.2106 ... # add a 'singular' column x[,1:10] - rep(1, 10) x.i - solve(x) Error in solve.default(x) : Lapack routine dgesv: system is exactly singular However, I would very much support if you would integrate an if loop to check if it will happen, and then mask the prediction location using NA. I mean, at the moment every time we want to run predictions, even if only at a single prediction location the model fails, we are not able to generate any output (this is sometimes really frustrating). Hence, I would support that you allow us to run predictions for all locations first, and then let the users judge if one should increase the search radius, remove some too-uniform predictors etc etc. I hope you agree with me, T. Hengl http://home.medewerker.uva.nl/t.hengl/ -Original Message- From: r-sig-geo-boun...@stat.math.ethz.ch [mailto:r-sig-geo-boun...@stat.math.ethz.ch] On Behalf Of Paul Hiemstra Sent: Tuesday, November 17, 2009 9:53 AM To: Edzer Pebesma Cc: r-sig-geo Subject: Re: [R-sig-Geo] 'LDLfactor' error in 'krige' function Edzer Pebesma wrote: My guess is that from constant data, the (co)variance is constant and zero, so the covariance matrix cannot be decomposed (hence: LDLfactor errors). Is this a case that autokrige should catch? I added a check in autoKrige. The output for the example below is now: kriging_result = autoKrige(zinc~1, meuse) Error in autoKrige(zinc ~ 1, meuse) : All data in attribute 'zinc' is identical and equal to 0 Can not interpolate this data A new version of automap with this feature has been uploaded to CRAN. cheers, Paul -- Edzer Mauricio Zambrano wrote: Dear List, During some OK interpolations of daily precipitation, with the 'automap' library, I got the following error: [using ordinary kriging] chfactor.c, line 130: singular matrix in function LDLfactor() Error en predict.gstat(g, newdata = newdata, block = block, nsim = nsim, : LDLfactor The code I'm using works fine for other days, so I assume that the distance between the measurement points is no the problem. When looking at the data that rose the error, I realized that all the measured values were equal to zero (I can not skip those days in which all the measured points have the same value in advance, because the measured value in those points change with time). According to a traceback that is given below, it seems that the cause is in the 'predict.gstat' function of the 'gstat' package. The same error can be risen with: # Data preparation data(meuse) coordinates(meuse) =~ x+y data(meuse.grid) gridded(meuse.grid) =~ x+y meuse$zinc - meuse$zinc*0 meuse$zinc # Ordinary kriging, no new_data object kriging_result = autoKrige(zinc~1, meuse) traceback() 6: .Call(gstat_predict, as.integer(nrow(as.matrix(new.X))), as.double(as.vector(raw$locations)), as.vector(new.X), as.integer(block.cols), as.vector(block), as.vector(bl_weights), as.integer(nsim), as.integer(BLUE)) 5: predict.gstat(g, newdata = newdata, block = block, nsim = nsim, indicators = indicators, na.action = na.action, debug.level = debug.level) 4: .local(formula, locations, ...) 3: krige(formula, input_data, new_data, variogram_object$var_model, block = block, ...) 2: krige(formula, input_data, new_data, variogram_object$var_model, block = block, ...) 1: autoKrige(zinc ~ 1, meuse) I would really appreciate any hint about the possible reason of this error and if it could be overcome in some way. Thanks in advance for any help. Mauricio -- 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 274 3113 Mon-Tue Phone: +3130 253 5773 Wed-Fri http://intamap.geo.uu.nl/~paul ___ 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] Spatial interpolation of river network observations
Hi Ultrich, Facundo Muñoz apparently made a GRASS function to derive distances along streams: https://stat.ethz.ch/pipermail/r-sig-geo/2009-November/006851.html 17 observations only? That is really tight for any geostatistical analysis (on top, you want to do 3 dimensions!). I would instead consider fitting a smooth surface (splines) / simulating the flow processes. all the best, T. Hengl http://home.medewerker.uva.nl/t.hengl/ -Original Message- From: r-sig-geo-boun...@stat.math.ethz.ch [mailto:r-sig-geo-boun...@stat.math.ethz.ch] On Behalf Of Ulrich Leopold Sent: Tuesday, November 17, 2009 10:53 AM To: R-sig-geo list; grass-u...@lists.osgeo.org Subject: [R-sig-Geo] Spatial interpolation of river network observations Dear all, I would like to interpolate 17 pollution observations in a storage lake in 3 dimensions (x,y,z). As I understand variogram analysis and kriging are not straightforward as we are dealing with non-euclidean (hydrologic) distances and down-stream direction. Could someone point me to some algorithms which can roughly estimate the 3d pollution body accounting for hydrologic distances? Attached is a file which shows the sampling locations and the storage lake. Thanks very much. Best regards, Ulrich -- __ Ulrich Leopold Resource Centre for Environmental Technologies, Public Research Centre Henri Tudor, Technoport Schlassgoart, 66 rue de Luxembourg, P.O. BOX 144, L-4002 Esch-sur-Alzette, Luxembourg tel: +352 42 5591 618 fax: +352 42 5591 555 mobile: +352 691 304813 http://www.crte.lu , http://www.tudor.lu ___ 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] 'LDLfactor' error in 'krige' function
Hi Edzer, Thanks for the info. I was not aware that I can simply set cn_max and the predictions will not break (and I still do read the gstat manual). I guess that this is then the solution to our problem. For me it enough to generate a map with NA's, then zoom into map to see where the problematic areas/points are, then either filter the maps, consider using different search radius or threshold. T. Hengl http://home.medewerker.uva.nl/t.hengl/ -Original Message- From: Edzer Pebesma [mailto:edzer.pebe...@uni-muenster.de] Sent: Tuesday, November 17, 2009 1:01 PM To: Tomislav Hengl Cc: 'r-sig-geo' Subject: Re: [R-sig-Geo] 'LDLfactor' error in 'krige' function Tom, you can already do this: library(gstat) data(meuse) coordinates(meuse)=~x+y data(meuse.grid) gridded(meuse.grid)=~x+y meuse=meuse[c(1,1:155),] # replicate first observation pr1 = predict(zinc~1,meuse,meuse.grid,vgm(1, Exp, 300)) # will break # the following will generate NA's for cells where the estimated condition number of the covariance matrix exceeds 1e10: pr1 = krige(zinc~1,meuse,meuse.grid,vgm(1, Exp, 300),nmax=30, set=list(cn_max=1e10)) # generate a full interpolation as comparison: pr1$idw = idw(zinc~1,meuse,meuse.grid)[[1]] # show side by side: spplot(pr1, c(var1.pred, idw), col.regions=bpy.colors()) ... missing values are generated for those cells that result in a singular covariance matrix, given their local neighbourhood of 30 nearest observations. cn_max referers to the maximum allowed condition number, see http://en.wikipedia.org/wiki/Condition_number Two issues are (i) what singularity means when we use floating point representations for real numbers, and (ii) that condition numbers of matrices are expensive to evaluate, and an estimate based on the LU decomposition is used. I threshold to 1e10 here, but that is purely for illustrational purposes. Note, for you of interest, that IIRC this thresholding is also done for singularity of the X matrix, holding the predictors. All this information didn't make it to the help pages of the krige function. This help page would cover tens of pages otherwise. For full documentation, still the old, gstat stand-alone manual on gstat.org is needed. I agree that this is not optimal. With best wishes, -- Edzer Tomislav Hengl wrote: Hi Paul, Edzer, I understand why the singular matrix problem happens and I know that there is not really a mathematical solution around it: x - matrix(runif(100), nrow=10) x.i - solve(x) str(x.i) num [1:10, 1:10] 0.8191 -1.0293 0.0826 1.068 -0.2106 ... # add a 'singular' column x[,1:10] - rep(1, 10) x.i - solve(x) Error in solve.default(x) : Lapack routine dgesv: system is exactly singular However, I would very much support if you would integrate an if loop to check if it will happen, and then mask the prediction location using NA. I mean, at the moment every time we want to run predictions, even if only at a single prediction location the model fails, we are not able to generate any output (this is sometimes really frustrating). Hence, I would support that you allow us to run predictions for all locations first, and then let the users judge if one should increase the search radius, remove some too-uniform predictors etc etc. I hope you agree with me, T. Hengl http://home.medewerker.uva.nl/t.hengl/ -Original Message- From: r-sig-geo-boun...@stat.math.ethz.ch [mailto:r-sig-geo-boun...@stat.math.ethz.ch] On Behalf Of Paul Hiemstra Sent: Tuesday, November 17, 2009 9:53 AM To: Edzer Pebesma Cc: r-sig-geo Subject: Re: [R-sig-Geo] 'LDLfactor' error in 'krige' function Edzer Pebesma wrote: My guess is that from constant data, the (co)variance is constant and zero, so the covariance matrix cannot be decomposed (hence: LDLfactor errors). Is this a case that autokrige should catch? I added a check in autoKrige. The output for the example below is now: kriging_result = autoKrige(zinc~1, meuse) Error in autoKrige(zinc ~ 1, meuse) : All data in attribute 'zinc' is identical and equal to 0 Can not interpolate this data A new version of automap with this feature has been uploaded to CRAN. cheers, Paul -- Edzer Mauricio Zambrano wrote: Dear List, During some OK interpolations of daily precipitation, with the 'automap' library, I got the following error: [using ordinary kriging] chfactor.c, line 130: singular matrix in function LDLfactor() Error en predict.gstat(g, newdata = newdata, block = block, nsim = nsim, : LDLfactor The code I'm using works fine for other days, so I assume that the distance between the measurement points is no the problem. When looking at the data that rose the error, I realized that all the measured values were equal to zero (I can not skip those days
Re: [R-sig-Geo] Filling in holes in DTM
Hi Cara, There is a very efficient function in SAGA called Close Gaps that does exactly that. What makes it especially efficient is that it allows you to set a mask map. See: rsaga.get.usage(grid_tools, 7) SAGA CMD 2.0.3 library path: C:/Progra~1/saga_vc/modules library name: grid_tools module name : Close Gaps Usage: 7 -INPUT str [-MASK str] [-RESULT str] [-THRESHOLD str] -INPUT:str Grid Grid (input) -MASK:str Mask Grid (optional input) -RESULT:str Changed Grid Grid (optional output) -THRESHOLD:str Tension Threshold BR, T. Hengl http://home.medewerker.uva.nl/t.hengl/ -Original Message- From: r-sig-geo-boun...@stat.math.ethz.ch [mailto:r-sig-geo-boun...@stat.math.ethz.ch] On Behalf Of Tobin Cara Sent: Thursday, November 19, 2009 3:16 PM To: r-sig-geo@stat.math.ethz.ch Subject: [R-sig-Geo] Filling in holes in DTM Hello, Would anyone know a good way to fill in holes within a DTM? There is no data inside these holes and it is affecting my calculations, so I prefer for an interpolation of the nearest neighbors to fill in a value or something similar. Thank you, Cara [[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
Re: [R-sig-Geo] A new ASTER Global DEM data set
Hi all, FYI: I've run a small comparison between ASTER DEM and LIDAR DEMs (say 'true' topography): http://geomorphometry.org/content/gdem-quick-assessment The 4 case studies can be downloaded from here: http://geomorphometry.org/system/files/GDEM_assessment.zip I'm not too happy with what I've got - GDEM shows very strange patterns in area of low relief and I definitively think that the 30 m resolution is overoptimistic; it should be degraded to 60-90 m (this would enhance the data sharing and save them a lot of trouble). The positional accuracy of GDEM is on the other hand pretty good. GDEM is a frankestein - once you zoom in and add some shading, you can see the stitches, so have this in mind (unlike SRTM DEM which is a complete and a consistent project; but also noisy in many areas). T. Hengl http://home.medewerker.uva.nl/t.hengl/ -Original Message- From: r-sig-geo-boun...@stat.math.ethz.ch [mailto:r-sig-geo-boun...@stat.math.ethz.ch] On Behalf Of Ashton Shortridge Sent: Friday, November 20, 2009 3:15 PM To: r-sig-geo@stat.math.ethz.ch Cc: Yong Li Subject: Re: [R-sig-Geo] A new ASTER Global DEM data set On Friday 20 November 2009 07:09:35 Barry Rowlingson wrote: 2009/11/20 Yong Li yong...@unimelb.edu.au: Dear all folks, I was informed in the 6th Digital Earth Conference that there is a better place to acquire high resolution of global DEM developed by ASTER, called GDEM with 30 m resolution, and fantastically free of charge. I tried some here (http://www.gdem.aster.ersdac.or.jp/) and it is really better than SRTM if you are outside USA. Hope you will enjoy the free meal. I just had a look at the tiles near me, and there's quite a bit of noise and obvious artefacts - unless there really is a 2000m tower about 180m across that I've not noticed in the middle of the countryside! I've not compared with SRTM yet... Barry I looked at some here around the North American Great Lakes, and I have to say I would be very leery about using it in low relief areas. Visually it looked a lot nicer in the California coastal range, so perhaps less vegetation and higher relief is important for the sensor. That said, it's terrific that alternative global medium-resolution DEMs are becoming available. Also, I think this ASTER-derived product has captured higher latitude locations than SRTM (which gets to about 60N and 60S), so it may be not simply the best but the only choice for many regions. Ashton -- Ashton Shortridge Associate Professor ash...@msu.edu Dept of Geography http://www.msu.edu/~ashton 235 Geography Buildingph (517) 432-3561 Michigan State University fx (517) 432-1671 ___ 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] Indicator Kriging R/gstat
Dear Pete, Edzer, If this is of any help, few years ago we have played with using auxiliary maps to interpolate categorical variables (multi-indicators?). The results are reported in: Hengl T., Toomanian N., Reuter H.I., Malakouti M.J. 2007. Methods to interpolate soil categorical variables from profile observations: lessons from Iran. Geoderma, 140(4): 417-427. http://dx.doi.org/10.1016/j.geoderma.2007.04.022 Experiences were frustrating. Many (most of) classes come in isolated patches; hence getting something out from the variograms was difficult (close to pure nugget). In addition, if you go ahead and interpolate the 0/1 values (ignoring binomial properties), you will produce values 0 and 1 i.e. nonsense values. Then we tried with using memberships and things changed. First, memberships you can simply convert to logits, then interpolate and back-transform and then none of the values will be outside the 0-1 range. It was also much easier to fit the variograms. But we had a luxury that soils are fuzzy classes per definition and we had extra data to 'fuzzify' the A, B, C values to memberships (in fact we discovered that surveyors should have designated memberships from the beginning - this solves all problems of transition classes etc). Here is a similar type of approach (fit a GLM with a binomial model, interpolate the residuals using OK) with the meuse dataset: data(meuse) coordinates(meuse) - ~x+y data(meuse.grid) coordinates(meuse.grid) - ~x+y gridded(meuse.grid) - TRUE fullgrid(meuse.grid) - TRUE meuse.ov - overlay(meuse.grid, meuse) meuse...@data - cbind(meuse...@data, meuse[c(zinc, lime)]...@data) # fit a GLM: glm.lime - glm(lime~dist+ffreq, meuse.ov, family=binomial(link=logit)) step.lime - step(glm.lime) # check if the predictions are within 0-1 range: summary(round(step.lime$fitted.values, 2)) # convert to logits: logits.lime - log(step.lime$fitted.values/(1-step.lime$fitted.values)) # predict at all locations: p.glm - predict(glm.lime, newdata=meuse.grid, type=link, se.fit=T) # predictions as logits str(p.glm) # attach spatial reference: lime.glm - as(meuse.grid, SpatialPointsDataFrame) lime.glm$lime - p.glm$fit gridded(lime.glm) - TRUE lime.glm - as(lime.glm, SpatialGridDataFrame) # fit a variogram for residuals: lime.ivgm - vgm(nugget=0, model=Exp, range=sqrt(diff(me...@bbox[x,])^2 + diff(me...@bbox[y,])^2)/4, psill=var(residuals(step.lime))) lime.rvgm - fit.variogram(variogram(residuals(step.lime)~1, meuse.ov), model=lime.ivgm) # does not work - singular model; fix by hand: lime.rvgm - vgm(nugget=0.6, Exp, psill=0.2, range=500) # interpolate residuals (logits): lime.rk - krige(residuals(step.lime)~1, meuse.ov, meuse.grid, lime.rvgm) # add predicted trend and residuals: lime.rk$var1.rk - lime.glm$lime + lime.rk$var1.pred # back-transform logits to the original response scale (0,1): lime.rk$var1.rko - exp(lime.rk$var1.rk)/(1+exp(lime.rk$var1.rk)) spplot(lime.rk[var1.rko], scales=list(draw=T), at=seq(0.05,1,0.05), col.regions=grey(rev(seq(0,0.95,0.05))), main=Liming requirements, sp.layout=list(sp.points, col=black, meuse)) write.asciigrid(lime.rk[var1.rko], lime_rk.asc, na.value=-1) Another thing you should look at is the geoRglm package (it has a function binom.krige to work with binomial data), but as far as I remember - I could not get it running. Maybe you should dig into that and then contact the authors if you need more help. HTH, T. Hengl http://home.medewerker.uva.nl/t.hengl/ -Original Message- From: r-sig-geo-boun...@stat.math.ethz.ch [mailto:r-sig-geo-boun...@stat.math.ethz.ch] On Behalf Of Edzer Pebesma Sent: Sunday, November 29, 2009 9:10 PM To: Pete Larson Cc: r-sig-geo@stat.math.ethz.ch Subject: Re: [R-sig-Geo] Indicator Kriging R/gstat Pete, for the universal kriging case I don't believe that indicator theory allows what you want to do. See also https://stat.ethz.ch/pipermail/r-sig-geo/2009-May/005701.html If you disagree, could you please let me know why, or where you found indications that this can be done? For ordinary kriging, usually values outside [0,1] are set back to their nearest allowed value. The degree to which they will be outside [0,1] also depends on the variogram used -- Gaussian variograms being notorously suspect. -- Edzer Pete Larson wrote: Hello, I would like to do indicator kriging in R/gstat. I have dichotomous 0/1 data and have performed ordinary kriging and universal kriging, but get predctions that are far from 0 and 1. Am I doing something wrong? Here is the code I have been using: Ordinary Kriging with krige function # # ordinary kriging: x - krige(z~1, ~x+y, model = v.fit, data = estand, newd = grd) ### universal block kriging: uk - krige(z~x+y+DistHF+RiverDist+RiverDist2, ~x+y, model = v.fit, data = estand, newdata = grd) Any help would be appreciated. Pete
Re: [R-sig-Geo] Spatial R Courses
David Rossiter runs this excellent distance education course: http://www.itc.nl/Pub/Study/Courses/C10-AES-DE-02 *Next one starts in 6 weeks! We will run a 5-day trainig course in June (which is more affordable, but shorter): 27 June - 4 July 2010, GEOSTAT 2010 Summer School for PhD students, Plasencia, Spain (more info soon) HTH T. Hengl http://home.medewerker.uva.nl/t.hengl/ Michael Denslow wrote: Dear r-sig-geo, I apologize if this is too off topic. I am interested in taking a online or other distance learning course in spatial analysis/programming. An R course would be great, but I am also interested in other things such as GRASS, Python or open web applications. The only caveats are that it must be a graduate course in Geography (or closely related field). Thanks in advance for any ideas, Michael ___ R-sig-Geo mailing list R-sig-Geo@stat.math.ethz.ch https://stat.ethz.ch/mailman/listinfo/r-sig-geo
[R-sig-Geo] Making web-repository of gridded maps: NetCDF or WKT Raster?
Dear R-sig-geo, As a part of our project (EcoGRID.nl) we have prepared some 60 thematic grids that we use as auxiliary predictors for species distribution modeling. At this stage, we would like to put the gridded maps (50/100 m base resolution) into some efficient sharable geo-database. We will most probably put the data into the NetCDF format (http://www.unidata.ucar.edu/software/netcdf/) because it can handle any-dimensional array data, and because it has been in continuous development and widespread use for many years. NetCDF grids can be read relatively easy into R using the RNetCDF package (e.g. http://spatial-analyst.net/DATA/readNCDF.zip). Another alternative is to use PostGIS WKT Raster format (http://trac.osgeo.org/postgis/wiki/WKTRaster), but this seems to be still rather experimental (?). Once we put the grids into NetCDF format. We plan to install OPeNDAP (http://www.opendap.org) server on top to make the files accessible through the web; then Geoserver (http://www.geoserver.org) or UMN Mapserver (http://mapserver.org) to feed a WMS from NetCDF files (raster data). Finally, we plan to add a simple OpenLayers interface on top of that (Geoserver has it built in) to allow direct browsing of the data and metadata (e.g. such as this one: http://africamap.harvard.edu/) Just to be clear, we want to put the data on a server because we would like to run a number of operations directly on the server (via rgdal?): 1. Overlay some point data and get the values of grids (without a need to download the grids locally); 2. Subset/mask and resample gridded data of interest (for a given bounding box and proj4 string; again without a need to download the data locally); 3. Upscale the grids from 100 m to 250, 500 m and 1 km resolution (then download the upscaled grids). 4. Write/upload new grids to the database (e.g. using WebDAV). 5. Browse the grids (and metadata) via the OpenLayers. These are only our wishes of course. We do not know if all this is really possible with the current software. Any examples or comments/suggestions/experiences are welcome (before we start installing and testing the functionality). Thanx! Tomislav Hengl and Lourens Veen ___ 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] Making web-repository of gridded maps: NetCDF or WKT Raster?
Dear Edzer, Barry, Servet, and Miguel, Thank you for your tips and for sharing your experiences with us (R-sig-geo keeps on surprising me considering the amount of support I get :). I will dig into the links you have forwarded. Most probably, we will try to implement two/three paths and then compare them to see which is 'better' (the main issue is that we want to do majority of processing on the server): (A) The rasdaman approach (http://kahlua.eecs.jacobs-university.de/trac/rasdaman/wiki/Rasdaman_8.0) 1. Install PostgreSQL support on server and download and install rasdaman File Set; 2. Create a database and put the grids in PostgreSQL database (raslib applications); 3. Test some examples with overlaying, subsetting, resampling, upscaling and visualizing data (WMS, KML); (B) The OPeNDAP/NetCDF approach (http://www.unidata.ucar.edu/software/netcdf/software.html#GDAL) 1. Install OPeNDAP server and UMN Mapserver to serve data; 2. Put all grids into a single (CF)NetCDF file and upload (using WebDAV) to the server; 3. Test some examples (R + GDAL) with overlaying, subsetting, resampling, upscaling and visualizing the data (WMS, KML); I will keep you informed once we produce some live examples. We might first start testing a WCPS with the worldmaps repository (http://spatial-analyst.net/worldmaps/); then develop some examples (R code) and test flexibility/speed of both approaches. cheers, T. Hengl http://home.medewerker.uva.nl/t.hengl/ -Original Message- From: r-sig-geo-boun...@stat.math.ethz.ch [mailto:r-sig-geo-boun...@stat.math.ethz.ch] On Behalf Of Edzer J. Pebesma Sent: Sunday, December 13, 2009 9:05 PM To: Stefano Casalegno Cc: r-sig-geo Subject: Re: [R-sig-Geo] Making web-repository of gridded maps: NetCDF or WKT Raster? Dear Tom, as your grids are spatial grids, you may want to look into the web service standards of the OGC (www.opengeospatial.org) for them, WCS (for reading) and WCS-T (for writing) in particular. We're experimenting with this currently, using the recently open sourced rasdaman software (www.rasdaman.org), which uses a data base in the back end. I don't think OPeNDAP has special infrastructure for spatial grids; CFNetCDF does, but by choosing this you may limit integration with other software that does assume spatial data to be served over OGC standards. CFNetCDF is one of the possible formats to be served over WCS (and GeoTIFF etc) -- if I'm correctly informed. Bests, -- Edzer Stefano Casalegno wrote: Dear Tom, my best acknowledgment for the Ecogrid challenging initiative. I guess the main limitation will be time and this is related to the size of grids. What is the extent of the grids ? What tool you plan to use for point 1. Overlay some point data and get the values of grids ? On a more theoretical basis, distribution models should rely upon data that covers the whole extent of the distribution of the biota under examination Please consider forest data from http://forest.jrc.ec.europa.eu/climate-change as supplementary predictors variables within an European biodiversity and nature conservation application. Cheers Stefano Stefano Casalegno, Ph.D. via Greppi 16, 21021 Angera - Italy email: stef...@casalegno.net web: http://www.casalegno.net/stefano On Dec 13, 2009, at 6:29 PM, Tomislav Hengl wrote: Dear R-sig-geo, As a part of our project (EcoGRID.nl) we have prepared some 60 thematic grids that we use as auxiliary predictors for species distribution modeling. At this stage, we would like to put the gridded maps (50/100 m base resolution) into some efficient sharable geo-database. We will most probably put the data into the NetCDF format (http:// www.unidata.ucar.edu/software/netcdf/) because it can handle any- dimensional array data, and because it has been in continuous development and widespread use for many years. NetCDF grids can be read relatively easy into R using the RNetCDF package (e.g. http:// spatial-analyst.net/DATA/readNCDF.zip). Another alternative is to use PostGIS WKT Raster format (http://trac.osgeo.org/postgis/wiki/ WKTRaster), but this seems to be still rather experimental (?). Once we put the grids into NetCDF format. We plan to install OPeNDAP (http://www.opendap.org) server on top to make the files accessible through the web; then Geoserver (http:// www.geoserver.org) or UMN Mapserver (http://mapserver.org) to feed a WMS from NetCDF files (raster data). Finally, we plan to add a simple OpenLayers interface on top of that (Geoserver has it built in) to allow direct browsing of the data and metadata (e.g. such as this one: http://africamap.harvard.edu/) Just to be clear, we want to put the data on a server because we would like to run a number of operations directly on the server (via rgdal?): 1. Overlay
[R-sig-Geo] GEOSTAT 2010 Summer School, 27 June - 4 July 2010; University of Plasencia, Spain
GEOSTAT 2010 Summer School 27 June - 4 July 2010; University of Plasencia, Spain This is the first call for participation in the GEOSTAT 2010 Summer School. GEOSTAT focuses on important aspects of statistical analysis of spatial and spatio-temporal data using open source / free GIS tools: R, SAGA GIS, GRASS GIS, Quantum GIS, GDAL, Google Earth and similar. The course participants learn how to move data back and forth between different environments, how to produce scripts and automate analysis. We welcome also R beginners and users needing refresh courses in programming. This year, we would also like to introduce/promote topics such as: web-based computing, WPS client-server environments, 3D and 4D geostatistics, combining R+SAGA/GRASS. This is a 5-day course with two parallel sessions, which means that there will be total 7 full-day blocks (three days with parallel sessions) of lectures; the last day of the summer school participants can present their research problems and ask for feedback from the whole summer school. For more info see: http://geostat2010.info/Info - LECTURERS: Roger Bivand (nhh.no), Edzer Pebesma (uni-muenster.de), Gerard Heuvelink (wur.nl), Olaf Conrad (geowiss.uni-hamburg.de), Tomislav Hengl (uva.nl) and Victor Olaya Ferrero (unex.es) - DEADLINE TO REGISTER: 15th February 2010; online at: http://geostat2010.info/deadline_registration - SELECTION OF CANDIDATES: The participants will be selected using the following criteria: 1. your publication record / academic position; 2. previous involvement in the open source activities (R-sig-geo, conferences and workshops) 3. distance from the venue (more distant applications have advantage); 4. time of application (first-come-first-served); - FEES: To be determined; usually between 200-300 EUR (excluding accommodation costs). 5-6 double bed rooms (in-house accommodation) available at subsidized price for participants who reside in, and are nationals of an ODA country (http://www.oecd.org/dataoecd/32/40/43540882.pdf). Subject to availability. - CONTACT (Local organizers): Victor Olaya Ferrero (vol...@unex.es), Tomislav Hengl (t.he...@uva.nl) and Juan Carlos Jiménez (jcfer...@unex.es) URL: www.geostat2010.info ___ 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] Kriging concept question
Tobin Cara wrote: Hello, I have recently read an interesting article about integrating Limited Area Models (LAMs) into kriging with external drift for temperature (Libert� et al. link below). www.wmo.int/pages/prog/www/IMOP/.../IOM.../P2(05)_Perini_Italy.doc This URL is incomplete. Please check. As I understand, it seems that the authors generated variograms with the differences between the closest LAM grid point and the ground point. Then, they used kriging with external drift where the ground observation stations were the primary data and the secondary data was the LAM. Conceptually, is it possible to generate a variogram with differences and then krig the data rather than the differences? I am not 100% sure what do you mean by differences. If you mean residuals, then the answer is positive (the variogram HAS to be of the residuals and not of the target variable). KED (or what I prefer to call regression-kriging) in a moving window is in fact probably more general approach than global analysis. For example you can be very relaxed about stationarity assumptions. The problem is that, for each prediction point, you will need to fit both the trend model (LAM?) and estimate a local variogram (using global variogram and local regression models means 'cheating'), which is probably very time consuming! In fact, I am not aware of any applications of KED with moving window (maybe I missed some?). For more info see also section 2.2 in my book: http://spatial-analyst.net/book/Regressionkriging HTH T. Hengl GEOSTAT 2010 - THE Spatial Data Analysis event of the year! http://geostat2010.info I would appreciate anyone's explanations or corrections to my understanding. Thank you, Cara [[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
Re: [R-sig-Geo] raster / polygon overlay
Dear Laure, If this is still actual, here are few tips on how to aggregate continuous information (rasters) given some polygon maps. - # obtain a polygon map e.g. World countries: load(url(http://spatial-analyst.net/DATA/WorldPolyCountries.Rdata;)) str(worldpolycountr...@data) 'data.frame': 246 obs. of 11 variables: $ FIPS : Factor w/ 243 levels AA,AC,AE,..: 2 5 6 7 8 10 11 12 13 17 ... $ ISO2 : Factor w/ 246 levels AD,AE,AF,..: 4 61 17 6 7 9 12 11 14 24 ... $ ISO3 : Factor w/ 246 levels ABW,AFG,AGO,..: 15 64 18 6 11 3 12 10 16 25 ... $ UN : int 28 12 31 8 51 24 16 32 36 48 ... $ NAME : Factor w/ 246 levels Afghanistan,..: 10 4 16 3 12 7 5 11 14 18 ... $ AREA : int 44 238174 8260 2740 2820 124670 20 273669 768230 71 ... $ POP2005 : int 83039 32854159 8352021 3153731 3017661 16095214 64051 38747148 20310208 724788 ... $ REGION : int 19 2 142 150 142 2 9 19 9 142 ... $ SUBREGION: int 29 15 145 39 145 17 61 5 53 145 ... $ LON : num -61.78 2.63 47.4 20.07 44.56 ... $ LAT : num 17.1 28.2 40.4 41.1 40.5 ... writeOGR(WorldPolyCountries, WorldPolyCountries.shp, WorldPolyCountries, ESRI Shapefile) # download some rasters, e.g. global elevations: download.file(http://spatial-analyst.net/worldmaps/globedem.zip;, destfile=paste(getwd(), /glodedem.zip, sep=)) unzip(paste(getwd(), /glodedem.zip, sep=)) rsaga.esri.to.sgrd(in.grids=globedem.asc, out.sgrd=globedem.sgrd, in.path=getwd()) # Get grid Statistics for Polygons using SAGA GIS: # rsaga.get.usage(shapes_grid, 2) rsaga.geoprocessor(lib=shapes_grid, module=2, param=list(GRID=globedem.sgrd, POLY=WorldPolyCountries.shp, RESULT=WorldPolyCountries.shp, QUANTILES=TRUE, QUANTILE_STEP=2)) # To have more control, read maps R: worldgrid - readGDAL(globedem.asc) names(worldgrid) - globedem pixelsize - worldg...@grid@cellsize[1] # convert Polygons to rasters: rsaga.geoprocessor(lib=grid_gridding, module=0, param=list(GRID=WorldCountries.sgrd, INPUT=WorldPolyCountries.shp, FIELD=3, LINE_TYPE=0, TARGET_TYPE=0, USER_CELL_SIZE=pixelsize, user_x_extent_min=worldg...@bbox[1,1]+pixelsize/2, user_x_extent_max=worldg...@bbox[1,2], user_y_extent_min=worldg...@bbox[2,1]+pixelsize/2, user_y_extent_max=worldg...@bbox[2,2]-pixelsize/2)) # read to R: worldgrid$UN - as.factor(readGDAL(WorldCountries.sdat)$band1) # aggregate: globedem.stats - boxplot(worldgrid$globedem ~ worldgrid$UN, plot=FALSE) # match the country names: UN.name - merge(data.frame(UN=as.integer(globedem.stats$names)), worldpolycountr...@data[,c(UN,NAME)], by=UN) globedem.stats$NAME - UN.name$NAME data.frame(elev=globedem.stats$stats[3,], country=globedem.stats$NAME)[order(globedem.stats$stats[3,], decreasing=TRUE)[1:10],] elev country 204 3396 Tajikistan 19 2802 Bhutan 114 2738 Kyrgyzstan 117 2192 Lesotho 6 2082 Andorra 16 1876 Armenia 85 1774 Greenland 174 1620 Rwanda 1 1588 Afghanistan 142 1566 Nepal # 10 highest countries in the World... - If you need to downscale gridded maps, then consider also using SAGA GIS, which has a very flexible down/up-scaling functionality: http://spatial-analyst.net/wiki/index.php?title=Export_maps_to_GE#Rescaling_maps_in_SAGA_GIS HTH T. Hengl http://home.medewerker.uva.nl/t.hengl/ -Original Message- From: r-sig-geo-boun...@stat.math.ethz.ch [mailto:r-sig-geo- boun...@stat.math.ethz.ch] On Behalf Of laure velez Sent: Tuesday, January 05, 2010 2:47 PM To: r-sig-geo@stat.math.ethz.ch Subject: Re: [R-sig-Geo] raster / polygon overlay Thanks again for the fast reply, In fact this problem is only a tiny part of the project and I'll need to run this overlay operation over more than 600 polygons (watersheds) to get mean values (and standard error) of ~10 environmental variables (rasters that are not at the same resolution, i.e 2.5 degrees or 0.1 degrees). Furthermore I will have to overlay the watershed polygons with polygons of landcover. The scale of this project is regional (the whole Mediterranean basin) thus the number (and size) of rasters preclude analyses to be run at the entire scale of the project. Thus I can imagine the following automated process : 1. calculate a buffer around the watershed polygon. 2. subset the grid area that overlay with the buffer (will it work with big grid cells ?) 3. spsample the raster to a finer resolution (the number of new samples will depend on the initial resolution of the raster), what about using the bb argument of the spsample function ? 3. overlay this resampled raster with the watershed polygon. 4. get statistics (mean, sd, ...) from the resulting dataset. As a newbie I would like to get some help (function names) for those different steps (specifically 1 (buffer), 2 (grid subset)). Thanks again for all your useful
Re: [R-sig-Geo] GEOSTAT 2010 Summer School, 27 June - 4 July 2010; University of Plasencia, Spain
Dear R-sig-geo, In relation with our preparations for the GEOSTAT 2010 summer school, it would be of great help if some of you would consider printing and posting the announcement for the summer school via e.g. news board at your department/institute. This way also people off this mailing list (but with potentially high motivation) might consider applying. For more info see: http://geostat2010.info/flyer A PDF of the flyer can be obtained from here: http://geostat2010.info/system/files/GEOSTAT2010_flyer.pdf [1.5 MB] Thank you! T. Hengl http://home.medewerker.uva.nl/t.hengl/ -Original Message- From: r-sig-geo-boun...@stat.math.ethz.ch [mailto:r-sig-geo- boun...@stat.math.ethz.ch] On Behalf Of Tomislav Hengl Sent: Monday, December 28, 2009 9:50 AM To: r-sig-geo Cc: Victor Olaya; Olaf Conrad Subject: [R-sig-Geo] GEOSTAT 2010 Summer School, 27 June - 4 July 2010; University of Plasencia, Spain GEOSTAT 2010 Summer School 27 June - 4 July 2010; University of Plasencia, Spain This is the first call for participation in the GEOSTAT 2010 Summer School. GEOSTAT focuses on important aspects of statistical analysis of spatial and spatio-temporal data using open source / free GIS tools: R, SAGA GIS, GRASS GIS, Quantum GIS, GDAL, Google Earth and similar. The course participants learn how to move data back and forth between different environments, how to produce scripts and automate analysis. We welcome also R beginners and users needing refresh courses in programming. This year, we would also like to introduce/promote topics such as: web-based computing, WPS client-server environments, 3D and 4D geostatistics, combining R+SAGA/GRASS. This is a 5-day course with two parallel sessions, which means that there will be total 7 full-day blocks (three days with parallel sessions) of lectures; the last day of the summer school participants can present their research problems and ask for feedback from the whole summer school. For more info see: http://geostat2010.info/Info - LECTURERS: Roger Bivand (nhh.no), Edzer Pebesma (uni-muenster.de), Gerard Heuvelink (wur.nl), Olaf Conrad (geowiss.uni-hamburg.de), Tomislav Hengl (uva.nl) and Victor Olaya Ferrero (unex.es) - DEADLINE TO REGISTER: 15th February 2010; online at: http://geostat2010.info/deadline_registration - SELECTION OF CANDIDATES: The participants will be selected using the following criteria: 1. your publication record / academic position; 2. previous involvement in the open source activities (R-sig-geo, conferences and workshops) 3. distance from the venue (more distant applications have advantage); 4. time of application (first-come-first-served); - FEES: To be determined; usually between 200-300 EUR (excluding accommodation costs). 5-6 double bed rooms (in-house accommodation) available at subsidized price for participants who reside in, and are nationals of an ODA country (http://www.oecd.org/dataoecd/32/40/43540882.pdf). Subject to availability. - CONTACT (Local organizers): Victor Olaya Ferrero (vol...@unex.es), Tomislav Hengl (t.he...@uva.nl) and Juan Carlos Jiménez (jcfer...@unex.es) URL: www.geostat2010.info ___ 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] Raster Sample Data
Here are several data sets that I try to maintain and improve always: DEMs for testing geomorphometry algorithms: http://geomorphometry.org/content/data-sets The Baranja hill data set includes many rasters of different type: http://geomorphometry.org/content/baranja-hill We used in in our book (http://geomorphometry.org/book). Data sets for testing geostatistical apps: http://spatial-analyst.net/book/data-sets 10 km resolution publicly available world maps: http://spatial-analyst.net/wiki/index.php?title=Global_datasets HTH T. Hengl http://home.medewerker.uva.nl/t.hengl/ Michele Tobias wrote: I'm looking for some raster data to use to with the Raster package (it's for students to work with in a class) and I was wondering if anyone had some suggestions of any sample datasets that come with other R packages. I've found the Meuse.grid data that comes with the SP Package which is a good start, but I'd also like to find some other data. Vegetation reflectance data would be ideal so we can do some NDVI calculations, but anything we can logically run some mathematical operations on would work. It would also be good if the files were of a relatively small size. Downloading full Landsat scenes, for example, will be pretty cumbersome. Thanks for your help! best, Michele PhD Candidate Geography Graduate Group University of California, Davis mmtob...@ucdavis.edu ___ 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] NDVI time-series filter
Dear Jonathan, Interesting topic. Please send us also your code example and if possible part of your dataset (I assume it is a SpatialGridDataFrame?). Operations on lists can be speed up by using lapply and by implementing your own functions, consider also running Rprof to check which operation is using the most time (read more in R inferno or e.g. here: http://manuals.bioinformatics.ucr.edu/home/programming-in-r#Progr_noloops), eventually if the operation is so complex, you simply cannot anticipate to have the results within few minutes (except if you try using super or grid computing; see http://cran.r-project.org/web/views/HighPerformanceComputing.html). HTH, T. Hengl http://home.medewerker.uva.nl/t.hengl/ -Original Message- From: r-sig-geo-boun...@stat.math.ethz.ch [mailto:r-sig-geo- boun...@stat.math.ethz.ch] On Behalf Of Jonathan Thayn Sent: Friday, January 29, 2010 5:02 PM To: r-sig-geo@stat.math.ethz.ch Subject: [R-sig-Geo] NDVI time-series filter I need to smooth a collection of annual NDVI time-series data for a fairly large image. Right now I am using an interruptive method that compares the difference between each point and the mean of its two neighbors. This process repeats until a threshold is reached. It works well, but it takes a long time. When I process the image without smoothing the data it takes about 2 hours – with the smoothing it looks like it will take up to 7 days. The method needs to elevate troughs in the data but retain the peaks – I assume that any troughs are really cloud contamination and that the peaks are accurate NDVI values. Does anyone know of such a smoother that isn't so time consuming? Jonathan B. Thayn, Ph.D. Illinois State University Department of Geography - Geology 200A Felmley Hall Normal, Illinois 61790-4400 (309) 438-8112 jth...@ilstu.edu my.ilstu.edu/~jthayn [[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
