[gdal-dev] [PATCH v2] Support Mercator_2SP in GeoTIFF
Antti Castrén
antti.castren at iki.fi
Wed Aug 6 05:52:32 PDT 2014
Hello again,
Because of few customer cases I had to look this GeoTIFF/CT_Mercator
issue once again, this time in depth.
I studied how three software packages write GeoTIFFs as well as how
they read different combination of keys. These software are
GeoMedia Professional 06.01.11.13
ArcGIS Desktop 10 SP4
FME Desktop 2013 SP1
Test case is a Finnish navigational chart, projected mercator with
standard parallel (ie. latitude of true scale) 60d20m on GRS80
ellipsoid.
Geomedia writes following projection related keys:
ProjNatOriginLongGeoKey = 0
ProjNatOriginLatGeoKey = 0
ProjFalseEastingGeoKey = 0
ProjFalseNorthingGeoKey = 0
ProjCenterLatGeoKey = 60.3333333
ProjScaleAtNatOriginGeoKey = 0.496208843
Geomedia clearly writes
http://www.remotesensing.org/geotiff/proj_list/mercator_1sp.html type
of keys with an odd addition of ProjCenterLatGeoKey having the value
of standard parallel.
ArcGIS writes following projection related keys:
ProjNatOriginLongGeoKey = 0
ProjNatOriginLatGeoKey = 60.3333333
ProjFalseEastingGeoKey = 0
ProjFalseNorthingGeoKey = 0
ProjScaleAtNatOriginGeoKey = 1
ArcGIS writes something more like
http://www.remotesensing.org/geotiff/proj_list/mercator_2sp.html type
of information, but in wrong keys. Standard parallel should be written
in ProjStdParallel1GeoKey and ProjNatOriginLatGeoKey should equal 0.
Also the scale factor is unnecessary and actually it is harmful since
it gives wrong scale value for correct origin latitude.
FME writes the same way as ArcGIS.
After trying some 20 different GeoKey combinations using geotifcp I
was able to find following behaviour in reading the projection
information:
1. ProjStdParallel1GeoKey is read by ArcGIS only, and it is primary
source of scaling information for ArcGIS.
2. If ProjStdParallel1GeoKey is not populated then ArcGIS tries to
find latitude of true scale in ProjNatOriginLatGeoKey
3. If neither one of previous keys are not populated then ArcGIS tries
to find latitude of true scale in ProjCenterLatGeoKey.
4. ArcGIS does not read ProjScaleAtNatOriginGeoKey at all.
5. Geomedia needs ProjScaleAtNatOriginGeoKey and either
ProjNatOriginLatGeoKey or ProjCenterLatGeoKey to open the file in the
first place.
6. Geomedia reads ProjNatOriginLatGeoKey and ProjCenterLatGeoKey as
origin information (which is the correct way IMHO)
7. FME reads otherwise as ArcGIS, but it does not read ProjStdParallel1GeoKey
Since I did not study ALL the possible combinations of relevant
GeoKeys there could be some errors, but I doubt it.
So my conclusion is that it is possible to write GeoTIFF in a way that
both Geomedia and ArcGIS read it the correct way. It is achieved by
populating following keys:
ProjStdParallel1GeoKey = 60.3333333 (this makes ArcGIS happy)
ProjNatOriginLongGeoKey = 0 (Geomedia needs this)
ProjNatOriginLatGeoKey = 0
ProjFalseEastingGeoKey = 0
ProjFalseNorthingGeoKey = 0
ProjScaleAtNatOriginGeoKey = 0.496208843 (Geomedia needs this)
This is also correct combination of mercator1sp and mercator2sp
Unfortunately FME does not fit in since it reads
ProjNatOriginLatGeoKey value as latitude of true scale the same time
as Geomedia treats it as latitude of origin. Therefore FME shows the
chart to be in Tunisia instead of Gulf of Finland.
Other bad combinations are the ones where ProjNatOriginLatGeoKey has
the value of latitude of true scale and the Geomedia uses that value
as latitude of origin and either scale factor 1 or 0.496208843.
Therefore Geomedia shows the chart to be in vicinity of Spitsbergen
either South or North of the islands respectively.
I don't know if these findings have any impact on development of gdal,
but at least they give some insight to the issue. Hopefully we see
some development on the above mentioned GIS software also.
Cheers,
Antti
2013-11-23 0:18 GMT+02:00 Antti Castrén <antti.castren at iki.fi>:
> Hi Trent et al.
>
> Thanks for your work on the issue.
>
> I tested the files, and results were not satisfactory. But at least we
> got more information on the issue.
>
> FME opened both files.
>
> The georeferencing of 18440_lowres_nol.tif was recognized by FME as:
> DESC_NM: Global Mercator
> DT_NAME: NAD83
> PARM2: 47.667
> PROJ: MRCAT
> UNIT: METER
>
> This should place the file at correct location.
>
> The georeferencing of 18440_lowres_sp1.tif was recognized by FME as:
> DESC_NM: Global Mercator
> DT_NAME: NAD83
> PROJ: MRCAT
> UNIT: METER
>
> To compare to the first file which you provided last week, its FME
> parameters were:
> DESC_NM: Global Mercator
> DT_NAME: NAD83
> PROJ: MRCAT
> UNIT: METER
>
> It seems FME doesn't read StdParallel1 or at least it doesn't use it.
>
> Geomedia did NOT recognize either of the files as valid GeoTIFF.
>
> To get projection (and of course location) right GeoMedia reads following keys:
> ProjCoordTransGeoKey (Short,1): CT_Mercator
> ProjLinearUnitsGeoKey (Short,1): Linear_Meter
> ProjNatOriginLongGeoKey (Double,1): 0
> ProjNatOriginLatGeoKey (Double,1): 0
> ProjFalseEastingGeoKey (Double,1): 0
> ProjFalseNorthingGeoKey (Double,1): 0
> ProjCenterLatGeoKey (Double,1): 60.7666667
> ProjScaleAtNatOriginGeoKey (Double,1): 0.48961699
>
> These keys are (in my humble opinion) logical and the most correct
> ones. Of course there is some redundancy. Unfortunately FME and ArcGIS
> don't read last two at all, and therefore they place the file as if it
> was projected from the Equator with scale 0.
>
> ArcGIS opened all your three files at correct location.
>
> By the way, I have seen that many sofware which are able to
> interchange latitude of true scale and scalefactor at the Equator fail
> to do the calculation with enough decimal places. My experience is
> that you need at least 8 decimal places in scale factor to align your
> raster data correctly. (As seen in the GeoMedia example above)
>
> Br,
>
> Antti
>
> 2013/11/20 Trent Piepho <tpiepho at gmail.com>:
>> Antii,
>>
>> Thank you for the testing. It looks like Arc/GIS works while FME and
>> Intergraph don't.
>>
>> That file had both StdParallel1 and ScaleAtNatOrigin defined.
>>
>> I've uploaded two new files. This one has only StdParallel1 defined with no
>> ScaleAtNatOrigin.
>> https://drive.google.com/file/d/0BybuTedE9CLxRnNlNUl2Zm5JRzg/edit?usp=sharing
>>
>> This file has the lat_ts stored in NatOriginLat and no ScaleAtNatOrigin
>> defined.
>> https://drive.google.com/file/d/0BybuTedE9CLxMzltYlBFWU51RzA/edit?usp=sharing
>>
>> I think the latter might have better compatibility. Yet it seems to clearly
>> be the incorrect way to do it.
>>
>>
>> On Mon, Nov 18, 2013 at 4:54 AM, Antti Castrén <antti.castren at iki.fi> wrote:
>>>
>>> Hi Trent,
>>>
>>> The chart opens in ArcMap (ArcGIS 10.0) well, and it is in right
>>> location (Seattle).
>>>
>>> Relevant properties of the file as seen by ArcGIS:
>>> Cell Size: 240.003787, 240.003787
>>> Extent
>>> Top:4143530.20898
>>> Left: -9278434.53415
>>> Right: -9165872.75807
>>> Bottom: 3984167.69444
>>> Spatial Reference: Global Mercator
>>> Linear Unit: Meter
>>> false_easting: 0
>>> false_northing: 0
>>> central_meridian: 0
>>> standard_parallel_1: 47.667
>>> Datum: D_North_American_1983
>>>
>>>
>>> I opened the file in FME Data Inspector (FME Desktop 2013 SP1) also.
>>> It did not read the latitude of true scale from the
>>> ProjStdParallel1GeoKey as you can see from the following parameters.
>>> Therefore it places the chart to NE Georgia/NW South Carolina as if
>>> the latitude of true scale was the Equator.
>>>
>>> <FME Parameters>
>>> Coordinate System Parameters
>>> CS_NAME: _FME_0
>>> DESC_NM: Global Mercator
>>> DT_NAME: NAD83
>>> PROJ: MRCAT
>>> UNIT: METER
>>>
>>> Datum Parameters
>>> DESC_NM: NAD 1983, Alaska, Canada,Continental US,Mexico,Central America
>>> ELLIPSOID: GRS1980
>>> SOURCE: US Defense Mapping Agency, TR-8350.2-B, August 1993
>>> USE: NAD83
>>>
>>> Ellipsoid Parameters
>>> DESC_NM: Geodetic Reference System of 1980
>>> E_RAD: 6378137
>>> P_RAD: 6356752.3141403478
>>> SOURCE: Stem, L.E., Jan 1989, State Plane Coordinate System of 1983
>>>
>>> OGC WKT Description
>>> PROJCS["Global Mercator",
>>> GEOGCS["NAD83",
>>> DATUM["North_American_Datum_1983",
>>> SPHEROID["Geodetic Reference System of
>>> 1980",6378137,298.2572221008916,
>>> AUTHORITY["EPSG","7019"]],
>>> AUTHORITY["EPSG","6269"]],
>>> PRIMEM["Greenwich",0],
>>> UNIT["degree",0.0174532925199433],
>>> AUTHORITY["EPSG","4269"]],
>>> PROJECTION["Mercator_1SP"],
>>> PARAMETER["central_meridian",0],
>>> PARAMETER["scale_factor",1],
>>> PARAMETER["false_easting",0],
>>> PARAMETER["false_northing",0],
>>> UNIT["METER",1]]
>>>
>>> Esri WKT Description
>>>
>>> PROJCS["Global_Mercator",GEOGCS["GCS_North_American_1983",DATUM["D_North_American_1983",SPHEROID["Geodetic_Reference_System_of_1980",6378137,298.2572221008916]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]],PROJECTION["Mercator"],PARAMETER["central_meridian",0],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["Meter",1],PARAMETER["standard_parallel_1",0.0]]
>>> </FME Parameters>
>>>
>>> Intergraph's Geomedia Professional recognises the file as GeoTIFF, but
>>> the chart appears to be in the same false place as with FME.
>>>
>>> If you want, I could update the file's georeferencing by FME and
>>> Geomedia so that location is correct. Then I could send listgeo
>>> information or even put the actual files available if desired.
>>>
>>> The goal is, of course, to get various software vendors to use the
>>> GeoTIFF-information the same way. Hopefully the common way would also
>>> be logically correct, but I guess that's too much to ask.
>>>
>>> Best regards,
>>>
>>> Antti
>>>
>>>
>>> Here is the listgeo info from the original file for those who are
>>> interested:
>>>
>>> Geotiff_Information:
>>> Version: 1
>>> Key_Revision: 1.0
>>> Tagged_Information:
>>> ModelTiepointTag (2,3):
>>> 0 0 0
>>> -9278434.53 4143530.21 0
>>> ModelPixelScaleTag (1,3):
>>> 240.003787 240.003787 0
>>> End_Of_Tags.
>>> Keyed_Information:
>>> GTModelTypeGeoKey (Short,1): ModelTypeProjected
>>> GTRasterTypeGeoKey (Short,1): RasterPixelIsArea
>>> GTCitationGeoKey (Ascii,16): "Global Mercator"
>>> GeographicTypeGeoKey (Short,1): GCS_NAD83
>>> GeogCitationGeoKey (Ascii,6): "NAD83"
>>> GeogAngularUnitsGeoKey (Short,1): Angular_Degree
>>> GeogSemiMajorAxisGeoKey (Double,1): 6378137
>>> GeogInvFlatteningGeoKey (Double,1): 298.257222
>>> Unknown-2062 (Double,3): 0 0 0
>>> ProjectedCSTypeGeoKey (Short,1): User-Defined
>>> ProjectionGeoKey (Short,1): User-Defined
>>> ProjCoordTransGeoKey (Short,1): CT_Mercator
>>> ProjLinearUnitsGeoKey (Short,1): Linear_Meter
>>> ProjStdParallel1GeoKey (Double,1): 47.667
>>> ProjNatOriginLongGeoKey (Double,1): 0
>>> ProjNatOriginLatGeoKey (Double,1): 0
>>> ProjFalseEastingGeoKey (Double,1): 0
>>> ProjFalseNorthingGeoKey (Double,1): 0
>>> ProjScaleAtNatOriginGeoKey (Double,1): 1
>>> End_Of_Keys.
>>> End_Of_Geotiff.
>>> Projection Method: CT_Mercator
>>> ProjNatOriginLatGeoKey: 0.000000 ( 0d 0' 0.00"N)
>>> ProjNatOriginLongGeoKey: 0.000000 ( 0d 0' 0.00"E)
>>> ProjScaleAtNatOriginGeoKey: 1.000000
>>> ProjFalseEastingGeoKey: 0.000000 m
>>> ProjFalseNorthingGeoKey: 0.000000 m
>>> GCS: 4269/NAD83
>>> Datum: 6269/North American Datum 1983
>>> Ellipsoid: 7019/GRS 1980 (6378137.00,6356752.31)
>>> Prime Meridian: 8901/Greenwich (0.000000/ 0d 0' 0.00"E)
>>> Projection Linear Units: 9001/metre (1.000000m)
>>> Corner Coordinates:
>>> Upper Left (-9278434.534, 4143530.209)
>>> Lower Left (-9278434.534, 3984167.694)
>>> Upper Right (-9165872.758, 4143530.209)
>>> Lower Right (-9165872.758, 3984167.694)
>>> Center (-9222153.646, 4063848.952)
>>>
>>>
>>> 2013/11/14 Trent Piepho <tpiepho at gmail.com>:
>>> > It looks like arcgis does not support Mercator_1SP in geotiff files.
>>> > It just ignores the scale factor and uses 1.0. Which means in GDAL
>>> > and libgeotiff as they are now, there is no way to export Mercator
>>> > projections that don't have a scale of 1.0 to arcgis.
>>> >
>>> > ArcGIS does support Mercator_2SP, but instead of placing the latitude
>>> > of true scale in ProjStdParallel1GeoKey, it uses
>>> > ProjNatOriginLatGeoKey. Which seems to be plainly wrong, yet that's
>>> > what it does.
>>> >
>>> > I wonder if ESRI used libgeotiff? With unpatched libgeotiff, the only
>>> > way to get a latitude from a Mercator geotiff it to use one of the
>>> > origin latitude GeoKeys. It's not correct, but it does allow one to
>>> > stick the number in the file and get it back out without having to
>>> > modify or understand any libgeotiff code.
>>> >
>>> > It's possible that ArcGIS supports getting the latitude of true scale
>>> > from ProjStdParallel1GeoKey, even it it doesn't store it there.
>>> > Antii, could you test this file with ArcGIS?
>>> >
>>> >
>>> > https://drive.google.com/file/d/0BybuTedE9CLxak9iakZRTC00cHM/edit?usp=sharing
>>> >
>>> > Pick file->download to get the geotiff file.
>>> >
>>> > It should have projection information like this:
>>> > PROJECTION["Mercator_2SP"],
>>> > PARAMETER["standard_parallel_1",47.667],
>>> > PARAMETER["central_meridian",0],
>>> > PARAMETER["false_easting",0],
>>> > PARAMETER["false_northing",0],
>>> > UNIT["metre",1,
>>> > AUTHORITY["EPSG","9001"]]]
>>> > Origin = (-9278434.534151820465922,4143530.208984711207449)
>>> > Pixel Size = (240.003786961970206,-240.003786961970206)
>>> > Corner Coordinates:
>>> > Upper Left (-9278434.534, 4143530.209) (123d32'26.22"W, 48d23'56.59"N)
>>> > Lower Left (-9278434.534, 3984167.694) (123d32'26.22"W, 46d57'58.64"N)
>>> > Upper Right (-9165872.758, 4143530.209) (122d 2'30.76"W, 48d23'56.59"N)
>>> > Lower Right (-9165872.758, 3984167.694) (122d 2'30.76"W, 46d57'58.64"N)
>>> > Center (-9222153.646, 4063848.952) (122d47'28.49"W, 47d41'15.43"N)
More information about the gdal-dev
mailing list