You are viewing a plain text version of this content. The canonical link for it is here.
Posted to dev@sis.apache.org by Thierry Danard <th...@geotoolkit.net.INVALID> on 2022/02/11 19:38:30 UTC
Support for GeoTIFF ProjCoordTransGeoKey values from 1 to 27
Hello,
I started using the org.apache.sis.storage.sis-geotiff module (version 1.1)
for my application.
I tested it with various public domain GeoTIFF files (
https://download.osgeo.org/geotiff/samples/), and got this exception a few
times:
org.opengis.util.NoSuchIdentifierException: No operation method found for
name or identifier “GeoTIFF:XX”, where XX varies from 1 to 27. Each one of
these values is an entry in the 6.3.3.3 section of this document:
http://geotiff.maptools.org/spec/geotiff6.html#6.3.3.3
This exception comes from the CRSBuilder class:
https://github.com/apache/sis/blob/master/storage/sis-geotiff/src/main/java/org/apache/sis/storage/geotiff/CRSBuilder.java
After initially sharing my findings with Martin, I understand that the
support for the ProjCoordTransGeoKey values from 1 to 27 is not complete
yet.
For full support, I also understand that Apache SIS would need to support
the following projection methods:
case "2": // CT_TransvMercator_Modified_Alaska
case "5": // CT_ObliqueMercator_Rosenmund
case "6": //CT_ObliqueMercator_Spherical
case "9": // CT_LambertConfConic_Helmert
case "13": // CT_EquidistantConic
case "14": // CT_Stereographic
case "19": // CT_Gnomonic
case "21": // CT_Orthographic
case "23": // CT_Robinson
case "25": // CT_VanDerGrinten
I am proposing to implement partial support of the following projection
methods:
case "1": // CT_TransverseMercator
case "3": // CT_ObliqueMercator
case "4": // CT_ObliqueMercator_Laborde
case "7": // CT_Mercator
case "8": // CT_LambertConfConic_2SP
case "10": // CT_LambertAzimEqualArea
case "11": // CT_AlbersEqualArea
case "12": // CT_AzimuthalEquidistant, using Modified
Azimuthal Equidistant
case "15": // CT_PolarStereographic
case "16": // CT_ObliqueStereographic
case "17": // CT_Equirectangular
case "18": // CT_CassiniSoldner
case "26": // CT_NewZealandMapGrid
case "27": // CT_TransvMercator_SouthOriented
case "20": // CT_MillerCylindrical
case "22": // CT_Polyconic
case "24": // CT_Sinusoidal
Using http://geotiff.maptools.org/proj_list/ as a resource, I came up with
this code as part of the CRSBuilder.createConversion method:
switch (type) { // based on
http://geotiff.maptools.org/proj_list/
case "1": // CT_TransverseMercator
method =
getCoordinateOperationFactory().getOperationMethod("EPSG:9807");
break;
case "3": // CT_ObliqueMercator
method =
getCoordinateOperationFactory().getOperationMethod("EPSG:9815");
break;
case "4": // CT_ObliqueMercator_Laborde
method =
getCoordinateOperationFactory().getOperationMethod("EPSG:9813");
break;
case "7": // CT_Mercator
method =
getCoordinateOperationFactory().getOperationMethod("EPSG:9804");
break;
case "8": // CT_LambertConfConic_2SP
method =
getCoordinateOperationFactory().getOperationMethod("EPSG:9802");
break;
case "10": // CT_LambertAzimEqualArea
method =
getCoordinateOperationFactory().getOperationMethod("EPSG:9820");
break;
case "11": // CT_AlbersEqualArea
method =
getCoordinateOperationFactory().getOperationMethod("EPSG:9822");
break;
case "12": // CT_AzimuthalEquidistant, using Modified
Azimuthal Equidistant
method =
getCoordinateOperationFactory().getOperationMethod("EPSG:9832");
break;
case "15": // CT_PolarStereographic
method =
getCoordinateOperationFactory().getOperationMethod("EPSG:9810");
break;
case "16": // CT_ObliqueStereographic
method =
getCoordinateOperationFactory().getOperationMethod("EPSG:9809");
break;
case "17": // CT_Equirectangular
method =
getCoordinateOperationFactory().getOperationMethod("EPSG:9823");
break;
case "18": // CT_CassiniSoldner
method =
getCoordinateOperationFactory().getOperationMethod("EPSG:9806");
break;
case "26": // CT_NewZealandMapGrid
method =
getCoordinateOperationFactory().getOperationMethod("EPSG:9811");
break;
case "27": // CT_TransvMercator_SouthOriented
method =
getCoordinateOperationFactory().getOperationMethod("EPSG:9808");
break;
case "20": // CT_MillerCylindrical
method = new MillerCylindrical();
break;
case "22": // CT_Polyconic
method = new Polyconic();
break;
case "24": // CT_Sinusoidal
method = new Sinusoidal();
break;
case "2": // CT_TransvMercator_Modified_Alaska (not
supported by Apache SIS)
case "5": // CT_ObliqueMercator_Rosenmund (not
supported by Apache SIS)
case "6": //CT_ObliqueMercator_Spherical (not supported
by Apache SIS)
case "9": // CT_LambertConfConic_Helmert (not supported
by Apache SIS)
case "13": // CT_EquidistantConic (not supported by
Apache SIS)
case "14": // CT_Stereographic (not supported by Apache
SIS)
case "19": // CT_Gnomonic (not supported by Apache SIS)
case "21": // CT_Orthographic (not supported by Apache
SIS)
case "23": // CT_Robinson (not supported by Apache SIS)
case "25": // CT_VanDerGrinten (not supported by Apache
SIS)
default:
method =
getCoordinateOperationFactory().getOperationMethod(Constants.GEOTIFF + ':'
+ type);
break;
}
Martin proposed that I share my findings on this forum. He also suggested
that there are potential issues with this approach, that I will summarize
here:
1/ http://geotiff.maptools.org/proj_list/ is not authoritative
2/ it's not certain whether Azimuthal Equidistant and Modified Azimuthal
Equidistant are a match in case #12
3/ having a hard-coded mapping might not be the best way of doing things.
External third parties might want to customize this to add their own
mapping. Or these mappings could be discovered from the classpath at
runtime.
I am opening this discussion to see what you guys think. I am volunteering
to submit this patch to the CRSBuilder class if there is a consensus that
this is the right change.
-Thierry
Re: Support for GeoTIFF ProjCoordTransGeoKey values from 1 to 27
Posted by Martin Desruisseaux <ma...@geomatys.com>.
Hello Thierry
Le 28/02/2022 à 22:20, Thierry Danard a écrit :
> When I call new MyCRSBuilder(....).getCoordinateOperationFactory().getOperationMethod("GeoTIFF:4");
>
> I get
>
> Exception in thread "main" org.opengis.util.NoSuchIdentifierException:
> No operation method found for name or identifier “GeoTIFF:4”.
"GeoTIFF:4" stands for CT_ObliqueMercator_Laborde, which is not yet
implemented [1]. If replacing by "GeoTIFF:3" (which stands for
CT_ObliqueMercator), then it works.
> I am unsure how the GeoTIFF:xxx operations get plugged in/recognized
> by the code.
>
There is no initializer to call or configuration file, it should be
automatic. All authority codes (GeoTIFF, EPSG, NetCDF, etc.) are
declared in this directory:
core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/provider/
Each file in this directory declares the identifier for a single
operation. For example:
* Mercator1SP.java declares EPSG:9804, GeoTIFF:7 and MapInfo:10.
* Mercator2SP.java declares EPSG:9805, MapInfo:26 and S-57:8.
The declarations appear in the Java code doing the initialization of the
PARAMETERS field in each class. Unfortunately we do not have a this time
a single place listing all GeoTIFF codes (for example) in a tabular
format. Something not too far (should be improved) is [2], which is
generated automatically from the declaration in above-cited Java codes.
> I don't mind the hard-coded logic so much. I'd rather contribute a
> standalone version of the CRSBuilder that essentially takes a map of
> <Short, Object> or <Integer, Object> in its constructor. This would
> make it easier to unit test, and directly reusable in my own use case.
> My use case doesn't need the TIFF file parsing part of the Apache SIS
> storage module, just the interpretation of the GeoTIFF meta data itself.
>
A possible approach could be that, after you got a version that
satisfies your needs, do a "diff" between that version and the current
CRSBuilder so we can see what could be the changes to apply?
Martin
[1] https://issues.apache.org/jira/browse/SIS-222
[2] https://sis.apache.org/tables/CoordinateOperationMethods.html
Re: Support for GeoTIFF ProjCoordTransGeoKey values from 1 to 27
Posted by Thierry Danard <th...@int.com.INVALID>.
I guess that's the part where I am confused.
When I call new
MyCRSBuilder(....).getCoordinateOperationFactory().getOperationMethod("GeoTIFF:4");
I get
Exception in thread "main" org.opengis.util.NoSuchIdentifierException: No
operation method found for name or identifier “GeoTIFF:4”.
at
org.apache.sis.referencing.operation.transform.DefaultMathTransformFactory.getOperationMethod(DefaultMathTransformFactory.java:474)
at
org.apache.sis.referencing.operation.DefaultCoordinateOperationFactory.getOperationMethod(DefaultCoordinateOperationFactory.java:287)
I am unsure how the GeoTIFF:xxx operations get plugged in/recognized by the
code. I am probably doing something wrong, maybe some initializer I should
be calling first? Or a configuration file I am missing?
I created my own MyCRSBuilder class for this test, based on the CRSBuilder
code, but taking only a map of <Integer, Object> in it's constructor .
I don't mind the hard-coded logic so much. I'd rather contribute a
standalone version of the CRSBuilder that essentially takes a map of
<Short, Object> or <Integer, Object> in its constructor. This would make
it easier to unit test, and directly reusable in my own use case. My use
case doesn't need the TIFF file parsing part of the Apache SIS storage
module, just the interpretation of the GeoTIFF meta data itself.
On Mon, Feb 14, 2022 at 5:44 AM Martin Desruisseaux <
martin.desruisseaux@geomatys.com> wrote:
> Hello Thierry
>
> Thanks for your tests and sharing the results. Regarding the support of
> map projections in GeoTIFF reader, the ones listed below should be
> already supported by current reader. Unless there is a bug (in which
> case I would be interested to have more information), it should not be
> necessary to add a switch statement for them.
>
> > case "1": // CT_TransverseMercator
> > case "3": // CT_ObliqueMercator
> > case "4": // CT_ObliqueMercator_Laborde
> > case "7": // CT_Mercator
> > case "8": // CT_LambertConfConic_2SP
> > case "10": // CT_LambertAzimEqualArea
> > case "11": // CT_AlbersEqualArea
> > case "15": // CT_PolarStereographic
> > case "16": // CT_ObliqueStereographic
> > case "17": // CT_Equirectangular
> > case "18": // CT_CassiniSoldner
> > case "26": // CT_NewZealandMapGrid
> > case "27": // CT_TransvMercator_SouthOriented
> > case "20": // CT_MillerCylindrical
> > case "22": // CT_Polyconic
> > case "24": // CT_Sinusoidal
>
>
> However this mapping is not obvious when reading the source code. This
> mapping exists, but is defined outside the GeoTIFF reader. The trick is
> that the following code:
>
> > switch (type) { // based on http://geotiff.maptools.org/proj_list/
> > case "1": // CT_TransverseMercator
> > method =
> getCoordinateOperationFactory().getOperationMethod("EPSG:9807");
> > break;
> > case "3": // CT_ObliqueMercator
> > method =
> getCoordinateOperationFactory().getOperationMethod("EPSG:9815");
> > break;
> > case "4": // CT_ObliqueMercator_Laborde
> > method =
> getCoordinateOperationFactory().getOperationMethod("EPSG:9813");
> > break;
> > case "7": // CT_Mercator
> > method =
> getCoordinateOperationFactory().getOperationMethod("EPSG:9804");
> > etc...
> > }
>
>
> Can also be written as below:
>
> switch (type) { // based on http://geotiff.maptools.org/proj_list/
> case "1": // CT_TransverseMercator
> method =
> getCoordinateOperationFactory().getOperationMethod("GeoTIFF:1");
> break;
> case "3": // CT_ObliqueMercator
> method =
> getCoordinateOperationFactory().getOperationMethod("GeoTIFF:3");
> break;
> case "4": // CT_ObliqueMercator_Laborde
> method =
> getCoordinateOperationFactory().getOperationMethod("GeoTIFF:4");
> break;
> case "7": // CT_Mercator
> method =
> getCoordinateOperationFactory().getOperationMethod("GeoTIFF:7");
> etc...
> }
>
>
> EPSG is one authority, but not the only one. ESRI, GeoTIFF, NetCDF and
> S57 are other authorities known to the SIS referencing module. So
> "GeoTIFF:3" is already considered by SIS as synonymous of "EPSG:9815".
> When using "GeoTIFF" instead of "EPSG" as the authority, the codes in
> the getOperationMethod(…) argument are identical to the switch argument.
> This is why the whole switch statement can be simplified by this single
> line:
>
> method = getCoordinateOperationFactory().getOperationMethod("GeoTIFF:"
> + type);
>
> The advantage of this approach is to make possible to plugin new map
> projections defined by external projects. For example the Geotk project
> defines the "New Zealand Map Grid" projection (that code can not be
> ported to Apache SIS as-is for licensing reason). If Geotk is on the
> classpath, the "GeoTIFF:26" code become automatically recognized. This
> extension capability would be lost if we replaced the current mechanism
> by a switch statement.
>
> An inconvenient of current approach is that it make more difficult to
> see what are the currently supported map projections, as the list does
> not appear explicitly in the code. However it is possible to write a
> small program that automatically generates this list, for example in
> HTML. That generated page could be included in the documentation of
> GeoTIFF module. Is it something you would like to contribute?
>
> Another consequence of this approach is that code like below will not
> really help. It will just replace the "No operation method found for
> name or identifier GeoTIFF:26" message by "No operation method found for
> name or identifier EPSG:9811".
>
> > case "26": // CT_NewZealandMapGrid
> > method =
> getCoordinateOperationFactory().getOperationMethod("EPSG:9811");
> > break;
>
>
> About the following projection:
>
> > case "12": // CT_AzimuthalEquidistant, using Modified Azimuthal
> Equidistant
> I do not have a clear answer about whether it would be correct or not.
> However when looking at the Snyder's book "Map Projection — a working
> manual", the section about Azimuthal Equidistant does not give formulas
> for the general case. Instead the most general formula that I saw was
> the approximation used for Micronesia, which is named "Modified
> Azimuthal Equidistant" by EPSG. My assumption is that a truly generic
> formula for "Azimuthal Equidistant" is not provided by Snyder's book
> because too complex. I also assume that most map projection libraries
> implemented that formula, and because it was given in a section named
> "Azimuthal Equidistant", the libraries used that name.
>
>
> About the non-supported projections:
>
> > case "2": // CT_TransvMercator_Modified_Alaska (not supported by Apache
> SIS)
> > case "5": // CT_ObliqueMercator_Rosenmund (not supported by Apache SIS)
> > case "6": // CT_ObliqueMercator_Spherical (not supported by Apache SIS)
> > case "9": // CT_LambertConfConic_Helmert (not supported by Apache SIS)
> > case "13": // CT_EquidistantConic (not supported by Apache SIS)
> > case "14": // CT_Stereographic (not supported by Apache SIS)
> > case "19": // CT_Gnomonic (not supported by Apache SIS)
> > case "21": // CT_Orthographic (not supported by Apache SIS)
> > case "23": // CT_Robinson (not supported by Apache SIS)
> > case "25": // CT_VanDerGrinten (not supported by Apache SIS)
>
> A problem with Apache SIS as currently designed is that the support of a
> map projection is a "all or nothing" behavior. We have not really
> explored mechanisms for partially implemented map projections, for
> example when only the projection and parameter names are known but the
> actual formula is not yet implemented. If we want to add this mechanism,
> it would actually be a work outside the GeoTIFF module (it would be in
> the referencing module instead). Is it something you would like to
> contribute?
>
> Regards,
>
> Martin
>
>
>
Re: Support for GeoTIFF ProjCoordTransGeoKey values from 1 to 27
Posted by Martin Desruisseaux <ma...@geomatys.com>.
Hello Thierry
Thanks for your tests and sharing the results. Regarding the support of
map projections in GeoTIFF reader, the ones listed below should be
already supported by current reader. Unless there is a bug (in which
case I would be interested to have more information), it should not be
necessary to add a switch statement for them.
> case "1": // CT_TransverseMercator
> case "3": // CT_ObliqueMercator
> case "4": // CT_ObliqueMercator_Laborde
> case "7": // CT_Mercator
> case "8": // CT_LambertConfConic_2SP
> case "10": // CT_LambertAzimEqualArea
> case "11": // CT_AlbersEqualArea
> case "15": // CT_PolarStereographic
> case "16": // CT_ObliqueStereographic
> case "17": // CT_Equirectangular
> case "18": // CT_CassiniSoldner
> case "26": // CT_NewZealandMapGrid
> case "27": // CT_TransvMercator_SouthOriented
> case "20": // CT_MillerCylindrical
> case "22": // CT_Polyconic
> case "24": // CT_Sinusoidal
However this mapping is not obvious when reading the source code. This
mapping exists, but is defined outside the GeoTIFF reader. The trick is
that the following code:
> switch (type) { // based on http://geotiff.maptools.org/proj_list/
> case "1": // CT_TransverseMercator
> method = getCoordinateOperationFactory().getOperationMethod("EPSG:9807");
> break;
> case "3": // CT_ObliqueMercator
> method = getCoordinateOperationFactory().getOperationMethod("EPSG:9815");
> break;
> case "4": // CT_ObliqueMercator_Laborde
> method = getCoordinateOperationFactory().getOperationMethod("EPSG:9813");
> break;
> case "7": // CT_Mercator
> method = getCoordinateOperationFactory().getOperationMethod("EPSG:9804");
> etc...
> }
Can also be written as below:
switch (type) { // based on http://geotiff.maptools.org/proj_list/
case "1": // CT_TransverseMercator
method = getCoordinateOperationFactory().getOperationMethod("GeoTIFF:1");
break;
case "3": // CT_ObliqueMercator
method = getCoordinateOperationFactory().getOperationMethod("GeoTIFF:3");
break;
case "4": // CT_ObliqueMercator_Laborde
method = getCoordinateOperationFactory().getOperationMethod("GeoTIFF:4");
break;
case "7": // CT_Mercator
method = getCoordinateOperationFactory().getOperationMethod("GeoTIFF:7");
etc...
}
EPSG is one authority, but not the only one. ESRI, GeoTIFF, NetCDF and
S57 are other authorities known to the SIS referencing module. So
"GeoTIFF:3" is already considered by SIS as synonymous of "EPSG:9815".
When using "GeoTIFF" instead of "EPSG" as the authority, the codes in
the getOperationMethod(…) argument are identical to the switch argument.
This is why the whole switch statement can be simplified by this single
line:
method = getCoordinateOperationFactory().getOperationMethod("GeoTIFF:" + type);
The advantage of this approach is to make possible to plugin new map
projections defined by external projects. For example the Geotk project
defines the "New Zealand Map Grid" projection (that code can not be
ported to Apache SIS as-is for licensing reason). If Geotk is on the
classpath, the "GeoTIFF:26" code become automatically recognized. This
extension capability would be lost if we replaced the current mechanism
by a switch statement.
An inconvenient of current approach is that it make more difficult to
see what are the currently supported map projections, as the list does
not appear explicitly in the code. However it is possible to write a
small program that automatically generates this list, for example in
HTML. That generated page could be included in the documentation of
GeoTIFF module. Is it something you would like to contribute?
Another consequence of this approach is that code like below will not
really help. It will just replace the "No operation method found for
name or identifier GeoTIFF:26" message by "No operation method found for
name or identifier EPSG:9811".
> case "26": // CT_NewZealandMapGrid
> method = getCoordinateOperationFactory().getOperationMethod("EPSG:9811");
> break;
About the following projection:
> case "12": // CT_AzimuthalEquidistant, using Modified Azimuthal Equidistant
I do not have a clear answer about whether it would be correct or not.
However when looking at the Snyder's book "Map Projection — a working
manual", the section about Azimuthal Equidistant does not give formulas
for the general case. Instead the most general formula that I saw was
the approximation used for Micronesia, which is named "Modified
Azimuthal Equidistant" by EPSG. My assumption is that a truly generic
formula for "Azimuthal Equidistant" is not provided by Snyder's book
because too complex. I also assume that most map projection libraries
implemented that formula, and because it was given in a section named
"Azimuthal Equidistant", the libraries used that name.
About the non-supported projections:
> case "2": // CT_TransvMercator_Modified_Alaska (not supported by Apache SIS)
> case "5": // CT_ObliqueMercator_Rosenmund (not supported by Apache SIS)
> case "6": // CT_ObliqueMercator_Spherical (not supported by Apache SIS)
> case "9": // CT_LambertConfConic_Helmert (not supported by Apache SIS)
> case "13": // CT_EquidistantConic (not supported by Apache SIS)
> case "14": // CT_Stereographic (not supported by Apache SIS)
> case "19": // CT_Gnomonic (not supported by Apache SIS)
> case "21": // CT_Orthographic (not supported by Apache SIS)
> case "23": // CT_Robinson (not supported by Apache SIS)
> case "25": // CT_VanDerGrinten (not supported by Apache SIS)
A problem with Apache SIS as currently designed is that the support of a
map projection is a "all or nothing" behavior. We have not really
explored mechanisms for partially implemented map projections, for
example when only the projection and parameter names are known but the
actual formula is not yet implemented. If we want to add this mechanism,
it would actually be a work outside the GeoTIFF module (it would be in
the referencing module instead). Is it something you would like to
contribute?
Regards,
Martin