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: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] On Behalf Of [EMAIL PROTECTED] Sent: 2008年12月3日 22:00 To: r-sig-geo@stat.math.ethz.ch Subject: R-sig-Geo Digest, Vol 64, Issue 3 Send R-sig-Geo mailing list submissions to r-sig-geo@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 [EMAIL PROTECTED] You can reach the person managing the list at [EMAIL PROTECTED] 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 <[EMAIL PROTECTED]> Subject: [R-sig-Geo] GSTAT: Optimize power value for IDW To: r-sig-geo@stat.math.ethz.ch Message-ID: <[EMAIL PROTECTED]> 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) [EMAIL PROTECTED] ------------------------------ Message: 2 Date: Tue, 02 Dec 2008 21:06:38 +0100 From: Paul Hiemstra <[EMAIL PROTECTED]> Subject: Re: [R-sig-Geo] GSTAT: Optimize power value for IDW To: Zev Ross <[EMAIL PROTECTED]> Cc: r-sig-geo@stat.math.ethz.ch Message-ID: <[EMAIL PROTECTED]> 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" <[EMAIL PROTECTED]> Subject: [R-sig-Geo] R: GSTAT: Optimize power value for IDW (Paul Hiemstra approach) To: "'Paul Hiemstra'" <[EMAIL PROTECTED]>, "'Zev Ross'" <[EMAIL PROTECTED]> Cc: r-sig-geo@stat.math.ethz.ch Message-ID: <[EMAIL PROTECTED]@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, [EMAIL PROTECTED],1], [EMAIL PROTECTED],2], [EMAIL PROTECTED],1], [EMAIL PROTECTED],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: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] Per conto di Paul Hiemstra Inviato: marted? 2 dicembre 2008 12.07 A: Zev Ross Cc: r-sig-geo@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@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 <[EMAIL PROTECTED]> Subject: Re: [R-sig-Geo] R: GSTAT: Optimize power value for IDW (Paul Hiemstra approach) To: Alessandro <[EMAIL PROTECTED]> Cc: r-sig-geo@stat.math.ethz.ch Message-ID: <[EMAIL PROTECTED]> 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 <[EMAIL PROTECTED]> Subject: Re: [R-sig-Geo] R: GSTAT: Optimize power value for IDW (Paul Hiemstra approach) To: Zev Ross <[EMAIL PROTECTED]> Cc: r-sig-geo@stat.math.ethz.ch Message-ID: <[EMAIL PROTECTED]> 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, >> [EMAIL PROTECTED],1], [EMAIL PROTECTED],2], >> [EMAIL PROTECTED],1], [EMAIL PROTECTED],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: [EMAIL PROTECTED] >> [mailto:[EMAIL PROTECTED] Per conto di Paul Hiemstra >> Inviato: marted? 2 dicembre 2008 12.07 >> A: Zev Ross >> Cc: r-sig-geo@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) > [EMAIL PROTECTED] > ------------------------------------------------------------------------ > > _______________________________________________ > 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 [EMAIL PROTECTED] ------------------------------ Message: 6 Date: Wed, 3 Dec 2008 09:29:29 +0000 From: Corrado <[EMAIL PROTECTED]> Subject: [R-sig-Geo] Raster file from ascii file and flattening Africa .... :) To: [EMAIL PROTECTED], r-sig-geo@stat.math.ethz.ch Message-ID: <[EMAIL PROTECTED]> 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: [EMAIL PROTECTED] ------------------------------ Message: 7 Date: Wed, 3 Dec 2008 10:26:14 -0000 From: "Kamran Safi" <[EMAIL PROTECTED]> Subject: Re: [R-sig-Geo] Raster file from ascii file and flattening Africa ....:) To: "Corrado" <[EMAIL PROTECTED]>, <[EMAIL PROTECTED]>, <r-sig-geo@stat.math.ethz.ch> Message-ID: <[EMAIL PROTECTED]> 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([EMAIL PROTECTED]) 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([EMAIL PROTECTED])) # 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([EMAIL PROTECTED])) { results[i] <- mean(<your spatial points>$values[which(overlay(as.SpatialPolygons.PolygonsList([EMAIL PROTECTED] 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: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] On Behalf Of Corrado Sent: 03 December 2008 09:29 To: [EMAIL PROTECTED]; r-sig-geo@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: [EMAIL PROTECTED] _______________________________________________ R-sig-Geo mailing list R-sig-Geo@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}} _______________________________________________ R-sig-Geo mailing list R-sig-Geo@stat.math.ethz.ch https://stat.ethz.ch/mailman/listinfo/r-sig-geo
