<div dir="ltr">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.<br><br><div>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?</div><br>Fort Islet lighthouse 18°33'31.33" -72°20'59.03", should be 18°33'18.08" -72°21'04.50"<br>Santo Domingo lighthouse 18°27'53.64" -69°52'52.74", should be 18°27'51.39" -69°52'34.90"<br>Dame Marie 18°36'47" -74°25'53", unknown but possibly 18°36'57.21" -74°25'26.94"<br><br>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?<br><br>The map in question is the 1st edition (<a href="https://www.loc.gov/resource/g4940m.gct00206/?sp=168" target="_blank">https://www.loc.gov/resource/g4940m.gct00206/?sp=168</a>). There are at least 4 subsequent editions (here is the 3rd <a href="https://www.loc.gov/resource/g4940m.gct00206/?sp=167" target="_blank">https://www.loc.gov/resource/g4940m.gct00206/?sp=167</a>). 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 <a href="https://www.loc.gov/resource/g4940m.gct00206/?sp=168&r=0.314,0.573,0.372,0.229,0" target="_blank">https://www.loc.gov/resource/g4940m.gct00206/?sp=168&r=0.314,0.573,0.372,0.229,0</a> to <a href="https://www.loc.gov/resource/g4940m.gct00206/?sp=167&r=0.284,0.495,0.238,0.147,0" target="_blank">https://www.loc.gov/resource/g4940m.gct00206/?sp=167&r=0.284,0.495,0.238,0.147,0</a>. What is going on here?? Is this a datum difference, projection difference, both??<br><br>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.<br><br>Again, I am indebted for any help. The code below only uses 8 GCPs of the UTM intersections offset and shifted described above.<br><br>```<br># download the tiff from Library of Congress<br>wget <a href="https://www.loc.gov/resource/g4940m.gct00206/?download=https%3A%2F%2Ftile.loc.gov%2Fstorage-services%2Fmaster%2Fgmd%2Fgmd4m%2Fg4940m%2Fg4940m%2Fgct00206%2Fcs000168.tif" target="_blank">https://www.loc.gov/resource/g4940m.gct00206/?download=https%3A%2F%2Ftile.loc.gov%2Fstorage-services%2Fmaster%2Fgmd%2Fgmd4m%2Fg4940m%2Fg4940m%2Fgct00206%2Fcs000168.tif</a> -O cs000168.tif<br># add 8 ground control points using the UTM grid described above<br>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<br># warp the map with thin plate spline<br>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<br>```</div><br><div class="gmail_quote"><div dir="ltr" class="gmail_attr">On Tue, Jul 2, 2024 at 4:40 AM Greg Troxel via PROJ <<a href="mailto:proj@lists.osgeo.org" target="_blank">proj@lists.osgeo.org</a>> wrote:<br></div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex">Landon Yarrington via PROJ <<a href="mailto:proj@lists.osgeo.org" target="_blank">proj@lists.osgeo.org</a>> writes:<br>
<br>
> I am trying to georeference an US Army map of northern Haiti made from<br>
> aerial photographs taken in 1942 and 1944 using `gdal_translate` and<br>
> `gdalwarp`. I need a sanity check because this is impossible, and so I'm<br>
> asking for help. Here is the map from Lib. of Congress (<br>
> <a href="https://www.loc.gov/resource/g4940m.gct00206/?sp=160&r=0.365,0.623,0.281,0.173,0" rel="noreferrer" target="_blank">https://www.loc.gov/resource/g4940m.gct00206/?sp=160&r=0.365,0.623,0.281,0.173,0</a><br>
> ).<br>
><br>
> The map info states,<br>
<br>
> Horizontal Datum is based on the following astronomic values:<br>
> Fort Islet Lighthouse Lat. 18°33'31.33", Long. 72°20'59.03"<br>
> Santo Domingo Lighthouse Lat. 18°27'53.64", Long. 69°52'52.74"<br>
> Cape Dame Marie Astro Lat. 18°36'47", Long. 74°25'53"<br>
> Extensions of control from Astronomic Positions have been made by multiplex<br>
> triangulation.<br>
<br>
This is remarkably clear and useful even if it is hard to understand.<br>
<br>
The mapmaker is saying that they made astronomic determinations of<br>
position at 3 named locations, presumably obvious on the map, and gives<br>
the values. These coordinates are not related to any particular datum.<br>
Instead, they define a local datum.<br>
<br>
Critically, there is absolutely no reason to expect that the coordinates<br>
above match WGS84, other than that the time signals are consistent with<br>
greenwhich, and latitude is physical. But that doesn't get you the last<br>
10s to 100s of meters.<br>
<br>
Almost certainly they had the benefit of radio time signals, so their<br>
longitudes should be pretty good. Latitudes should be pretty good<br>
without needing to think about time. Without really thinking, 1000m for<br>
a single determination strikes me as pretty good. Before NAD83, which<br>
used satellite/etc. data, NAD27 was different from modern ITRF by 40m<br>
around me. And that's a national datum that got to average multiple<br>
observations, even if it was held to be somewhat consistent with<br>
mid/late 1800s work.<br>
<br>
Our current datums derive from averages of many astronomic observations,<br>
straightened out with classical and satellite observations. This "us<br>
army northern haiti datum of 1942" to name it, is local only and not<br>
connected to other datums.<br>
<br>
Surely if the people making the map had the technical ability to connect<br>
to NAD27, they would have done that. However Haiti is a long way from<br>
Florida :-)<br>
<br>
> One thousand meter Universal Transverse Mercator grid, Clarke 1866<br>
> spheroid, zone 18.<br>
<br>
This is a projection, not a datum, and simply describes how to get<br>
projected grid coordinates from geodetic (lat/lon) coordinates.<br>
Indeed it is a projection typically used with NAD27. NAD27 was the<br>
standard datum for the mainland.<br>
<br>
The "Thousand meter" pretty obviously refers to grid lines drawn on the<br>
map. UTM does not really have a scale like that.<br>
<br>
> The very best accuracy to WGS84 I've been able to get with the 1st edition<br>
> LOC map is ~1,000 meters offset. I think there must be some transformation<br>
> involved, but it isn't clear to me at what point in the process. Should the<br>
> Eastings and Northings in the map be transformed and use those transformed<br>
> values as GCPs? Or should the GCPs use the grid values and some transform<br>
> happen later? Or should is this just a really elusive horizontal datum...<br>
<br>
There is no reason to expect this datum to be aligned with WGS84.<br>
<br>
> I've tried everything I can think of. Here is a starter using lat long GCP<br>
> points on the LOC map image.<br>
<br>
You really aren't explaining where those points came from.<br>
<br>
What I would do is:<br>
<br>
First, I'd probably try to do this in qgis, and use the georeferencer.<br>
<br>
Examine the map and find the 3 points where they have published<br>
coordinates. Realize that these coordinates are in the map's (local)<br>
datum. Guess at where they measured, but they only gave integer<br>
seconds.<br>
<br>
For those points, find their coordinates in WGS84/ITRF/NAD83 by using<br>
modern georeferenced imagery. Much harder, go there and make 48h<br>
static dual-frequency GNSS observations.<br>
<br>
Compute distances from modern coordinates and also from local<br>
coordinates. See if the distances match. A sanity check.<br>
<br>
Understand how different these are. I'm guessing there will be<br>
disagreement at the 0.01 degree level.<br>
<br>
Perhaps, compute a transform from the local to the modern coordinates.<br>
But, you'll need to transform in UTM, which means from haiti-local-UTM<br>
to NAD27-UTM.<br>
<br>
Accept that you really aren't going to be able to get an accurate<br>
transform from the above. Assume that other than orientation and<br>
offset the map is pretty good quality (seems fair to me). Assume that<br>
the distortion from nad27-utm to nad83-utm due to the differing<br>
ellipsoids is small compared to errors in the map (also seems fair to<br>
me, at least for a first step).<br>
<br>
Find about 10 features on the map that are still present today and<br>
that you can identify on modern georeferenced imagery.<br>
<br>
Add the 10 GCPs in the qgis georeferencer, based on modern<br>
WGS84/ITRF/NAD83 of the features and clicking on the map.<br>
<br>
Let the georeferencer estimate the transform, and see what it looks<br>
like.<br>
<br>
<br>
<br>
<br>
Most importantly, after asking for others to help you on the list,<br>
write up careful notes about what you did and what you found and<br>
report back to the list. This is in my view the reciprocal obligation<br>
undertaken by asking.<br>
<br>
<br>
Hope this helps!<br>
Greg<br>
_______________________________________________<br>
PROJ mailing list<br>
<a href="mailto:PROJ@lists.osgeo.org" target="_blank">PROJ@lists.osgeo.org</a><br>
<a href="https://lists.osgeo.org/mailman/listinfo/proj" rel="noreferrer" target="_blank">https://lists.osgeo.org/mailman/listinfo/proj</a><br>
</blockquote></div><br clear="all"><br><span class="gmail_signature_prefix">-- </span><br><div dir="ltr" class="gmail_signature"><div dir="ltr"><div><span style="color:rgb(0,0,0)">Landon Yarrington, Ph.D.<br></span></div><div><span style="color:rgb(0,0,0)">Colorado State University</span></div><div><span style="color:rgb(0,0,0)">Department of Anthropology and Geography</span></div><div><span style="color:rgb(0,0,0)">Andrew G. Clark Bldg, 301 University Ave #219<br></span></div><div><span style="color:rgb(0,0,0)">Fort Collins, CO 80521</span></div><div><span style="color:rgb(0,0,0)">USA<br></span></div></div></div>