[PROJ] Fwd: Re: Help with impossibly elusive horizontal datum

[molnar]

Dear Landon,

This is not your fault. There is no simple (7-parameter) definition 
between the applied AMS-Haiti datum and any modern datum.
We have to admit, that our predecessors not always
knew, what they did (or more probably, they did not have enough time, to 
finish their work properly).
Fortunatly, other predecessors invented the correction grid method, and 
this seems to be working.

I really wanted to analyze, what might be the problem with the Haiti 
maps.
I downloaded some 25 sheets of the 1st editions of the western part of 
Haiti from the website.
I georefenced them in QGIS in their native map projection, by clicking 5 
grid intersection points per sheet, to have an affine transformation 
(plus quality control) between image coordinates and printed utm 
coordinates.
Next I restarted geocoding by collecting 15-20 control points (GCPs) per 
sheet, so I clicked the same surface object on raw images (not using the 
previous transformation) and clicked the same point on google satellite 
map.
Using this latter dataset, for some hundred GCPs, I generated 
Clrk66-WGS84 coordinate pairs.
This is, how I got the Clrk66 coordinets of the GCPs: I converted the 
GCPs image coordinates, with affine transformation to overprinted UTM 
coordinates, and with inverse proj transformation I transformed them to 
Clrk66 ellipsoidal coordinate)
This is how I got WGS854 coordinates: From the webmercator coordinates 
of the GCPs on google satellite, with inverse proj  I got WGS84 
ellipsoidal coordinates.
(Certainly the GCPs obtained with this method might have as much as 30 
meter errors, but it is half a millimetre on map, and Landon's problem 
is the kilometer shift, that was not explained.)

Analyzing these coordinate pairs some conclusions:
-I could not find a single 7 parameter BursaWolf, (towgs84) 
transformation,
- For a single map sheet, finding and clicking GCPs on a map sheet, 
usually the RMS error of a 2D Helmert transformation between raw image 
coordinates and the google earth webmercator coordinates was better than 
30 meter, that means the original mapping material used for this map was 
quite ok. (Any map projection on a 1:50000 map on a 25x18km map sheet 
works like this, if the mapping material is ok)
- I generated a coordinate difference chart, (the attached 
"latlon_diff_vectors.jpg") the coordinate differences are plotted as 
vectors. The vector field seems to be continuous, but not uniform. On 
top left corner of the area, the for Dame Marie map sheet, the 
differences are
negligible, but for other map sheets (like the Jerome sheet, next to it) 
the ellipsoidal coordinate differences represent several hundred meter 
differences. This implies, that even with correction grid, at the 
vicinity of the borders of Dame Marie map sheet, there might be errors 
as much as 200 meters.

Using the ellipsoidal coordinate difference dataset, I generated a GSB 
(Grid Shift Binary) for horizontal grid shift with some GDAL commands. 
(Say thank to GDAL developers.) It took me a while, (This might be 
another post :-)), but finally I got it. After getting through some QGIS 
version problem (where to store the GSB files...) it worked as expected.

And finally:
For using the coordinate correction grid (attachment to this email), in 
QGIS 3.24, this is what to do:
1. Copy the attached "haiti_ams_wgs84.gsb" to
C:\Users\USER\AppData\Roaming\QGIS\QGIS3\profiles\default\proj (USER is 
your username) This is for Windows, for any other OS, check Qgis manual.
2. Define a user defined map projection in QGIS, using the deprecated, 
not recommended proj definition:
"+proj=utm +zone=18 +ellps=clrk66 +nadgrids=haiti_ams_wgs84.gsb +units=m 
+no_defs"
3. Geocode your map sheets by some utm grid intersections, choosing the 
"To" coordinate system the user-defined coordinate system.
4. You can add the geocoded sheet to any project (that has 
transformation to WGS84), QGIS will warp your mapsheets on-the-fly, so 
it might be slow.
5. If it is too slow, you can resample your mapsheet into another 
coordinate system with gdalwarp, and use the warped dataset in QGIS. For 
QGIS it is less CPU consuming to handle a dataset that has a 7-parameter 
transformation definition to present day map datums (instead of grid 
shift definition) .

It should work for western Haiti. If you need a higher accuracy grid, or 
the whole country, you have to provide me GCPs, and I can generate the 
GSB file.
Best wishes,
Gabor

The mentioned files:
latlon_diff_vectors.jpg:
https://pasteboard.co/rEwS1oORV8i7.jpg

I could not upload "haiti_ams_wgs84.gsb" on the pasteboard.co site, 
neither in tif format, so I send it to you in a private email...




> Landon Yarrington <lyarrington at gmail.com> writes:
> 
>> First, the three control extension points listed in the map info are 
>> quite
>> far away and not on the map. And, the coordinates given for the points 
>> do
>> not appear to align with the actual locations when viewed on Google 
>> Earth.
>> Maybe this is a clue to the offset?
> 
> This is exactly a clue.  The lat/lon given in the map data are labeled
> astronomic, and they are the results of some number of measurements
> which have considerable error.  They decided to fix those stations with
> those values; this is the definition of a datum.  In the 1940s I am
> aware of no better way to connect coordinates in Haiti to coordinates 
> in
> CONUS.  Surely the Army Map Service was world class in doing things 
> like
> that.
> 
>> Fort Islet lighthouse 18°33'31.33" -72°20'59.03", should be 
>> 18°33'18.08" > -72°21'04.50"
>> Santo Domingo lighthouse 18°27'53.64" -69°52'52.74", should be 
>> 18°27'51.39" -69°52'34.90"
>> Dame Marie 18°36'47" -74°25'53", unknown but possibly 18°36'57.21" 
>> -74°25'26.94"
> 
> You say "should be", but you are really saying "In some particular
> modern coordinate system".  The definition of 0 longitude even in ITRF
> is totally arbitrary.  It is just that the values adopted for mapping
> Haiti and the values measured in a modern frame are different.  Yes, 
> for
> latitude, one can say they are wrong, but gravity variation makes this
> hard.
> 
> Looking only at minutes, and local - wgs84
> 
> 		lat		lon
> Fort Islet	13"		5.5"
> Santo Domingo	2.2"		-18"
> Dame Marie	-10"		-26"
> 
> so this is not really that consistent.  This leads me to think there is
> rotation too.
> 
> But really, you have a map which you hope is internally consistent and
> you can identify points on the modern ground that match points on the
> map.
> 
>> I took your advice and used the QGIS georeferencer. This is actually 
>> what I
>> initially did before reaching out here. I put google and bing satelite
>> imagery on the canvas to pin ground control points. I chose 19 GCPs 
>> and set
> 
> There is the question of whether those imagery layers are aligned to
> recent WGS84, recent NAD83, something else, or are just off.  But
> probably they are decent.
> 
> I have checked MassGIS imagery against NAD83(2011) epoch 2010.0 (which
> the imagery is said to be in) and it matches to a pixel.  But it was
> taken recently, in an area where there is great ground control and a
> statewide RTK network.
> 
>> the algo to thin plate spline. The result has excellent accuracy. 
>> However,
>> the result is not easily reproducible. I have personally been to all 
>> the
>> locations I chose as GCPs and I know that---for the purpose of
>> georeferencing---those sites are unchanged in any meaningful way since 
>> the
>> 1942 aerial photography used to make the map was taken.
> 
> It's great you have enough local knowledge to know what's undisturbed,
> and this is basically necessary to do what you are doing correctly.
> 
>> Not everyone is going to have this specialized knowledge, and why
>> should they when there is already a UTM grid on the map? Can the grid
>> intersections make the georeferencing process more reproducible
>> without the need of ethnographic fieldwork?
> 
> I think you are fundamentally misunderstanding the nature of UTM.  UTM
> is a projection, not a datum.  It takes latitude and longitude -- in
> some datum, relative to something -- and turns it into easting and
> northing.  This is just math, and has no bearing on or input from
> whether the latitude/longitude pairs are consistent with modern WGS84.
> 
> Further, this is NAD27-style UTM with the Clarke ellipsoid, not
> NAD83/WGS84-style UTM with the GRS80 ellipsoid.
> 
> Assuming that the UTM coordinates are in NAD27 is exactly as wrong an
> assumption as assuming that the map lat/lon are in NAD27.
> 
>> The map in question is the 1st edition (
>> https://www.loc.gov/resource/g4940m.gct00206/?sp=168). There are at 
>> least 4
>> subsequent editions (here is the 3rd
>> https://www.loc.gov/resource/g4940m.gct00206/?sp=167). The 1st edition 
>> has
>> the map info I quoted in my first message but is dropped from all 
>> later
>> editions. All later editions explicitly state the datum is NAD27 (or 
>> NAD83
>> for 4th and 5th editions).
> 
> The link you sent is the 2nd, but that's better, as it is the first 
> that
> says NAD27, and the compilation date is 1963.
> 
> I am not aware of any technology that would enable extending NAD27 from
> CONUS to Haiti in 1963.  That seems early for TRANSIT, but the dates in
> wikipedia do not 100% preclude AMS trying to tie in Haiti with it:
>    https://en.wikipedia.org/wiki/Transit_(satellite)
> 
> If anyone knows any other techniques for tying Florida to Haiti in 
> 1963,
> please post.  I know this isn't really Surveying-L but I bet a lot of
> readers would be interested.
> 
>> Curiously, the geographic representation in the
>> 1st edition is more expansive than the others; even though all the 
>> editions
>> have their neatlines within 19°50', -72°15' and 19°40', -72°00', the 
>> 1st
>> edition shows considerably more area. Example: the town of Limonade at 
>> the
>> bottom center of the map. The road that leads southeast from Limonade
>> extends about 1.29 km further than what is shown in later editions. In
>> fact, if you look at the grid intersection 2178000m N and 801000m E on 
>> the
>> 1st edition and the same point on the later editions, the grid is 
>> nearly
>> 1000 meters off. Compare
>> https://www.loc.gov/resource/g4940m.gct00206/?sp=168&r=0.314,0.573,0.372,0.229,0
>> to
>> https://www.loc.gov/resource/g4940m.gct00206/?sp=167&r=0.284,0.495,0.238,0.147,0.
>> What is going on here?? Is this a datum difference, projection 
>> difference,
>> both??
> 
> If the neatline are the same lat/lon and they show different content,
> then either they have fixed errors or the datum is different.  Same 
> with
> projected coordinates.  But the datum is different; they say so.
> 
> I would get the NAD27 coordinates for their 3 control points and then
> compute distance/bearing from "AMS Haiti Datum of 1948" (I'm making 
> that
> up) to NAD27.  Those aren't real distances, but a datum shift.
> 
>> So... the second thing I tried was selecting GCPs at each grid 
>> intersection
>> on the 1st edition but subtracting 1000 meters from the northings. (I 
>> did
>> 141 GCPs.) This got the map's accuracy much better, but there was 
>> still a
>> systematic offset.
> 
> This doesn't really make sense.  You are assuming that the original
> coordinate system is the same as some modern one and it just
> isn't. Think about if you were given a map in NAD27 and you are trying
> to georeference it into WGS84, but nobody will given you a transform
> between the two datums.
> 
> Or how about you had an army of surveyors to do an experiment.  Go to
> modern place that is 20km by 20 km.  Pick 3 places and drive pins.  For
> each, get really good GPS-derived coordinates.  Then, generate 6 random
> numbers and add them into the 3 pairs of coordinates.  Write those down
> and hand them to the surveying teams.  The teams get to use total
> stations, with 1 arc second theodolites, and 1-2 ppm electronic 
> distance
> measurement, but are forbidden from using any satellites, any outside
> info, and from making any astronomic observations.  They are to survey
> the town.
> 
> The map they make will be highly consistent with itself.  But there 
> will
> be some offset and rotation because the assumed coordinates of the 3
> control points are in error.
> 
> This is, as I understand it, essentially what happened, except that
> there was no oracle with random numbers.  There was instead the errors
> in astronomic determination of latitude and longitude.  Besides
> measuremnt errors, the direction of gravity is not aligned with the
> normal to the ellipsoid, and this results in latitude errors.
> 
>> The best accuracy was in the center of the map. In the
>> 4th edition of the map, there is a note instructing how to convert 
>> from
>> NAD27 to NAD83, which reads, "Grid: Add 35mE, add 198mN. Geographic:
>> Subtract 1.4" Long, Add 2.2" Lat." OK. So, I also did these shifts to 
>> the
>> GCPs that I had already offset by 1000m N. For each easting, I added 
>> 35m
>> and each northing I added 198m. Then, I ran gdalwarp (code below). The
>> result still has an offset but it's the best I've been able to get 
>> using
>> the UTM grid, but sadly it is not nearly as accurate as the QGIS 
>> visual
>> approach.
> 
> You are just fudging numbers by treating the UTM grid as true, and
> basically computing a single E/N shift by iteratively adjusting and
> looking.  But you are not estimating a rotation.  When you use GCPs 
> with
> map coordinates in pixels and accurate WGS84 from imagery (if that is
> indeed acccurate), the georeferencer is doing least squares among all 
> the
> points.
> 
>> Again, I am indebted for any help. The code below only uses 8 GCPs of 
>> the
>> UTM intersections offset and shifted described above.
>> 
>> ```
>> # download the tiff from Library of Congress
>> wget
>> https://www.loc.gov/resource/g4940m.gct00206/?download=https%3A%2F%2Ftile.loc.gov%2Fstorage-services%2Fmaster%2Fgmd%2Fgmd4m%2Fg4940m%2Fg4940m%2Fgct00206%2Fcs000168.tif
>> -O cs000168.tif
>> # add 8 ground control points using the UTM grid described above
>> gdal_translate -of GTiff -gcp 833.408 5291.81 789035 2179198 -gcp 
>> 1146.904
>> 5296.971 790035 2179198 -gcp 3963.057 6290.092 799035 2176198 -gcp 
>> 843.374
>> 4660.743 789035 2181198 -gcp 1763.341 5935.218 792035 2177198 -gcp 
>> 2077.227
>> 5941.007 793035 2177198 -gcp 8369.034 5735.588 813035 2178198 -gcp 
>> 8684.098
>> 5741.097 814035 2178198 cs000168.tif cs000168_GCPs.tif
>> # warp the map with thin plate spline
>> gdalwarp -r near -tps -co COMPRESS=JPEG -co TILED=YES -s_srs 
>> EPSG:32618
>> -t_srs EPSG:32618 cs000168_GCPs.tif cs000168_TO_epsg32618.tif
>> ```
> 
> Why are you running gdal_translate rather than just using the
> georeferencer and letting it make a new tif?
> 
> Another consideration is that the map datum is a local one, but
> presumably distances are correct.   So there should only be one
> translation and one rotation to go from map datum to NAD27.   Or, one
> translation, one rotation, and one scale factor to go from map datum 
> UTM
> to WGS84 UTM.  And, the scale factor can be computed since the
> translation is small.    By using a complicated mapping fit, you are
> allowing distortions in the map datum to be leveled out.   If the map 
> is
> distorted locally because of survey errors, that's good.  But if the 
> map
> is internally consistent and the individual points have random errors,
> that's bad.
> _______________________________________________
> PROJ mailing list
> PROJ at lists.osgeo.org
> https://lists.osgeo.org/mailman/listinfo/proj