Re: [gdal-dev] [PATCH v2] Support Mercator_2SP in GeoTIFF

[Antti Castrén]

Hello Trent, Even et al.



Your work and comments on GeoTIFF and Mercator 2sp hopefully addresses the
problem I ran into few years ago. I wrote a description to ESRI's forums
with extracts of header information by ArcGIS and Geomedia. The postings
are at:

http://forums.arcgis.com/threads/24850-Projection-problem-with-GeoTIFFs-and-Mercator



Since then I have seen the similar behavior in other software also. It
confuses my geodesist-engineer-brain a little to see the terms latitude of
origin, latitude of true scale, standard parallel and natural origin
latitude used more or less interchangeably...



If you need volunteer for testing the patch or to give more insight to the
problem I am more than willing. Hopefully the software vendors follow the
resulting changes in specs.



To circumvent the whole problem I have advised my colleagues to reproject
their future Mercator_2sp GeoTIFFs to the Equator. This way the scale
problem disappears since both latitude parameter values are in fact 0. In
addition they can then use predefined CRS-code EPSG::3395, because nautical
charts are nowadays by international standards published using WGS84 as
horizontal datum.



Cheers,



  Antti
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev

Re: [gdal-dev] [PATCH v2] Support Mercator_2SP in GeoTIFF

[Antti Castrén]

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.


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]]


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):
 000
 -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): 000
  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.00 (  0d 0' 0.00"N)
   ProjNatOriginLongGeoKey: 0.00 (  0d 0' 0.00"E)
   ProjScaleAtNatOriginGeoKey: 1.00
   ProjFalseEastingGeoKey: 0.00 m
   ProjFalseNorthingGeoKey: 0.00 m
GCS: 4269/NAD83
Datum: 6269/North American Datum 1983
Ellipsoid: 7019/GRS 1980 (6378137.00,6356752.31)
Prime Meridian: 8901/Greenwich (0.00/  0d 0' 0.00"E)
Projection Linear Units: 9001/metre (1.00m)
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 :
> 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 li

Re: [gdal-dev] [PATCH v2] Support Mercator_2SP in GeoTIFF

[Antti Castrén]

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.767
  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 :
> 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  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.
>>
>> 
>> 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["Me

Re: [gdal-dev] gdalwarp EPSG:32662 problem

[Antti Castrén]

Hello Steve,

You asked for any thoughts, so here are some.

The 20km shift could easily be caused by using spherical version in
one implementation of equirectangular projection and ellipsoidal
version in the other. Amount of shift depends on the latitude.

You can find some information on these pages:

http://www.remotesensing.org/geotiff/proj_list/equirectangular.html

http://comments.gmane.org/gmane.comp.gis.proj-4.devel/981

Read especially the comments of the latter.

Hope these help,

  Antti

2014-02-08 8:11 GMT+02:00 Stephen Woodbridge :
> Hi,
>
> Sorry for cross posting. This seems like a proj4 issue but might be gdal
> related.
>
> I've run into a problem that I think is related to a change in the EPSG
> definition for EPSG:32662
>
> The symptom is that the image after gdalwarp reprojection is about 20km
> south of where it should be.
>
> old system is correctly aligned:
>
> GDAL 1.4.2.0, released 2007/06/27
> proj Rel. 4.5.0, 22 Oct 2006
>
> # WGS 84 / Plate Carree
> <32662> +proj=eqc +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84
> +units=m +no_defs  <>
>
> new system is shifted about 20km south:
>
> GDAL 1.9.2, released 2012/10/08
> proj Rel. 4.8.0, 6 March 2012
>
> # WGS 84 / Plate Carree (deprecated)
> <32662> +proj=eqc +lat_ts=0 +lat_0=0 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84
> +units=m +no_defs  <>
>
>
> http://spatialreference.org/ref/epsg/32662/
> This indicates that the definition in proj 4.8.0 is wrong. Would the
> +ellps=WGS84 account for this shift?
>
> Any thoughts on this?
>
> Thanks,
>   -Steve
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev

Re: [gdal-dev] gdalwarp EPSG:32662 problem

[Antti Castrén]

2014-02-09 8:40 GMT+02:00 Andre Joost :
> Am 09.02.2014 00:42, schrieb Even Rouault:
>>
>> But Antti guess seems right. Instead of +ellps=WGS84 (or +datum=WGS84), if
>> you
>> play with the a (semi-major axis) and b (semi-minor axis) parameters, you
>> can
>> see that only +a has an influence, so latest proj version seems to use a
>> spherical version of eqc.
>
>
> If you look at
> 
>
> and the chapter 12 of Snyders manual, you will only find formulas for the
> sphere. So I guess there is no other way to calculate eqc.
>
> Maybe older versions calculated another radius for the sphere when an
> ellipsoid was given.

Stephen's "shift" was about 20km south, which correlates quite well if
you use semi-minor axis of WGS84 as radius of sphere while calculating
forward, and semi-major axis as radius of sphere while calculating the
inverse. At latitude of 55 degrees the difference is ca. 20 530 meters
(55 degrees -> 54.8156 degrees).

There are several different radii of the Earth, and some of them could
arguably be used in this context in place of semi-major axis.

All radii ar not created equal: http://en.wikipedia.org/wiki/Earth_radius

You have to bear in mind that this projection is intended for small
scale mapping, for example mapping the whole world. In that scale 20
km is nothing. If you need better representation of the Earth you have
to use a projection which takes ellipsoidal properties into account.
Of course the beef in this thread is not about choosing a projection,
but the change/difference in formulae used, which can create problems
as Stephen pointed out.

Cheers,

  Antti
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev

Re: [gdal-dev] [PATCH v2] Support Mercator_2SP in GeoTIFF

[Antti Castrén]

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.333
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.333
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.333 (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 :
> 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
>   Proj

Re: [gdal-dev] S57 and area reconstruction

[Antti Castrén]

Hello Frits,

It is possible to (re)merge the objects across cell boundaries if they
share common Feature Object IDentifier (FOID), which is a combination
of AGEN, FIND and FINS subfields. Some ECDISes can hide the boundary
depending on masking. Usage of FOIDs vary, so it might be possible for
your datasets - if you are lucky.

See Clause 3.1 of Edition 2.0 (November 2000) of the ENC Product
Specification (S-57 Appendix B.1) and ENC Encoding Bulletin Number 33.

Cheers,

  Antti


2014-09-11 8:33 GMT+03:00 vwf :
> Hello,
>
> My question is about S57, not gdal, but one of you may be able to help
> me.
>
> When an area (like "prcare") is at the edge of a cell (map) there will
> be a vertex on the edge of the cel. With several cells the full area can
> have several lines going through the area. Is there a trick/technique to
> reconstruct the outer edge or remove the extra lines?
>
> Thank you,
> Frits
>
> ___
> gdal-dev mailing list
> gdal-dev@lists.osgeo.org
> http://lists.osgeo.org/mailman/listinfo/gdal-dev



-- 
Antti Castrén
+358505997538
antti.cast...@iki.fi
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev