Projected data in R-Gstat

Dear R-Sig-Geo,


I am trying to interpolate air pollution from air monitoring thanks to
ordinary kriging in order to assign the pollution exposure for a sample of
individual for the city of Phoenix, Arizona, using the Gstat package.


However, I have some trouble to understand how my data should be projected
in order to fit the variogram. My data are in longitude/latitude:


id    co1h    latitude    longitude
1    390    33.49462    -112.1310
2    470    33.48385    -112.1426
3    170    33.41045    -111.8651
4    210    33.56033    -112.0663
5    210    33.56936    -112.1915
6    360    33.45793    -112.0460
7    200    33.47968    -111.9172
8    300    33.46093    -112.1175
9    370    33.40316    -112.0753
10    180    33.29898    -111.8843
11    240    33.41240    -111.9347
12    150    33.63713    -112.3418
13    70    33.37005    -112.6207
14    310    33.50318    -112.0956


In order to fit the variogram I use the following commands with:

 co<-read.dta("co1h.dta")

coordinates(co) = ~longitude+latitude

v <- variogram(log(co1h) ~ 1, co)

plot(v)

vgm<-vgm(0.59638453, "Hol", 0.1490525,0.02910648, anis = c(90, 0.8))

v.fit<-fit.variogram(v, vgm)

v.fit

plot(v,fit.variogram(v, vgm))



Does anyone know if my projection is correct or whether I need to change
it? If it is necessary to change, do you know how can I do it and what
should be the correct projection?

Thank you very much for your help!

Kind regards,

Marguerit David

PhD in economics,

University Paris Dauphine, France

	[[alternative HTML version deleted]]

On 05/11/2012 12:35 PM, David Marguerit wrote:

Dear R-Sig-Geo,


I am trying to interpolate air pollution from air monitoring thanks to
ordinary kriging in order to assign the pollution exposure for a sample of
individual for the city of Phoenix, Arizona, using the Gstat package.


However, I have some trouble to understand how my data should be projected
in order to fit the variogram. My data are in longitude/latitude:


id    co1h    latitude    longitude
1    390    33.49462    -112.1310
2    470    33.48385    -112.1426
3    170    33.41045    -111.8651
4    210    33.56033    -112.0663
5    210    33.56936    -112.1915
6    360    33.45793    -112.0460
7    200    33.47968    -111.9172
8    300    33.46093    -112.1175
9    370    33.40316    -112.0753
10    180    33.29898    -111.8843
11    240    33.41240    -111.9347
12    150    33.63713    -112.3418
13    70    33.37005    -112.6207
14    310    33.50318    -112.0956


In order to fit the variogram I use the following commands with:

 co<-read.dta("co1h.dta")

coordinates(co) = ~longitude+latitude

v <- variogram(log(co1h) ~ 1, co)

plot(v)

vgm<-vgm(0.59638453, "Hol", 0.1490525,0.02910648, anis = c(90, 0.8))

v.fit<-fit.variogram(v, vgm)

v.fit

plot(v,fit.variogram(v, vgm))



Does anyone know if my projection is correct or whether I need to change
it? If it is necessary to change, do you know how can I do it and what
should be the correct projection?

Yes, it is not correct. You didn't set it, so

proj4string(co)

will give you NA (missing). This means that gstat will assume they are
metric (in m, or km), and compute distances using regular (Euclidian)
distances -- which is wrong!

You have to do the following:

proj4string(co) =  "+proj=longlat +datum=WGS84"

(assuming this is the correct datum setting!)

and then rerun the variogram fitting. In this case, gstat will use great
circle distances instead of Euclidian distances.

An alternative would be to reproject your data, e.g. to the appropriate
UTM zone, and then re-run the variogram estimation and fitting.

http://en.wikipedia.org/wiki/Universal_Transverse_Mercator_coordinate_system

might point you to the appropriate zone.

With best regards,

Thank you very much for your help!

Kind regards,

Marguerit David

PhD in economics,

University Paris Dauphine, France

	[[alternative HTML version deleted]]