[R-sig-Geo] georeferencing point shape

[Roger Bivand]

On Tue, 15 Jun 2010, Tom Gottfried wrote:

> Hi Tim,
>
>> I have a point shape containing field measurements in a small scale studyarea. Unfortunately these
> points didn`t have geographic or projected coordinates. The spatial reference for these points are
> theodolite measurements. For four corner points gps measurements of projected coordinates are
> available. I want to georeference my point shape in these projected coordinate system and have no
> idea how.
>> Here is a subset of my data:
>>
>> x <- c(5.786, -3.766, -12.613, -21.836, -26.340, 3.958, -11.120, -17.266, -18.938,
>>   16.751, -21.507, 26.420, -19.916, 23.184, 9.660, -5.711, -18.256, 27.888, 12.634, -33.510)
>> y <- c(4.470, 0.797, -3.130, -7.656, -9.313, 8.700, 6.923, 8.545, 12.468, 4.748,
>>   -20.920, -5.396, -24.830, -11.636, 22.462, 19.210, -28.851, -9.510, 27.421, 8.035)
>> z <- c(1100.493, 1100.867, 1101.798, 1102.559, 1103.703, 1102.366, 1107.249, 1110.781,
>>  1113.920, 1096.967, 1095.284, 1088.956, 1092.869, 1086.786, 1101.300, 1110.197, NA, NA, NA, NA)
>> gps_x <- c(NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA, 4560136, 4560184, 4560172, 4560121)
>> gps_y <- c(NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA, 5289465, 5289482, 5289521, 5289502)
>>
>> d <- data.frame(x,y,z,gps_x,gps_y)
>> d
>> coordinates(d) = ~x+y
>> plot(d, axes=T)
>>
>> gps_points <- d[17:20,]
>> points(gps_points, pch=19, col="red")
>>
>> Because both measurements - the theodolite values and the gps values - are measured in meters, my
>> first idea was to calculate the coordinates for all points by sum the difference between a gps
>> point and every single theodolite point in x- and y-direction. Unfortunately there are little
>> differences depending which gps point is used as basis.
>
> If you can assume the little differences being due to the error in gps-measurements, you could
> calculate your coordinates with all your gps points as a basis and then take the mean as final value.
>
> Tom

There is a problem with the known points too, in that their surveyed and 
GPS interpoint distances differ somewhat:

dist(cbind(gps_x, gps_y)[17:20,])
dist(cbind(x, y)[17:20,])
dist(cbind(gps_x, gps_y)[17:20,]) - dist(cbind(x, y)[17:20,])

This means that with only 4 GPS points, it is hard to model the 
transformation with more terms:

df <- data.frame(gps_x, gps_y, x, y)
lmx <- lm(gps_x ~ x + y, data=df[17:20,])
lmy <- lm(gps_y ~ x + y, data=df[17:20,])
nx <- predict(lmx, data.frame(x, y))
ny <- predict(lmy, data.frame(x, y))
dist(cbind(nx, ny)[17:20,])
dxy <- dist(cbind(x, y))
dnxny <- dist(cbind(nx, ny))
all.equal(dxy, dnxny, check.attributes=FALSE)

More GPS measurements would have helped, but at this scale, GPS are always 
going to be approximate. You cannot measure the centres of the GPS and 
surveyed points accurately either, I'm afraid. To fix things, you would 
need GPS measurements in the precision of the data set, so in most cases 
like this DGPS rather than GPS. If one of the surveyed points is clearly 
visible on a registered high resolution image, you could use that instead 
of the GPS - something like the corner of a building (measured to about 
1cm and known in the projection of interest). However, planning this 
before doing the fieldwork was the only effective remedy.

Roger

PS. Do you know that the GPS was using the ellipse you specify? Which 
datum are you assuming (WGS84)?


>
>> Is there a better way to transform the point pattern into the projected 
>> coordinate system? The proj4string for my coordinate system is: 
>> +proj=tmerc +lat_0=0 +lon_0=12 +k=1 +x_0=4500000 +y_0=0 +ellps=bessel 
>> +units=m +no_defs
>>
>> Thank you very much.
>>
>> TIM
>>
>>
>
>

-- 
Roger Bivand
Economic Geography Section, Department of Economics, Norwegian School of
Economics and Business Administration, Helleveien 30, N-5045 Bergen,
Norway. voice: +47 55 95 93 55; fax +47 55 95 95 43
e-mail: Roger.Bivand at nhh.no