[PROJ] Help with impossibly elusive horizontal datum

[Greg Troxel]

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.