<div dir="ltr">Even,<div>Thanks for the explanation and tips. I'm not clear on how setting to 26988 (Michigan North meters) would help. If I do a gdal_translate with that, it shows as being in New Zealand, just as Arc had done. It would certainly make lots more sense to actually put the data in a standard projection that has an EPSG code (meters or international feet), but that gets back to the embarrassing tale. It may well be that I'll have to force GeoTiff 1.0. I think there are other software packages that haven't caught up. Even when they do catch up, I've come across people running 15 year old CAD software and complaining about the DXF version GDAL produces.</div><div><br></div><div>Thanks to the whole team for the work you do to maintain GDAL.</div><div><br clear="all"><div><div dir="ltr" class="gmail_signature" data-smartmail="gmail_signature"><div dir="ltr">Kirk</div></div></div><br></div></div><br><div class="gmail_quote"><div dir="ltr" class="gmail_attr">On Thu, Jan 25, 2024 at 10:47 AM Even Rouault <<a href="mailto:even.rouault@spatialys.com">even.rouault@spatialys.com</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"><u></u>
<div>
<p>Kirk,</p>
<p>perhaps you could try using -a_srs EPSG:26988+6360 instead of a
full WKT. That way GDAL will *not* write the projected parameter
definitions, but just reference the horizontal and vertical CRS
codes:</p>
<p>Geotiff_Information:<br>
Version: 1<br>
Key_Revision: 1.1<br>
Tagged_Information:<br>
ModelTiepointTag (2,3):<br>
0 0 0 <br>
440720 3751320 0 <br>
ModelPixelScaleTag (1,3):<br>
60 60 1 <br>
End_Of_Tags.<br>
Keyed_Information:<br>
GTModelTypeGeoKey (Short,1): ModelTypeProjected<br>
GTRasterTypeGeoKey (Short,1): RasterPixelIsArea<br>
GTCitationGeoKey (Ascii,46): "NAD83 / Michigan North +
NAVD88 height (ftUS)"<br>
ProjectedCSTypeGeoKey (Short,1): PCS_NAD83_Michigan_North<br>
VerticalCSTypeGeoKey (Short,1): Code-6360 (NAVD88 height
(ftUS))<br>
End_Of_Keys.<br>
End_Of_Geotiff.<br>
</p>
<p>Even<br>
</p>
<div>Le 25/01/2024 à 16:17, Even Rouault via
gdal-dev a écrit :<br>
</div>
<blockquote type="cite">
<p>Kirk,</p>
<p>this is a complicated story indeed... What changed since <a href="https://trac.osgeo.org/gdal/changeset/31405" target="_blank">https://trac.osgeo.org/gdal/changeset/31405</a>
is that OGC GeoTIFF 1.1 was published, which mostly fixes issues
with vertical CRS, which your CRS uses. So by default GDAL 3.1
or later, when seeing a CompoundCRS write it as GeoTIFF 1.1, and
this new GeoTIFF version was the opportunity to remove writing
the ProjLinearUnitsInterpCorrectGeoKey=3059 hack, since using
GeoTIFF 1.1 is a way to indicate that you're implementing
correctly the standard.</p>
<p>If you force writing GeoTIFF 1.0 by adding -co
GEOTIFF_VERSION=1.0, GDAL will write
the ProjLinearUnitsInterpCorrectGeoKey when needed, and listgeo
will show the "Unknown-3059 (Short,1): Unknown-1" geotiff key.</p>
<p>To know if a Geotiff is 1.0 or 1.1 with listgeo , look at the
Key_Revision in the 3rd line which will be 1.0 or 1.1</p>
<p>So I assume ArcPro doesn't specifically implement 1.1. Either
they implement GeoTIFF themselves, or they use GDAL < 3.1
that hasn't GeoTIFF 1.1 explicit support</p>
<p>Even<br>
</p>
<div>Le 25/01/2024 à 15:24, Kirk Waters -
NOAA Federal via gdal-dev a écrit :<br>
</div>
<blockquote type="cite">
<div dir="ltr">Hi,
<div>I've run across an odd issue with GeoTIFFS using custom
projections written by GDAL and their interpretation in
ArcPro. The specific case is doing a projection for Michigan
State Plane North using U.S. survey feet. [insert
embarrassing tale of another US agency wanting to always use
survey feet regardless of state legislation or that the
federal government has been officially metric for a long
time]. For a file that should land in the Michigan upper
peninsula, ArcPro puts it in New Zealand. It appears the
issue is that Arc is interpreting the false easting as
meters instead of survey feet. If I do a define projection
in ArcPro and make a custom projection, it lands in the
right place. </div>
<div><br>
</div>
<div>I used an old listgeo from the geotiff library to look at
the tags both in the original file and a copy that had
define projection applied to see what was different. They
looked essentially the same, including having the same value
for the false easting. However, the Arc version had an extra
tag (3059) that wasn't in the GDAL produced one. An internet
search turned up <a href="https://trac.osgeo.org/gdal/ticket/3901" target="_blank">this 13-year old GDAL issue</a>. If
I understand it correctly, a bug in GDAL was found that
always used meters for the false easting and northing. In
addition to fixing that, a tag (3059) was added to indicate
this was fixed so software could differentiate broken GDAL
output from new output. It looks like Arc is looking for
this tag and if it's not there, false easting is assumed to
be meters.</div>
<div><br>
</div>
<div>A few questions, assuming I interpreted that issue
correctly:</div>
<div>1) Why is gdal not adding that 3059 tag? (GDAL command
used is below)</div>
<div>2) If it's only supposed to be relevant to GDAL written
files, why is Arc writing it?</div>
<div>3) How do I get GDAL to add the tag? </div>
<div><br>
</div>
<div>Command used in file creation:</div>
<div>gdal_translate -a_srs 'COMPD_CS[PROJCS["NAD83 / Michigan
North",GEOGCS["NAD83",DATUM["North_American_Datum_1983",SPHEROID["GRS
1980",6378137,298.257222101,AUTHORITY["EPSG","7019"]],TOWGS84[0,0,0,0,0,0,0],AUTHORITY["EPSG","6269"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4269"]],PROJECTION["Lambert_Conformal_Conic_2SP"],PARAMETER["standard_parallel_1",47.08333333333334],PARAMETER["standard_parallel_2",45.48333333333333],PARAMETER["latitude_of_origin",44.78333333333333],PARAMETER["central_meridian",-87],PARAMETER["false_easting",26246666.6666667],PARAMETER["false_northing",0],UNIT["US
survey
foot",0.304800609601219],AXIS["X",EAST],AXIS["Y",NORTH]],VERTCRS["NAVD88
height (ftUS)", VDATUM["North American Vertical Datum
1988"], CS[vertical,1], AXIS["gravity-related height
(H)",up], LENGTHUNIT["US survey foot",0.304800609601],
ID["EPSG",6360]]]' -co TILED=YES -co COMPRESS=DEFLATE -co
PREDICTOR=3 -co BIGTIFF=IF_SAFER input.tif output.tif <br>
</div>
<div><br clear="all">
<div>
<div dir="ltr" class="gmail_signature">
<div dir="ltr">
<div><font face="arial, helvetica, sans-serif">Kirk
Waters, PhD </font></div>
<div><font face="arial, helvetica, sans-serif">NOAA
Office for Coastal Management<br>
</font></div>
<div><font face="arial, helvetica, sans-serif">Applied
Sciences Program </font></div>
<div>
<div><font face="arial, helvetica, sans-serif"><a href="http://coast.noaa.gov/digitalcoast" target="_blank">coast.noaa.gov/digitalcoast</a></font></div>
</div>
<div><br>
</div>
<div><br>
</div>
</div>
</div>
</div>
</div>
</div>
<br>
<fieldset></fieldset>
<pre>_______________________________________________
gdal-dev mailing list
<a href="mailto:gdal-dev@lists.osgeo.org" target="_blank">gdal-dev@lists.osgeo.org</a>
<a href="https://lists.osgeo.org/mailman/listinfo/gdal-dev" target="_blank">https://lists.osgeo.org/mailman/listinfo/gdal-dev</a>
</pre>
</blockquote>
<pre cols="72">--
<a href="http://www.spatialys.com" target="_blank">http://www.spatialys.com</a>
My software is free, but my time generally not.</pre>
<br>
<fieldset></fieldset>
<pre>_______________________________________________
gdal-dev mailing list
<a href="mailto:gdal-dev@lists.osgeo.org" target="_blank">gdal-dev@lists.osgeo.org</a>
<a href="https://lists.osgeo.org/mailman/listinfo/gdal-dev" target="_blank">https://lists.osgeo.org/mailman/listinfo/gdal-dev</a>
</pre>
</blockquote>
<pre cols="72">--
<a href="http://www.spatialys.com" target="_blank">http://www.spatialys.com</a>
My software is free, but my time generally not.</pre>
</div>
</blockquote></div>