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

[Paul Hiemstra]

On 05/14/2012 09:51 PM, Edzer Pebesma wrote:
>
> 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.

...in addition, for reprojection of data you can use the spTransform
function from the rgdal package.

regards,
Paul

> 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


-- 
Paul Hiemstra, Ph.D.
Global Climate Division
Royal Netherlands Meteorological Institute (KNMI)
Wilhelminalaan 10 | 3732 GK | De Bilt | Kamer B 3.22
P.O. Box 201 | 3730 AE | De Bilt
tel: +31 30 2206 494

http://nl.linkedin.com/pub/paul-hiemstra/20/30b/770