[R-sig-Geo] Projected data in R-Gstat

[Edzer Pebesma]

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]]
> 
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at r-project.org
> 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.pebesma at wwu.de