[PROJ] Help with impossibly elusive horizontal datum
Landon Yarrington
lyarrington at gmail.com
Fri Jul 12 13:55:50 PDT 2024
Thank you very much for your response!! Please note I've included minimal
reproducible code at the end of the email, and I would be very grateful to
anyone who would care to take a look and help.
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?
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"
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
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. 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?
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). 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??
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. 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.
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
```
On Tue, Jul 2, 2024 at 4:40 AM Greg Troxel via PROJ <proj at lists.osgeo.org>
wrote:
> Landon Yarrington via PROJ <proj at lists.osgeo.org> writes:
>
> > I am trying to georeference an US Army map of northern Haiti made from
> > aerial photographs taken in 1942 and 1944 using `gdal_translate` and
> > `gdalwarp`. I need a sanity check because this is impossible, and so I'm
> > asking for help. Here is the map from Lib. of Congress (
> >
> https://www.loc.gov/resource/g4940m.gct00206/?sp=160&r=0.365,0.623,0.281,0.173,0
> > ).
> >
> > The map info states,
>
> > Horizontal Datum is based on the following astronomic values:
> > Fort Islet Lighthouse Lat. 18°33'31.33", Long. 72°20'59.03"
> > Santo Domingo Lighthouse Lat. 18°27'53.64", Long. 69°52'52.74"
> > Cape Dame Marie Astro Lat. 18°36'47", Long. 74°25'53"
> > Extensions of control from Astronomic Positions have been made by
> multiplex
> > triangulation.
>
> This is remarkably clear and useful even if it is hard to understand.
>
> The mapmaker is saying that they made astronomic determinations of
> position at 3 named locations, presumably obvious on the map, and gives
> the values. These coordinates are not related to any particular datum.
> Instead, they define a local datum.
>
> Critically, there is absolutely no reason to expect that the coordinates
> above match WGS84, other than that the time signals are consistent with
> greenwhich, and latitude is physical. But that doesn't get you the last
> 10s to 100s of meters.
>
> Almost certainly they had the benefit of radio time signals, so their
> longitudes should be pretty good. Latitudes should be pretty good
> without needing to think about time. Without really thinking, 1000m for
> a single determination strikes me as pretty good. Before NAD83, which
> used satellite/etc. data, NAD27 was different from modern ITRF by 40m
> around me. And that's a national datum that got to average multiple
> observations, even if it was held to be somewhat consistent with
> mid/late 1800s work.
>
> Our current datums derive from averages of many astronomic observations,
> straightened out with classical and satellite observations. This "us
> army northern haiti datum of 1942" to name it, is local only and not
> connected to other datums.
>
> Surely if the people making the map had the technical ability to connect
> to NAD27, they would have done that. However Haiti is a long way from
> Florida :-)
>
> > One thousand meter Universal Transverse Mercator grid, Clarke 1866
> > spheroid, zone 18.
>
> This is a projection, not a datum, and simply describes how to get
> projected grid coordinates from geodetic (lat/lon) coordinates.
> Indeed it is a projection typically used with NAD27. NAD27 was the
> standard datum for the mainland.
>
> The "Thousand meter" pretty obviously refers to grid lines drawn on the
> map. UTM does not really have a scale like that.
>
> > The very best accuracy to WGS84 I've been able to get with the 1st
> edition
> > LOC map is ~1,000 meters offset. I think there must be some
> transformation
> > involved, but it isn't clear to me at what point in the process. Should
> the
> > Eastings and Northings in the map be transformed and use those
> transformed
> > values as GCPs? Or should the GCPs use the grid values and some transform
> > happen later? Or should is this just a really elusive horizontal datum...
>
> There is no reason to expect this datum to be aligned with WGS84.
>
> > I've tried everything I can think of. Here is a starter using lat long
> GCP
> > points on the LOC map image.
>
> You really aren't explaining where those points came from.
>
> What I would do is:
>
> First, I'd probably try to do this in qgis, and use the georeferencer.
>
> Examine the map and find the 3 points where they have published
> coordinates. Realize that these coordinates are in the map's (local)
> datum. Guess at where they measured, but they only gave integer
> seconds.
>
> For those points, find their coordinates in WGS84/ITRF/NAD83 by using
> modern georeferenced imagery. Much harder, go there and make 48h
> static dual-frequency GNSS observations.
>
> Compute distances from modern coordinates and also from local
> coordinates. See if the distances match. A sanity check.
>
> Understand how different these are. I'm guessing there will be
> disagreement at the 0.01 degree level.
>
> Perhaps, compute a transform from the local to the modern coordinates.
> But, you'll need to transform in UTM, which means from haiti-local-UTM
> to NAD27-UTM.
>
> Accept that you really aren't going to be able to get an accurate
> transform from the above. Assume that other than orientation and
> offset the map is pretty good quality (seems fair to me). Assume that
> the distortion from nad27-utm to nad83-utm due to the differing
> ellipsoids is small compared to errors in the map (also seems fair to
> me, at least for a first step).
>
> Find about 10 features on the map that are still present today and
> that you can identify on modern georeferenced imagery.
>
> Add the 10 GCPs in the qgis georeferencer, based on modern
> WGS84/ITRF/NAD83 of the features and clicking on the map.
>
> Let the georeferencer estimate the transform, and see what it looks
> like.
>
>
>
>
> Most importantly, after asking for others to help you on the list,
> write up careful notes about what you did and what you found and
> report back to the list. This is in my view the reciprocal obligation
> undertaken by asking.
>
>
> Hope this helps!
> Greg
> _______________________________________________
> PROJ mailing list
> PROJ at lists.osgeo.org
> https://lists.osgeo.org/mailman/listinfo/proj
>
--
Landon Yarrington, Ph.D.
Colorado State University
Department of Anthropology and Geography
Andrew G. Clark Bldg, 301 University Ave #219
Fort Collins, CO 80521
USA
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/proj/attachments/20240712/c133a728/attachment.htm>
More information about the PROJ
mailing list