Thanks for the quick reply! My replies are in-line below...
On Apr 14, 2010, at 5:14 PM, Edzer Pebesma wrote: > Hi, > > telling us what the drop-down menus tell you leaves me a bit of guess > work. Universal probably means universal kriging; linear drift a first > order trend in the coordinates (pass degree = 1 to krige, and something > like value~x+y to variogram, when long/lat are named x and y). Does the > first "Linear" in "Linear with linear drift" refer to a fitted linear > variogram model? I think so; see the copy-paste of the documentation below... > > Prof. Jeffrey Cardille wrote: >> Hello all, >> >> This is my first posting to R-sig-geo; I looked for the answer in the >> archives, but have not had any luck for my particular question. >> >> I have a newbie question about the correspondence between gstat and ArcGIS. >> I have a problem that gives very very satisfying results in ArcGIS, and I >> would like to write a public function in R for others to use that does the >> same thing. The choice I used in Arc was "Universal" and in the dropdown >> menu it says "Linear with Linear Drift". My data matrix is pretty standard: >> 2000 x 2000, with some NA, and values between 0 and 255. >> >> So I have a few questions. Maybe I'll enumerate them. >> 1. Has anyone made a comprehensive list between ArcGIS functionality and how >> to do it in gstat? > > Great idea! > >> 2. If not: does anyone know a simple correspondence between what ArcGIS is >> doing and how to set up a gstat call for "Universal/Linear with Linear >> Drift"? I have looked online for precise details about Arc's behavior, but >> haven't seen anything detailed enough. > > Well, you have a license, doesn't that come with access to some sort of > documentation? Here is what I can see in the help text from ArcGIS: It's a bit fragmentary because there are two distinct ways to execute what seems to be the same code-- once through the basic ArcGIS functionality, and one through the add-on package "Geostatistical Analyst". So this text is spliced from a few different web pages. Does this help to clarify the correspondence to the gstat parameters? > The Semivariogram Properties dialog box has several models to choose from. > When the Kriging method is set to Ordinary, the available models are > Spherical, Circular, Exponential, Gaussian, and Linear. When the Kriging > method is set to Universal, the available models are Linear with linear drift > and Linear with quadratic drift. > > Universal kriging assumes that there is a structural component present and > that the local trend varies from one location to another. Plotting a > semivariogram for each input sample point is not practical. > > Universal Kriging assumes the model: Z(s) = µ(s) + ε(s), where > µ(s) is some deterministic function. <Imagine a figure here> > For example, in the following figure, which has the same data that was used > for Ordinary Kriging concepts, the observed data is given by the solid > circles. A second-order polynomial is the trendâlong dashed lineâwhich > is µ(s). If you subtract the second-order polynomial from the original > data, you obtain the errors, ε(s), which are assumed to be random. The mean > of all ε(s) is 0. Conceptually, the autocorrelation is now modeled from the > random errors ε(s). Of course, you could have fit a linear trend, a cubic > polynomial, or any number of other functions. The figure above looks just > like a polynomial regression from any basic statistics course. In fact, that > is what Universal Kriging is. You are doing regression with the spatial > coordinates as the explanatory variables. However, instead of assuming the > errors ε(s) are independent, you model them to be autocorrelated. The advice > is the same as Ordinary Kriging: there is no way to decide, based on the data > alone, on the proper decomposition. > >> >> 3. Are there any tricks in gstat that make interpolation with gstat >> especially fast or slow? For example, if the numbers were treated faster if >> they were between 0 and 1, or if NA should be recoded, or if I should never >> use a data frame but always a matrix. That sort of thing. >> > Not the things you mention; defining a local neighbourhood matters -- > small neighbourhoods are faster, but at some stage also produce worse > interpolations. > > For filling in missing strips in images I'd try to use segmented > neighbourhoods, to guarantee values are taken from all sides. Quadrant > search is available in the gstat C code, but not interfaced through the > R package. I will implement it if you volunteer to test it. Thanks! I will definitely volunteer to test that if it ends up being necessary or a big advantage. Let's see how gstat works out as-is and I will let you know. > >> If not, I'll resort to trial and error and diving deeper into the gstat >> documentation. But it seems like a good idea to check in first. For anyone >> interested, I am trying to repair Landsat satellite photos, which have big >> broken strips of nodata due to a mechanical failure. The strips are up to >> 14 pixels wide. I need to do this 2000x2000 interpolation about 2000 >> times-- so speed differences of even a few seconds or minutes are quite >> important.. > > I'm interested in your findings! > > BTW, just out of curiosity, why don't you do the whole thing with ArcGIS? > For two reasons, ArcGIS is not the best for this. First, there is a lot of other R code that happens before this interpolation. So if I can write something that only requires one program, that would be much cleaner. The second is that ArcGIS is expensive, and so if it can all be freely distributed and accessible, that will be better. A third reason is that I produce a lot of these images, in parallel, and so I need multiple processes making these little images at once. Thanks for helping me along; i am a little over my head here- not with the programming, but in understanding the tiny details of the geostats that make a big difference. Thanks ! >> >> Thanks! >> Jeff >> >> ------------------------------------------ >> Prof. Jeffrey Cardille >> jeffrey.cardi...@umontreal.ca >> >> ****************************************************************************************** >> ** http://www.geog.umontreal.ca : Web >> ** >> ** http://sites.google.com/site/jeffcardille : Plans de Cours et >> Recherche ** >> ** >> ** >> ** "Petit à petit l'oiseau fait son nid ..." N'hésitez pas à corriger >> mon français ** >> ********************************************************************************************************* >> ** Département de Géographie ** Bureau: >> ** >> ** assistant professor / professeur adjoint ** Salle 440 >> ** >> ** Université de Montréal ** Pavillon >> Strathcona ** >> ** C.P. 6128 ** 520, chemin >> de la Côte-Ste-Catherine ** >> ** Succursale Centre-ville ** Montreal, QC >> H2V 2B8 ** >> ** Montréal, QC, H3C 3J7 ** Télé: >> (514) 343-8003 ** >> ********************************************************************************************************* >> ** >> ** Pour aider à Haiti avec votre intérêt de géographie, visitez: >> ** http://groups.google.com/group/udem-geo-aide/ >> ** >> ********************************************************************************************************* >> >> >> >> >> >> [[alternative HTML version deleted]] >> >> >> >> ------------------------------------------------------------------------ >> >> _______________________________________________ >> R-sig-Geo mailing list >> R-sig-Geo@stat.math.ethz.ch >> https://stat.ethz.ch/mailman/listinfo/r-sig-geo > > -- > 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.52north.org/geostatistics e.pebe...@wwu.de Jeff ------------------------------------------ Prof. Jeffrey Cardille jeffrey.cardi...@umontreal.ca ****************************************************************************************** ** http://www.geog.umontreal.ca : Web ** ** http://sites.google.com/site/jeffcardille : Plans de Cours et Recherche ** ** ** ** "Petit à petit l'oiseau fait son nid ..." N'hésitez pas à corriger mon français ** ********************************************************************************************************* ** Département de Géographie ** Bureau: ** ** assistant professor / professeur adjoint ** Salle 440 ** ** Université de Montréal ** Pavillon Strathcona ** ** C.P. 6128 ** 520, chemin de la Côte-Ste-Catherine ** ** Succursale Centre-ville ** Montreal, QC H2V 2B8 ** ** Montréal, QC, H3C 3J7 ** Télé: (514) 343-8003 ** ********************************************************************************************************* ** ** Pour aider à Haiti avec votre intérêt de géographie, visitez: ** http://groups.google.com/group/udem-geo-aide/ ** ********************************************************************************************************* [[alternative HTML version deleted]]
_______________________________________________ R-sig-Geo mailing list R-sig-Geo@stat.math.ethz.ch https://stat.ethz.ch/mailman/listinfo/r-sig-geo
