[R-sig-Geo] gstat variogram & great circle distance
Timothy W. Hilton
hilton at meteo.psu.edu
Wed Apr 30 17:44:45 CEST 2008
Hi Edzer,
I am still having some trouble with the great distance calculations in 'variogram'. Your suggestion below works, but the distances are not correct (at least, not in kilometers, meters, or miles). I do not have proj.4 or gdal libraries installed, nor do I have the R packages proj4 or rgdal. It seems like I should not need them, as I am not doing a projection. Given that I am setting the 'proj4string' attribute in order to achieve a great cicle distance calculation, though, I wonder if not having those packages installed is a problem.
An example (dump output for foo is below):
> foo
z lon lat
1 -8.2920352 -68.74028 45.20407
2 0.3962574 -157.40894 70.46961
3 -5.3976371 -89.97919 46.08268
> coordinates(foo) <- ~lon+lat
> proj4string(foo)=CRS("+longlat")
> variogram(z~1, foo, cloud=T, cutoff=10000)
dist gamma dir.hor dir.ver id left right
1 9880.046 37.74321 0 0 var1 2 1
2 2366.200 4.18877 0 0 var1 3 1
The distance from 1 to 2 is 5299 km, and from 3 to 1 is 1650 km. I'm not sure what the 9880.046 and 2366.200 represent. The ratio of the variogram$dist values to the correct distances in km are not the same, so those values cannot both be correct distances in any units.
Perhaps I am missing a package that gstat needs?
Many thanks for your help.
Cheers,
Tim
=======
foo <-
structure(list(z = c(-8.29203519866722, 0.396257381218808, -5.39763713302683), lon = c(-68.740278, -157.408944, -89.97919), lat = c(45.20407,70.469611, 46.08268)), .Names = c("z", "lon", "lat"), class = "data.frame", row.names = c(NA,3L))
=======
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
> >
>
More information about the R-sig-Geo
mailing list