georeferencing point shape
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