[R-sig-Geo] gstat variogram & great circle distance

[Timothy W. Hilton]

This works great, Edzer.  I expected there was a simple solution.  Many thanks!

-Tim

On Wed, Apr 2008, 30 at 07:42:36AM +0200, Edzer Pebesma wrote:
> Timothy, for some reason the projected argument was not meant to be set 
> by users at this level of abstraction; I'll look into it. The following 
> seems to work:
> 
> > proj4string(foo)=CRS("+longlat")
> > proj4string(foo)
> [1] "+longlat"
> > variogram(z~1,foo)
>  np      dist        gamma dir.hor dir.ver   id
> 1  1  177.7815   0.48861387       0       0 var1
> 2  2  614.0040   1.11639574       0       0 var1
> 3  3  715.8402 115.50300578       0       0 var1
> 4  1  829.6047   0.01873678       0       0 var1
> 5  1  893.0483   3.66061567       0       0 var1
> 6  1  992.9436 268.46555052       0       0 var1
> 7  4 1095.7359  83.25299050       0       0 var1
> 8  1 1274.4497  12.33954872       0       0 var1
> 9  1 1357.4796   3.01500925       0       0 var1
> 
> --
> Edzer
> 
> 
> Timothy W. Hilton wrote:
> >Hello,
> >
> >I am trying to use gstat to compute a semivariogram for data whose 
> >coordinates are latitude/longitude pairs.  I would like to use the great 
> >circle distance between pairs.  The documentation implies that gstat can 
> >do this, but I am not having any success.  If anyone could suggest the 
> >correct syntax, I would greatly appreciate it.
> >
> >Here is a sample of my data (see output from dump below):
> >
> >
> >  
> >>foo
> >>    
> >            z        lon      lat
> >1  -1.9582483 -125.29228 49.87217
> >2  -1.9158902  -82.15560 48.21670
> >3   4.2221176  -98.52472 55.90583
> >4   3.2335693  -99.94833 56.63583
> >5   1.1203839 -104.69174 53.91626
> >6   0.3461385  -79.42083 39.06333
> >7   1.1258993 -105.10053 48.30788
> >8  23.5179123  -88.29187 40.00610
> >9   3.0519159  -72.17148 42.53776
> >10  3.2026143 -121.55694 44.44889
> >11 -2.1094711  -89.34765 46.24202
> >
> >I can calculate a variogram:
> >
> >  
> >>coordinates(foo) <- ~lon+lat
> >>variogram(z~1, foo)
> >>    
> >   np      dist       gamma dir.hor dir.ver   id
> >1   1  1.599865   0.4886139       0       0 var1
> >2   2  5.545490   1.1163957       0       0 var1
> >3   4  6.712018  86.6319381       0       0 var1
> >4   1  8.038953   3.6606156       0       0 var1
> >5   3  9.422337  91.0816908       0       0 var1
> >6   2 10.149322 164.1162183       0       0 var1
> >7   2 11.868366   7.6772788       0       0 var1
> >8   1 13.326965  20.0445076       0       0 var1
> >9   1 14.846073  14.2740402       0       0 var1
> >10  1 15.887767   5.2338108       0       0 var1
> >11  3 16.792331  72.2669527       0       0 var1
> >12  2 17.828085  16.0787636       0       0 var1
> >
> >The distances are clearly not great circle distances, though.  Setting the 
> >"projected" flag to "false" gives me this error:
> >
> >  
> >>variogram(z~1, foo, projected=FALSE)
> >>    
> >Error in variogram.default(y, locations, X, trend.beta = beta, grid = 
> >grid,  :
> >  formal argument "projected" matched by multiple actual arguments
> >
> >Thanks in advance for any help,
> >Tim
> >
> >==================
> >`foo` <-
> >structure(list(z = c(-1.95824831109744, -1.91589016435630, 
> >4.22211761150161,
> >3.23356929459598, 1.12038389231868, 0.34613850821113, 1.12589932643631,
> >23.5179122516170, 3.05191586902680, 3.20261431141517, -2.10947106854739
> >), lon = c(-125.29228, -82.1556, -98.524722, -99.948333, -104.691741,
> >-79.420833, -105.100533, -88.291867, -72.171478, -121.556944,
> >-89.34765), lat = c(49.87217, 48.2167, 55.905833, 56.635833,
> >53.916264, 39.063333, 48.307883, 40.0061, 42.537756, 44.448889,
> >46.242017)), .Names = c("z", "lon", "lat"), row.names = c(NA,
> >-11L), class = "data.frame")
> >
> >_______________________________________________
> >R-sig-Geo mailing list
> >R-sig-Geo at stat.math.ethz.ch
> >https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> >  
>