[R-sig-Geo] How to pick up cell values of grids using points?
Yong Li
yong.li at unimelb.edu.au
Wed Dec 3 13:13:34 CET 2008
Hi ALL,
I want to pick up cells values of a grid (say a RS image) using points' positions. For example, I have an image and a point shape file with the same coordinate system, then I try to overlay points on the image, and want to get the cell values of the on-point cell and 8 surrounding cells. Could any expert give some tips to solve my problem?
Regards
Yong
-----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: 2008年12月3日 22:00
To: r-sig-geo at stat.math.ethz.ch
Subject: R-sig-Geo Digest, Vol 64, Issue 3
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. GSTAT: Optimize power value for IDW (Zev Ross)
2. Re: GSTAT: Optimize power value for IDW (Paul Hiemstra)
3. R: GSTAT: Optimize power value for IDW (Paul Hiemstra
approach) (Alessandro)
4. Re: R: GSTAT: Optimize power value for IDW (Paul Hiemstra
approach) (Zev Ross)
5. Re: R: GSTAT: Optimize power value for IDW (Paul Hiemstra
approach) (Edzer Pebesma)
6. Raster file from ascii file and flattening Africa .... :)
(Corrado)
7. Re: Raster file from ascii file and flattening Africa ....:)
(Kamran Safi)
----------------------------------------------------------------------
Message: 1
Date: Tue, 02 Dec 2008 13:59:47 -0500
From: Zev Ross <zev at zevross.com>
Subject: [R-sig-Geo] GSTAT: Optimize power value for IDW
To: r-sig-geo at stat.math.ethz.ch
Message-ID: <493585A3.2040909 at zevross.com>
Content-Type: text/plain; charset=ISO-8859-1; format=flowed
Hi All,
ArcGIS has a nice little button that calculates the optimal power value
to use for inverse distance weighting based on cross-validation and
RMSPE. Just wondering if anyone had written something similar in R --
I'm using GSTAT and I'd like to avoid back and forth with ArcGIS (and
obviously I'd like to avoid writing it myself as well!).
Zev
--
Zev Ross
ZevRoss Spatial Analysis
120 N Aurora, Suite 3A
Ithaca, NY 14850
607-277-0004 (phone)
866-877-3690 (fax, toll-free)
zev at zevross.com
------------------------------
Message: 2
Date: Tue, 02 Dec 2008 21:06:38 +0100
From: Paul Hiemstra <p.hiemstra at geo.uu.nl>
Subject: Re: [R-sig-Geo] GSTAT: Optimize power value for IDW
To: Zev Ross <zev at zevross.com>
Cc: r-sig-geo at stat.math.ethz.ch
Message-ID: <4935954E.7070705 at geo.uu.nl>
Content-Type: text/plain; charset=ISO-8859-1; format=flowed
Zev Ross schreef:
> Hi All,
>
> ArcGIS has a nice little button that calculates the optimal power
> value to use for inverse distance weighting based on cross-validation
> and RMSPE. Just wondering if anyone had written something similar in R
> -- I'm using GSTAT and I'd like to avoid back and forth with ArcGIS
> (and obviously I'd like to avoid writing it myself as well!).
>
> Zev
>
Hi,
I don't have any code lying around, but a brute force optimization
approach should be quite easy. Also because the range of idw-powers is
relatively small. The speed would ofcourse depend on the size of the
dataset. In code it would look something like (actually works :)):
library(gstat)
data(meuse)
coordinates(meuse) = ‾x+y
# make function to do the cv
do_cv = function(idp) {
meuse_idw = gstat(id = "zn", formula = zinc‾1, data = meuse, set =
list(idp = idp))
out = gstat.cv(meuse_idw)
return(sqrt(mean(out$residual^2))) # RMSE as optimization criterium
}
idw_pow = seq(0.1,5, by = 0.2) # the idwpower values that will be checked
cv_vals = sapply(idw_pow, do_cv) # calculate the rmse
# List of outcomes
print(data.frame(idp = idw_pow, cv_rmse = cv_vals))
After this you select the idw value with the smallest RMSE.
cheers and hth,
Paul
--
Drs. Paul Hiemstra
Department of Physical Geography
Faculty of Geosciences
University of Utrecht
Heidelberglaan 2
P.O. Box 80.115
3508 TC Utrecht
Phone: +31302535773
Fax: +31302531145
http://intamap.geo.uu.nl/‾paul
------------------------------
Message: 3
Date: Tue, 2 Dec 2008 12:44:08 -0800
From: "Alessandro" <alessandro.montaghi at unifi.it>
Subject: [R-sig-Geo] R: GSTAT: Optimize power value for IDW (Paul
Hiemstra approach)
To: "'Paul Hiemstra'" <p.hiemstra at geo.uu.nl>, "'Zev Ross'"
<zev at zevross.com>
Cc: r-sig-geo at stat.math.ethz.ch
Message-ID: <010401c954be$baacdc60$30069520$@montaghi at unifi.it>
Content-Type: text/plain; charset="iso-8859-1"
Hi
Normally I use the R+SAGA to calculate the IDW and create a raster, with
this follow code. I change the radius input with a loop formula to create
several raster and after check the best result (I am studying a oak forest
with low density)
radii.list <- as.list(c(5, 10, 15, 20, 25, 30))
for(i in 1:length(radii.list)){
rsaga.geoprocessor(lib="grid_gridding", module=0,
param=list(GRID=set.file.extension(paste("DEM_iw",radii.list[[i]],sep=""),".
sgrd"),
SHAPES="ground.shp", FIELD=0, TARGET= 0, POWER=1.0,
RADIUS=radii.list[[i]], NPOINTS=20, USER_CELL_SIZE=dem.pixelsize,
USER_X_EXTENT_MIN=ground at bbox[1,1], USER_X_EXTENT_MAX=ground at bbox[1,2],
USER_Y_EXTENT_MIN=ground at bbox[2,1], USER_Y_EXTENT_MAX=ground at bbox[2,2]))
}
After I have 6 raster (DEM_idw_5.sgrd, DEM_idw_10.sgrd, DEM_idw_15.sgrd,
DEM_idw_20.sgrd, etc. etc.). I check hand-made the best IDW. I don't know is
It possible to improve this formula with RMSE in gstat and calculate the
best power.
Ale
-----Messaggio originale-----
Da: r-sig-geo-bounces at stat.math.ethz.ch
[mailto:r-sig-geo-bounces at stat.math.ethz.ch] Per conto di Paul Hiemstra
Inviato: marted? 2 dicembre 2008 12.07
A: Zev Ross
Cc: r-sig-geo at stat.math.ethz.ch
Oggetto: Re: [R-sig-Geo] GSTAT: Optimize power value for IDW
Zev Ross schreef:
> Hi All,
>
> ArcGIS has a nice little button that calculates the optimal power
> value to use for inverse distance weighting based on cross-validation
> and RMSPE. Just wondering if anyone had written something similar in R
> -- I'm using GSTAT and I'd like to avoid back and forth with ArcGIS
> (and obviously I'd like to avoid writing it myself as well!).
>
> Zev
>
Hi,
I don't have any code lying around, but a brute force optimization
approach should be quite easy. Also because the range of idw-powers is
relatively small. The speed would ofcourse depend on the size of the
dataset. In code it would look something like (actually works :)):
library(gstat)
data(meuse)
coordinates(meuse) = ‾x+y
# make function to do the cv
do_cv = function(idp) {
meuse_idw = gstat(id = "zn", formula = zinc‾1, data = meuse, set =
list(idp = idp))
out = gstat.cv(meuse_idw)
return(sqrt(mean(out$residual^2))) # RMSE as optimization criterium
}
idw_pow = seq(0.1,5, by = 0.2) # the idwpower values that will be checked
cv_vals = sapply(idw_pow, do_cv) # calculate the rmse
# List of outcomes
print(data.frame(idp = idw_pow, cv_rmse = cv_vals))
After this you select the idw value with the smallest RMSE.
cheers and hth,
Paul
--
Drs. Paul Hiemstra
Department of Physical Geography
Faculty of Geosciences
University of Utrecht
Heidelberglaan 2
P.O. Box 80.115
3508 TC Utrecht
Phone: +31302535773
Fax: +31302531145
http://intamap.geo.uu.nl/‾paul
_______________________________________________
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: Tue, 02 Dec 2008 15:54:35 -0500
From: Zev Ross <zev at zevross.com>
Subject: Re: [R-sig-Geo] R: GSTAT: Optimize power value for IDW (Paul
Hiemstra approach)
To: Alessandro <alessandro.montaghi at unifi.it>
Cc: r-sig-geo at stat.math.ethz.ch
Message-ID: <4935A08B.6090101 at zevross.com>
Content-Type: text/plain; charset="us-ascii"
An HTML attachment was scrubbed...
URL: <https://stat.ethz.ch/pipermail/r-sig-geo/attachments/20081202/1fcf93f3/attachment-0001.html>
------------------------------
Message: 5
Date: Tue, 02 Dec 2008 22:29:31 +0100
From: Edzer Pebesma <edzer.pebesma at uni-muenster.de>
Subject: Re: [R-sig-Geo] R: GSTAT: Optimize power value for IDW (Paul
Hiemstra approach)
To: Zev Ross <zev at zevross.com>
Cc: r-sig-geo at stat.math.ethz.ch
Message-ID: <4935A8BB.3010900 at uni-muenster.de>
Content-Type: text/plain; charset=ISO-8859-1; format=flowed
Hi Zev,
There's no need to brute force, as optimize is there to help you -- my
guess is that the function is convex. The following takes a while:
> f = function(idp, formula, data,...)
sum(krige.cv(formula,data,set=list(debug=0,idp=idp),...)$residual**2)
> optimize(f, interval=c(0.01,4), formula=log(zinc)‾1, data=meuse)
$minimum
[1] 3.532051
$objective
[1] 32.09126
but that is probably due to the (my?) hopelessly inefficient (though
flexible) setup of krige.cv. Speeding up can be done by allowing a
larger tolerance, or passing e.g. nfold=5 to optimize(). This nfold will
also result in somewhat random output, as it randomly folds the data in
5 partitions.
--
Edzer
Zev Ross wrote:
> Hi All,
>
> Thanks to Paul and Alessandro for their suggestions. Paul's code
> (brute force) worked well for me and the results match up well with
> ArcGIS. I'm not using a large dataset so the speed isn't an issue but
> with a larger dataset it would be. In ArcGIS the optimization is
> instantaneous so clearly the software is doing something different
> than testing out all possible values.
>
> I used Paul's code with no substantive changes then it's
> straightforward to use:
>
> plot(idw_pow, cv_vals)
> idw_pow[which.min(cv_vals)]
>
> And pull out the min.
>
> Thanks for the advice! Zev
>
>
>
>
> Alessandro wrote:
>> Hi
>>
>> Normally I use the R+SAGA to calculate the IDW and create a raster, with
>> this follow code. I change the radius input with a loop formula to create
>> several raster and after check the best result (I am studying a oak forest
>> with low density)
>>
>> radii.list <- as.list(c(5, 10, 15, 20, 25, 30))
>> for(i in 1:length(radii.list)){
>> rsaga.geoprocessor(lib="grid_gridding", module=0,
>> param=list(GRID=set.file.extension(paste("DEM_iw",radii.list[[i]],sep=""),".
>> sgrd"),
>> SHAPES="ground.shp", FIELD=0, TARGET= 0, POWER=1.0,
>> RADIUS=radii.list[[i]], NPOINTS=20, USER_CELL_SIZE=dem.pixelsize,
>> USER_X_EXTENT_MIN=ground at bbox[1,1], USER_X_EXTENT_MAX=ground at bbox[1,2],
>> USER_Y_EXTENT_MIN=ground at bbox[2,1], USER_Y_EXTENT_MAX=ground at bbox[2,2]))
>> }
>>
>>
>> After I have 6 raster (DEM_idw_5.sgrd, DEM_idw_10.sgrd, DEM_idw_15.sgrd,
>> DEM_idw_20.sgrd, etc. etc.). I check hand-made the best IDW. I don't know is
>> It possible to improve this formula with RMSE in gstat and calculate the
>> best power.
>>
>> Ale
>>
>>
>>
>> -----Messaggio originale-----
>> Da: r-sig-geo-bounces at stat.math.ethz.ch
>> [mailto:r-sig-geo-bounces at stat.math.ethz.ch] Per conto di Paul Hiemstra
>> Inviato: marted? 2 dicembre 2008 12.07
>> A: Zev Ross
>> Cc: r-sig-geo at stat.math.ethz.ch
>> Oggetto: Re: [R-sig-Geo] GSTAT: Optimize power value for IDW
>>
>> Zev Ross schreef:
>>
>>> Hi All,
>>>
>>> ArcGIS has a nice little button that calculates the optimal power
>>> value to use for inverse distance weighting based on cross-validation
>>> and RMSPE. Just wondering if anyone had written something similar in R
>>> -- I'm using GSTAT and I'd like to avoid back and forth with ArcGIS
>>> (and obviously I'd like to avoid writing it myself as well!).
>>>
>>> Zev
>>>
>>>
>> Hi,
>>
>> I don't have any code lying around, but a brute force optimization
>> approach should be quite easy. Also because the range of idw-powers is
>> relatively small. The speed would ofcourse depend on the size of the
>> dataset. In code it would look something like (actually works :)):
>>
>> library(gstat)
>> data(meuse)
>> coordinates(meuse) = ‾x+y
>>
>> # make function to do the cv
>> do_cv = function(idp) {
>> meuse_idw = gstat(id = "zn", formula = zinc‾1, data = meuse, set =
>> list(idp = idp))
>> out = gstat.cv(meuse_idw)
>> return(sqrt(mean(out$residual^2))) # RMSE as optimization criterium
>> }
>>
>> idw_pow = seq(0.1,5, by = 0.2) # the idwpower values that will be checked
>> cv_vals = sapply(idw_pow, do_cv) # calculate the rmse
>> # List of outcomes
>> print(data.frame(idp = idw_pow, cv_rmse = cv_vals))
>>
>> After this you select the idw value with the smallest RMSE.
>>
>> cheers and hth,
>> Paul
>>
>>
>
> --
> Zev Ross
> ZevRoss Spatial Analysis
> 120 N Aurora, Suite 3A
> Ithaca, NY 14850
> 607-277-0004 (phone)
> 866-877-3690 (fax, toll-free)
> zev at zevross.com
> ------------------------------------------------------------------------
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at 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.pebesma at wwu.de
------------------------------
Message: 6
Date: Wed, 3 Dec 2008 09:29:29 +0000
From: Corrado <ct529 at york.ac.uk>
Subject: [R-sig-Geo] Raster file from ascii file and flattening Africa
.... :)
To: r-sig-ecology at r-project.org, r-sig-geo at stat.math.ethz.ch
Message-ID: <200812030929.29575.ct529 at york.ac.uk>
Content-Type: text/plain; charset="utf-8"
Dear friends,
I am a kind of advanced newbie, if that makes sense.
I have a text file of the form
coordinate x,coordinate y,cat={real number between 250 and 450}
where coordinate are expressed in latitude and longitude. The files represents
measurements of the size of a skulls on sites all over Africa.
>From it, I would like to build a raster file, 100 km by 100km. There are 2
problems:
1) Unfortunately, in some 100km x 100km squares, there is one of the points
whilst in others there are maybe 20. How do I average, so that in each square
I only have 1 value representing the average?
2) How do we "flatten" Africa so that we may use 100km x 100km squares instead
of 1 degree x 1 degree, without committing a geographical crime? What we need
is to respect the areas ....
Best regards and apologies for the silliness of the questions.
--
Corrado Topi
Global Climate Change & Biodiversity Indicators
Area 18,Department of Biology
University of York, York, YO10 5YW, UK
Phone: + 44 (0) 1904 328645, E-mail: ct529 at york.ac.uk
------------------------------
Message: 7
Date: Wed, 3 Dec 2008 10:26:14 -0000
From: "Kamran Safi" <Kamran.Safi at ioz.ac.uk>
Subject: Re: [R-sig-Geo] Raster file from ascii file and flattening
Africa ....:)
To: "Corrado" <ct529 at york.ac.uk>, <r-sig-ecology at r-project.org>,
<r-sig-geo at stat.math.ethz.ch>
Message-ID: <41E1ED29E5E8E34BBDD8B82CFA1A9D04062E6357 at zsl26>
Content-Type: text/plain; charset="us-ascii"
Hi Corrado,
Being a advanced newbie myself, I first of all understand what you mean
by that and secondly ask you to qualify my answer.
I would, tackling your problem, create a raster polygon in a metric
equal area projection, such as Mollweide. Then you use overlay() and get
for each polygon the set of points that are within each raster polygon.
You need to import the xyz file in R and convert it into a Spatialpoints
data frame.
Here's the first bit
This reads a coast line shapefile and extracts africa from it. Then uses
the boundings boxes to produce a grid at the extent of africa. Then it
projects that raster back to longlat for the overlay() procedure.
map <- readShapePoly("E:/Science/continent.shp", ID="CONTINENT",
proj4string = CRS("+proj=longlat"))
africa <- as.SpatialPolygons.PolygonsList(map at polygons[1])
africa.proj <- spTransform(africa, CRS("+proj=moll"))
grd <- GridTopology(c(bbox(africa.proj)[1,1]+5000,
bbox(africa.proj)[2,1]+50000), c(100000,100000),
c(ceiling((bbox(africa.proj)[1,2] - bbox(africa.proj)[1,1]) /
100000),ceiling((bbox(africa.proj)[1,2] - bbox(africa.proj)[1,1]) /
100000)))
# if you should not have a coastline of africa:
# these are the values you'll need to produce the raster you need to
proceed
# bbox(africa.proj)
# min max
# r1 -2472164 6131319
# r2 -4202811 4490010
polys.proj <- as.SpatialPolygons.GridTopology(grd)
proj4string(polys.proj) <- CRS("+proj=moll")
polys <- spTransform(polys.proj, CRS("+proj=longlat"))
# now you have a spatialpolygon in longlat proj that has equal area and
a
# cell size of 100km2
#prepare a list for you results
results <- rep(NA, length(polys at polygons))
# use something like this to calculate the values per raster grid
# this here is not working probably, as I didn't think about it all too
much
# I have though somewhere a code lying around doing exactly this step,
so if
# you don't succeed, let me know and I dig that out
for(i in 1:length(polys at polygons))
{
results[i] <- mean(<your spatial
points>$values[which(overlay(as.SpatialPolygons.PolygonsList(map.Proj at po
lygons[i]), <your Spatialpoints>) != NA, ])
}
Apart from the final bits, I tested the start and it worked. The last
bit should not be too difficult to solve. The whole thing could be more
elegant by excluding those polygons that are in the sea. But I think
that is something others can try to get their heads around it. Shouldn't
be too difficult: get the centroid coordinates, overlay it with the
costlines of Africa and convert it back into a grid...
Hope this helped.
Kamran
------------------------
Kamran Safi
Postdoctoral Research Fellow
Institute of Zoology
Zoological Society of London
Regent's Park
London NW1 4RY
http://www.zoo.cam.ac.uk/ioz/people/safi.htm
http://spatialr.googlepages.com
http://asapi.wetpaint.com
-----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 Corrado
Sent: 03 December 2008 09:29
To: r-sig-ecology at r-project.org; r-sig-geo at stat.math.ethz.ch
Subject: [R-sig-Geo] Raster file from ascii file and flattening Africa
....:)
Dear friends,
I am a kind of advanced newbie, if that makes sense.
I have a text file of the form
coordinate x,coordinate y,cat={real number between 250 and 450}
where coordinate are expressed in latitude and longitude. The files
represents
measurements of the size of a skulls on sites all over Africa.
>From it, I would like to build a raster file, 100 km by 100km. There
are 2
problems:
1) Unfortunately, in some 100km x 100km squares, there is one of the
points
whilst in others there are maybe 20. How do I average, so that in each
square
I only have 1 value representing the average?
2) How do we "flatten" Africa so that we may use 100km x 100km squares
instead
of 1 degree x 1 degree, without committing a geographical crime? What we
need
is to respect the areas ....
Best regards and apologies for the silliness of the questions.
--
Corrado Topi
Global Climate Change & Biodiversity Indicators
Area 18,Department of Biology
University of York, York, YO10 5YW, UK
Phone: + 44 (0) 1904 328645, E-mail: ct529 at york.ac.uk
_______________________________________________
R-sig-Geo mailing list
R-sig-Geo at stat.math.ethz.ch
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Click
https://www.mailcontrol.com/sr/wQw0zmjPoHdJTZGyOCrrhg==
rtsaKomrEVtGO6LLtLGhCXg+F32PftV4uyzpFBU8KFm0g== to report this email as
spam.
The Zoological Society of London is incorporated by Royal Charter
Principal Office England. Company Number RC000749
Registered address:
Regent's Park, London, England NW1 4RY
Registered Charity in England and Wales no. 208728
_________________________________________________________________________
This e-mail has been sent in confidence to the named add...{{dropped:15}}
More information about the R-sig-Geo
mailing list