[R-sig-Geo] georeferencing point shape

[Tom Gottfried]

Am 15.06.2010 15:48, schrieb 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,])

Yes! And the surveyed interpoint distances should be more accurate. Thus, I think it's better to
rely on the surveyed (relative) coordinates, rather than transforming them to obtain the best fit
with the erroneous GPS-measurements.
But "snapping" the relative coordinates to the projected coordinate system by something like

x_new <- rep(0, times=length(x))
y_new <- x_new
for (i in 17:20){x_new <- x_new + gps_x[i]+x-x[i]; y_new <- y_new + gps_y[i]+y-y[i]}
x_new <- x_new/4; y_new <- y_new/4

(as I suggested last time), assumes that the error of GPS is random in both dimensions and has an
expected value of 0. Might this assumption be unrealistic?

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

Orthophotos (I think with a resolution of 40 cm) are available. This should be already much better
than the mean of the GPS measurements.

Tom

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

-- 
Technische Universität München
Department für Pflanzenwissenschaften
Lehrstuhl für Grünlandlehre
Am Hochanger 1
85350 Freising / Germany
Phone: ++49 (0)8161 715324
Fax:   ++49 (0)8161 713243
email: tom.gottfried at wzw.tum.de
http://www.wzw.tum.de/gruenland