[R-sig-Geo] FW: Interpolation option: IDW or OK?
Yong Li
yong.li at unimelb.edu.au
Wed Feb 11 12:50:04 CET 2009
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
More information about the R-sig-Geo
mailing list