Interpolation option: IDW or OK?
Hi Edzer,
I found a problem with gstat "krige" prediction if I am correct this
time. For example, with OLSENP variable in my dataset, I developed a
theoretical model for log(OLSENP) variogram: model=exp, c0=0.2019,
c1=0.2218 and Range=162 meters. By using this model, I predicted the
OLSENP surface using gstat krige and ArcGIS Spatial Analyst with the
same nmax and found the results are significantly different. gstat
significantly underestimated values (here is where I had an idea to use
IDW), but Spatial Analyst gave very reliable predictions. I would like
you to double check the algorithms behind gstat.
By the way, I start to puzzle because I used gstat to calculate other
variables in my dataset, such as total nitrogen and extractable
potassium and also I used gstat to run a series regression-kriging to
predict soil organic matter distribution and submitted to a journal
already.
I appreciate you could give me your response as soon as possible.
Regards
Yong Li
-----Original Message-----
From: r-sig-geo-bounces at stat.math.ethz.ch
[mailto:r-sig-geo-bounces at stat.math.ethz.ch] On Behalf Of
r-sig-geo-request at stat.math.ethz.ch
Sent: Tuesday, 10 February 2009 10:00 PM
To: r-sig-geo at stat.math.ethz.ch
Subject: R-sig-Geo Digest, Vol 66, Issue 10
Send R-sig-Geo mailing list submissions to
r-sig-geo at stat.math.ethz.ch
To subscribe or unsubscribe via the World Wide Web, visit
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
or, via email, send a message with subject or body 'help' to
r-sig-geo-request at stat.math.ethz.ch
You can reach the person managing the list at
r-sig-geo-owner at stat.math.ethz.ch
When replying, please edit your Subject line so it is more specific
than "Re: Contents of R-sig-Geo digest..."
Today's Topics:
1. error rgdal library if used from within Grass plugin
(Windows) (Agustin Lobo)
2. Re: error rgdal library if used from within Grass plugin
(Windows) (Roger Bivand)
3. Re: FW: Interpolcation option: IDW or OK? (Tomislav Hengl)
4. Calculating descriptive stats from many maps (Rainer M Krug)
5. Re: Calculating descriptive stats from many maps (Tomislav Hengl)
6. Re: Calculating descriptive stats from many maps (Roger Bivand)
7. Interpolcation option: IDW or OK? (Yong Li)
8. Re: Calculating descriptive stats from many maps (Robert Hijmans)
9. Re: FW: Interpolcation option: IDW or OK? (Robert Hijmans)
10. Re: Calculating descriptive stats from many maps (Rainer M Krug)
----------------------------------------------------------------------
Message: 1
Date: Mon, 09 Feb 2009 15:06:24 +0100
From: Agustin Lobo <aloboaleu at gmail.com>
Subject: [R-sig-Geo] error rgdal library if used from within Grass
plugin (Windows)
To: r-sig-geo at stat.math.ethz.ch
Message-ID: <49903860.8060406 at gmail.com>
Content-Type: text/plain; charset=iso-8859-1; format=flowed
Hi!
The following error only occurs if the R (2.7.2) session is started
from a GRASS shell opened through the QGIS GRASS plugin in windows.
The error does not occur If R is started form its own icon
or by double click in the .RData object, but then the spgrass6
package would not find the GRASS environment.
It also works in linux.
The involved commands are:
1. Start QGIS
2. Star GRASS plugin and open mapset
3. Open Grass Shell
4. Run R and execute require(spgrass6)
Loading required package: spgrass6
Loading required package: sp
Loading required package: rgdal
Error in fun(...) :
GDAL Error 1: Can't load requested DLL:
C:\OSGeo4W\bin\gdalplugins\gdal_ECW_JP2ECW.dll
126: N?o foi poss?vel encontrar o m?dulo especificado.
Error : .onLoad failed in 'loadNamespace' for 'rgdal'
Erro: package 'rgdal' could not be loaded
Any help appreciated. Agus -- Dr. Agustin Lobo Institut de Ciencies de la Terra "Jaume Almera" (CSIC) LLuis Sole Sabaris s/n 08028 Barcelona Spain Tel. 34 934095410 Fax. 34 934110012 email: Agustin.Lobo at ija.csic.es http://www.ija.csic.es/gt/obster ------------------------------ Message: 2 Date: Mon, 9 Feb 2009 16:21:25 +0100 (CET) From: Roger Bivand <Roger.Bivand at nhh.no> Subject: Re: [R-sig-Geo] error rgdal library if used from within Grass plugin (Windows) To: Agustin.Lobo at ija.csic.es Cc: r-sig-geo at stat.math.ethz.ch Message-ID: <alpine.LRH.2.00.0902091618410.12628 at reclus.nhh.no> Content-Type: text/plain; charset="iso-8859-1"; Format="flowed"
On Mon, 9 Feb 2009, Agustin Lobo wrote:
Hi! The following error only occurs if the R (2.7.2) session is started from a GRASS shell opened through the QGIS GRASS plugin in windows. The error does not occur If R is started form its own icon or by double click in the .RData object, but then the spgrass6 package would not find the GRASS environment. It also works in linux. The involved commands are: 1. Start QGIS 2. Star GRASS plugin and open mapset 3. Open Grass Shell 4. Run R and execute require(spgrass6)
Loading required package: spgrass6
Loading required package: sp
Loading required package: rgdal
Error in fun(...) :
GDAL Error 1: Can't load requested DLL:
C:\OSGeo4W\bin\gdalplugins\gdal_ECW_JP2ECW.dll
126: N?o foi poss?vel encontrar o m?dulo especificado.
Error : .onLoad failed in 'loadNamespace' for 'rgdal'
Erro: package 'rgdal' could not be loaded
Which rgdal binary from where? This really is a question for osgeo4w, which is bleeding edge at the moment. Are you using the early draft rgdal binary built against a specific osgeo4w version and available from my website? Can we settle this off-list until osgeo4w stabilises? Roger
Any help appreciated. Agus
-- 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: Roger.Bivand at nhh.no ------------------------------ Message: 3 Date: Mon, 9 Feb 2009 16:55:51 +0100 From: "Tomislav Hengl" <T.Hengl at uva.nl> Subject: Re: [R-sig-Geo] FW: Interpolcation option: IDW or OK? To: <r-sig-geo at stat.math.ethz.ch> Message-ID: <ECA350872DC649129EF25F1C51EBCFFA at pcibed193> Content-Type: text/plain; charset="windows-1250" 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.1111/j.1467-9671.2006.01015.x
-----Original Message----- From: r-sig-geo-bounces at stat.math.ethz.ch
[mailto:r-sig-geo-bounces at stat.math.ethz.ch] On Behalf
Of Edzer Pebesma Sent: Monday, February 09, 2009 9:08 AM To: Yong Li Cc: r-sig-geo at 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.pebesma at wwu.de
_______________________________________________ R-sig-Geo mailing list R-sig-Geo at stat.math.ethz.ch https://stat.ethz.ch/mailman/listinfo/r-sig-geo
------------------------------ Message: 4 Date: Mon, 9 Feb 2009 17:58:26 +0200 From: Rainer M Krug <r.m.krug at gmail.com> Subject: [R-sig-Geo] Calculating descriptive stats from many maps To: R-sig-Geo at stat.math.ethz.ch Message-ID: <fb7c7e870902090758n66f987fk122de94146c13fe4 at mail.gmail.com> Content-Type: text/plain; charset=ISO-8859-1 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 ------------------------------ Message: 5 Date: Mon, 9 Feb 2009 17:33:51 +0100 From: "Tomislav Hengl" <T.Hengl at uva.nl> Subject: Re: [R-sig-Geo] Calculating descriptive stats from many maps To: <R-sig-Geo at stat.math.ethz.ch> Message-ID: <D0BA11FFBCAD490B9D79F65EAC48FDE9 at pcibed193> Content-Type: text/plain; charset="windows-1250" 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(maps at data, na.rm=T, dims=1) maps$avg <- maps$Nsum/(length(names(meuse.grid at data))-1)
You could also loop the sd, mean and quantile function over a range of cells:
for(i in length(names(maps at data))) {
maps at data$sd[i] <- sd(maps at data[i,])
maps at data$mean[i] <- mean(maps at 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-bounces at stat.math.ethz.ch
[mailto:r-sig-geo-bounces at stat.math.ethz.ch] On Behalf
Of Rainer M Krug Sent: Monday, February 09, 2009 4:58 PM To: R-sig-Geo at 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 at stat.math.ethz.ch https://stat.ethz.ch/mailman/listinfo/r-sig-geo
------------------------------ Message: 6 Date: Mon, 9 Feb 2009 18:05:25 +0100 (CET) From: Roger Bivand <Roger.Bivand at nhh.no> Subject: Re: [R-sig-Geo] Calculating descriptive stats from many maps To: Tomislav Hengl <T.Hengl at uva.nl> Cc: R-sig-Geo at stat.math.ethz.ch Message-ID: <alpine.LRH.2.00.0902091800170.15636 at reclus.nhh.no> Content-Type: TEXT/PLAIN; charset=US-ASCII; format=flowed
On Mon, 9 Feb 2009, Tomislav Hengl wrote:
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(maps at data, na.rm=T, dims=1) maps$avg <- maps$Nsum/(length(names(meuse.grid at data))-1)
You could also loop the sd, mean and quantile function over a range of
cells:
for(i in length(names(maps at data))) {
maps at data$sd[i] <- sd(maps at data[i,])
maps at data$mean[i] <- mean(maps at data[i,])
...
}
This could take a lot of time!
Tom, Rainer, Yes, using sapply(slot(maps, "data"), summary) or similar, you get the map-wise statistics. But have I misunderstood, or are the statistics in question cell-wise? This would involve stacking subset areas for all 25' maps, wouldn't it? Brutally, a loop in readGDAL() from rgdal with offset= and region.dim= shifted? Is there a canned way to do this in the R-forge raster package (by the way, regularly one of the R-forge packages showing most developer activity)? Roger
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-bounces at stat.math.ethz.ch
[mailto:r-sig-geo-bounces at stat.math.ethz.ch] On Behalf
Of Rainer M Krug Sent: Monday, February 09, 2009 4:58 PM To: R-sig-Geo at 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 at stat.math.ethz.ch https://stat.ethz.ch/mailman/listinfo/r-sig-geo
_______________________________________________ R-sig-Geo mailing list R-sig-Geo at stat.math.ethz.ch https://stat.ethz.ch/mailman/listinfo/r-sig-geo
-- 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: Roger.Bivand at nhh.no ------------------------------ Message: 7 Date: Tue, 10 Feb 2009 10:40:40 +1100 From: Yong Li <yong.li at unimelb.edu.au> Subject: [R-sig-Geo] Interpolcation option: IDW or OK? To: Edzer Pebesma <edzer.pebesma at uni-muenster.de> Cc: r-sig-geo at stat.math.ethz.ch Message-ID: <86DBA0678E017341B449A62F258E29561549C6 at IS-EX-BEV3.unimelb.edu.au> Content-Type: text/plain; charset=us-ascii Many thanks for the message, Edzer. -----Original Message----- From: Edzer Pebesma [mailto:edzer.pebesma at uni-muenster.de] Sent: Monday, 9 February 2009 7:08 PM To: Yong Li Cc: r-sig-geo at stat.math.ethz.ch Subject: Re: FW: [R-sig-Geo] 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. YONGLI================================================================== ==== Could you have a look at the predictions generated by my nested model and compare to the observations? The values were severely underestimated. YONGLI================================================================== ==== 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?
YONGLI==================================================================
====
Could you have a look at "Mueller et al. (2004) Map Quality for Ordinary
Kriging and Inverse Distance Weighted Interpolation. Soil Sci. Soc. Am.
J. 68:2042-2047."?
YONGLI==================================================================
====
Regards,
Yong Li
------------------------------
Message: 8
Date: Tue, 10 Feb 2009 10:07:11 +0800
From: Robert Hijmans <r.hijmans at gmail.com>
Subject: Re: [R-sig-Geo] Calculating descriptive stats from many maps
To: R-sig-Geo at stat.math.ethz.ch
Message-ID:
<dc22b2570902091807m4698caf7xbe79e562da35bce2 at mail.gmail.com>
Content-Type: text/plain; charset=ISO-8859-1
Dear Rainer,
This is how can you can do it with the raster package
# install.packages("raster", repos="http://R-Forge.R-project.org")
require(raster)
# Try it for a few files first..
n <- 10
# create a list (or vector) of file names, e.g. :
fn <- list()
for (i in 1:n) { fn[i] <- paste('myfile', i, '.tif', sep='') }
# make a RasterStack
s <- stack(fn)
r1 <- mCalc(s, fun=mean)
r2 <- mCalc(s, fun=sd)
#r can be plotted, coerced to sp objects, etc.
plot(r1)
# or saved to file
r1 <- setFilename(r1, 'cellmeans.tif')
r1 <- writeRaster(r1, format='GTiff')
Robert
On Tue, Feb 10, 2009 at 1:05 AM, Roger Bivand <Roger.Bivand at nhh.no>
wrote:
On Mon, 9 Feb 2009, Tomislav Hengl wrote:
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(maps at data, na.rm=T, dims=1) maps$avg <- maps$Nsum/(length(names(meuse.grid at data))-1)
You could also loop the sd, mean and quantile function over a range
of
cells:
for(i in length(names(maps at data))) {
maps at data$sd[i] <- sd(maps at data[i,])
maps at data$mean[i] <- mean(maps at data[i,])
...
}
This could take a lot of time!
Tom, Rainer, Yes, using sapply(slot(maps, "data"), summary) or similar, you get the map-wise statistics. But have I misunderstood, or are the statistics
in
question cell-wise? This would involve stacking subset areas for all
25'
maps, wouldn't it? Brutally, a loop in readGDAL() from rgdal with
offset=
and region.dim= shifted? Is there a canned way to do this in the
R-forge
raster package (by the way, regularly one of the R-forge packages
showing
most developer activity)? Roger
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-bounces at stat.math.ethz.ch [mailto:r-sig-geo-bounces at stat.math.ethz.ch] On Behalf Of Rainer M Krug Sent: Monday, February 09, 2009 4:58 PM To: R-sig-Geo at 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 at stat.math.ethz.ch https://stat.ethz.ch/mailman/listinfo/r-sig-geo
_______________________________________________ R-sig-Geo mailing list R-sig-Geo at stat.math.ethz.ch https://stat.ethz.ch/mailman/listinfo/r-sig-geo
-- 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: Roger.Bivand at nhh.no
_______________________________________________ R-sig-Geo mailing list R-sig-Geo at stat.math.ethz.ch https://stat.ethz.ch/mailman/listinfo/r-sig-geo
------------------------------ Message: 9 Date: Tue, 10 Feb 2009 12:44:57 +0800 From: Robert Hijmans <r.hijmans at gmail.com> Subject: Re: [R-sig-Geo] FW: Interpolcation option: IDW or OK? To: Tomislav Hengl <T.Hengl at uva.nl> Cc: r-sig-geo at stat.math.ethz.ch Message-ID: <dc22b2570902092044w29078ad6yad92a971f376b662 at mail.gmail.com> Content-Type: text/plain; charset=ISO-8859-1 Why not use cross-validation to empirically determine which method performs best for this dataset (in addition to asking if they are better than a random draw)? Robert 2009/2/9 Tomislav Hengl <T.Hengl at uva.nl>:
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.1111/j.1467-9671.2006.01015.x
-----Original Message----- From: r-sig-geo-bounces at stat.math.ethz.ch
[mailto:r-sig-geo-bounces at stat.math.ethz.ch] On Behalf
Of Edzer Pebesma Sent: Monday, February 09, 2009 9:08 AM To: Yong Li Cc: r-sig-geo at 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.pebesma at wwu.de
_______________________________________________ R-sig-Geo mailing list R-sig-Geo at stat.math.ethz.ch https://stat.ethz.ch/mailman/listinfo/r-sig-geo
_______________________________________________ R-sig-Geo mailing list R-sig-Geo at stat.math.ethz.ch https://stat.ethz.ch/mailman/listinfo/r-sig-geo
------------------------------ Message: 10 Date: Tue, 10 Feb 2009 10:29:25 +0200 From: Rainer M Krug <r.m.krug at gmail.com> Subject: Re: [R-sig-Geo] Calculating descriptive stats from many maps To: Robert Hijmans <r.hijmans at gmail.com> Cc: R-sig-Geo at stat.math.ethz.ch Message-ID: <fb7c7e870902100029kf6e436du542afd36c4067f2f at mail.gmail.com> Content-Type: text/plain; charset=ISO-8859-1 Thanks a lot to all of you. You are right Roger, I need cell-wise statistics I like the idea of the "raster" package, and I will try it out just now. Concerning SAGA: I'll look into that if "raster" does not work (or is to slow). I'll report back Rainer On Tue, Feb 10, 2009 at 4:07 AM, Robert Hijmans <r.hijmans at gmail.com> wrote:
Dear Rainer,
This is how can you can do it with the raster package
# install.packages("raster", repos="http://R-Forge.R-project.org")
require(raster)
# Try it for a few files first..
n <- 10
# create a list (or vector) of file names, e.g. :
fn <- list()
for (i in 1:n) { fn[i] <- paste('myfile', i, '.tif', sep='') }
# make a RasterStack
s <- stack(fn)
r1 <- mCalc(s, fun=mean)
r2 <- mCalc(s, fun=sd)
#r can be plotted, coerced to sp objects, etc.
plot(r1)
# or saved to file
r1 <- setFilename(r1, 'cellmeans.tif')
r1 <- writeRaster(r1, format='GTiff')
Robert
On Tue, Feb 10, 2009 at 1:05 AM, Roger Bivand <Roger.Bivand at nhh.no>
wrote:
On Mon, 9 Feb 2009, Tomislav Hengl wrote:
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(maps at data, na.rm=T, dims=1) maps$avg <- maps$Nsum/(length(names(meuse.grid at data))-1)
You could also loop the sd, mean and quantile function over a range
of
cells:
for(i in length(names(maps at data))) {
maps at data$sd[i] <- sd(maps at data[i,])
maps at data$mean[i] <- mean(maps at data[i,])
...
}
This could take a lot of time!
Tom, Rainer, Yes, using sapply(slot(maps, "data"), summary) or similar, you get
the
map-wise statistics. But have I misunderstood, or are the statistics
in
question cell-wise? This would involve stacking subset areas for all
25'
maps, wouldn't it? Brutally, a loop in readGDAL() from rgdal with
offset=
and region.dim= shifted? Is there a canned way to do this in the
R-forge
raster package (by the way, regularly one of the R-forge packages
showing
most developer activity)? Roger
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-bounces at stat.math.ethz.ch [mailto:r-sig-geo-bounces at stat.math.ethz.ch] On Behalf Of Rainer M Krug Sent: Monday, February 09, 2009 4:58 PM To: R-sig-Geo at 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 at stat.math.ethz.ch https://stat.ethz.ch/mailman/listinfo/r-sig-geo
_______________________________________________ R-sig-Geo mailing list R-sig-Geo at stat.math.ethz.ch https://stat.ethz.ch/mailman/listinfo/r-sig-geo
-- 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: Roger.Bivand at nhh.no
_______________________________________________ R-sig-Geo mailing list R-sig-Geo at stat.math.ethz.ch https://stat.ethz.ch/mailman/listinfo/r-sig-geo
_______________________________________________ R-sig-Geo mailing list R-sig-Geo at stat.math.ethz.ch https://stat.ethz.ch/mailman/listinfo/r-sig-geo
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 at stat.math.ethz.ch https://stat.ethz.ch/mailman/listinfo/r-sig-geo End of R-sig-Geo Digest, Vol 66, Issue 10