gstat variogram & great circle distance

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")

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")

foo

coordinates(foo) <- ~lon+lat
proj4string(foo)=CRS("+longlat")
variogram(z~1, foo, cloud=T, cutoff=10000)

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")