sis-commits mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From desruisse...@apache.org
Subject [sis] branch geoapi-4.0 updated: Add inference of "grid to CRS" transform from netCDF attributes associated to "grid_mapping" variable.
Date Fri, 26 Apr 2019 10:36:35 GMT
This is an automated email from the ASF dual-hosted git repository.

desruisseaux pushed a commit to branch geoapi-4.0
in repository https://gitbox.apache.org/repos/asf/sis.git


The following commit(s) were added to refs/heads/geoapi-4.0 by this push:
     new 53ade5d  Add inference of "grid to CRS" transform from netCDF attributes associated
to "grid_mapping" variable.
53ade5d is described below

commit 53ade5da47e913c1edd84e615ba811b695f4a3f9
Author: Martin Desruisseaux <martin.desruisseaux@geomatys.com>
AuthorDate: Fri Apr 26 12:36:06 2019 +0200

    Add inference of "grid to CRS" transform from netCDF attributes associated to "grid_mapping"
variable.
---
 .../internal/referencing/provider/Sinusoidal.java  |  7 +-
 .../apache/sis/internal/earth/netcdf/GCOM_C.java   | 94 +++++++++++++++++++---
 .../org/apache/sis/internal/netcdf/Convention.java | 77 +++++++++++++-----
 .../apache/sis/internal/netcdf/GridMapping.java    |  7 +-
 4 files changed, 153 insertions(+), 32 deletions(-)

diff --git a/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/provider/Sinusoidal.java
b/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/provider/Sinusoidal.java
index 7c00ea9..3335093 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/provider/Sinusoidal.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/provider/Sinusoidal.java
@@ -48,6 +48,11 @@ public final class Sinusoidal extends MapProjection {
     private static final long serialVersionUID = -3236247448683326299L;
 
     /**
+     * Name of this projection.
+     */
+    public static final String NAME = "Sinusoidal";
+
+    /**
      * The operation parameter descriptor for the <cite>Longitude of projection center</cite>
(λ₀) parameter value.
      * Valid values range is [-180 … 180]° and default value is 0°.
      */
@@ -71,7 +76,7 @@ public final class Sinusoidal extends MapProjection {
     private static final ParameterDescriptorGroup PARAMETERS;
     static {
         PARAMETERS = builder().setCodeSpace(Citations.OGC, Constants.OGC)
-                .addName      ("Sinusoidal")
+                .addName      (NAME)
                 .addName      ("Sanson-Flamsteed")
                 .addName      (Citations.GEOTIFF,  "CT_Sinusoidal")
                 .addIdentifier(Citations.GEOTIFF,  "24")
diff --git a/storage/sis-earth-observation/src/main/java/org/apache/sis/internal/earth/netcdf/GCOM_C.java
b/storage/sis-earth-observation/src/main/java/org/apache/sis/internal/earth/netcdf/GCOM_C.java
index a541688..0a88e7e 100644
--- a/storage/sis-earth-observation/src/main/java/org/apache/sis/internal/earth/netcdf/GCOM_C.java
+++ b/storage/sis-earth-observation/src/main/java/org/apache/sis/internal/earth/netcdf/GCOM_C.java
@@ -22,6 +22,9 @@ import java.util.HashMap;
 import java.util.Collections;
 import java.util.LinkedHashSet;
 import java.util.regex.Pattern;
+import org.opengis.referencing.crs.ProjectedCRS;
+import org.opengis.referencing.operation.MathTransform;
+import org.opengis.referencing.operation.TransformException;
 import org.apache.sis.storage.netcdf.AttributeNames;
 import org.apache.sis.internal.netcdf.Convention;
 import org.apache.sis.internal.netcdf.Decoder;
@@ -29,9 +32,14 @@ import org.apache.sis.internal.netcdf.Variable;
 import org.apache.sis.internal.netcdf.VariableRole;
 import org.apache.sis.internal.netcdf.Linearizer;
 import org.apache.sis.internal.netcdf.Node;
+import org.apache.sis.internal.referencing.provider.Sinusoidal;
+import org.apache.sis.internal.referencing.provider.Equirectangular;
+import org.apache.sis.internal.referencing.provider.PolarStereographicA;
 import org.apache.sis.referencing.operation.transform.TransferFunction;
-import org.apache.sis.measure.NumberRange;
+import org.apache.sis.referencing.operation.transform.MathTransforms;
+import org.apache.sis.referencing.operation.matrix.Matrix3;
 import org.apache.sis.referencing.CommonCRS;
+import org.apache.sis.measure.NumberRange;
 
 
 /**
@@ -369,16 +377,80 @@ public final class GCOM_C extends Convention {
      */
     @Override
     public Map<String,Object> projection(final Node node) {
-        String method = node.getAttributeAsString("Image_projection");
-        if (method != null) {
-            if (method.matches("EQA\\b.*")) {
-                method = "Sinusoidal";
-                final Map<String,Object> definition = new HashMap<>(4);
-                definition.put(BASE_CRS, CommonCRS.WGS84.geographic());
-                definition.put("grid_mapping_name", method);
-                return definition;
-            }
+        final String name = node.getAttributeAsString("Image_projection");
+        if (name == null) {
+            return super.projection(node);
+        }
+        final String method;
+        final int s = name.indexOf(' ');
+        final String code = (s >= 0) ? name.substring(0, s) : name;
+        if (code.equalsIgnoreCase("EQA")) {
+            method = Sinusoidal.NAME;
+        } else if (code.equalsIgnoreCase("EQR")) {
+            method = Equirectangular.NAME;
+        } else if (code.equalsIgnoreCase("PS")) {
+            method = PolarStereographicA.NAME;
+        } else {
+            return super.projection(node);
+        }
+        final Map<String,Object> definition = new HashMap<>(4);
+        definition.put(BASE_CRS, CommonCRS.SPHERE.geographic());
+        definition.put("grid_mapping_name", method);
+        return definition;
+    }
+
+    /**
+     * The attributes storing values of the 4 corners in degrees with (latitude, longitude)
axis order.
+     * This is used by {@link #projection(Node)} for inferring a "grid to CRS" transform.
+     */
+    private static final String[] CORNERS = {
+        "Upper_left_latitude",
+        "Upper_left_longitude",
+        "Upper_right_latitude",
+        "Upper_right_longitude",
+        "Lower_left_latitude",
+        "Lower_left_longitude",
+        "Lower_right_latitude",
+        "Lower_right_longitude"
+    };
+
+    /**
+     * Returns the <cite>grid to CRS</cite> transform for the given node.
+     * This method is invoked after call to {@link #projection(Node)} resulted in creation
of a projected CRS.
+     * The {@linkplain ProjectedCRS#getBaseCRS() base CRS} shall have (latitude, longitude)
axes in degrees.
+     *
+     * @param  node  the same node than the one given to {@link #projection(Node)}.
+     * @param  crs   the projected coordinate reference system created from the information
given by {@code node}.
+     * @return the "grid corner to CRS" transform, or {@code null} if none or unknown.
+     * @throws TransformException if a coordinate operation was required but failed.
+     */
+    @Override
+    public MathTransform gridToCRS(final Node node, final ProjectedCRS crs) throws TransformException
{
+        final double[] corners = new double[CORNERS.length];
+        for (int i=0; i<corners.length; i++) {
+            corners[i] = node.getAttributeAsNumber(CORNERS[i]);
+        }
+        crs.getConversionFromBase().getMathTransform().transform(corners, 0, corners, 0,
corners.length / 2);
+        /*
+         * Compute spans of data (typically in metres) as the average of the spans on both
sides
+         * (width as length of top and bottom edges, height as length of left and right edges).
+         * This code assumes (easting, northing) axes — this is currently not verified.
+         */
+        double sx = ((corners[2] - corners[0]) + (corners[6] - corners[4])) / 2;
+        double sy = ((corners[1] - corners[5]) + (corners[3] - corners[7])) / 2;
+        /*
+         * Transform the spans into pixel sizes (resolution), then build the transform.
+         */
+        sx /= (node.getAttributeAsNumber("Number_of_pixels") - 1);
+        sy /= (node.getAttributeAsNumber("Number_of_lines")  - 1);
+        if (Double.isFinite(sx) && Double.isFinite(sy)) {
+            final Matrix3 m = new Matrix3();
+            m.m00 =  sx;
+            m.m11 = -sy;
+            m.m02 = corners[0];
+            m.m12 = corners[1];
+            return MathTransforms.linear(m);
         }
-        return super.projection(node);
+        return super.gridToCRS(node, crs);
     }
 }
diff --git a/storage/sis-netcdf/src/main/java/org/apache/sis/internal/netcdf/Convention.java
b/storage/sis-netcdf/src/main/java/org/apache/sis/internal/netcdf/Convention.java
index 55ffd26..b1aaeae 100644
--- a/storage/sis-netcdf/src/main/java/org/apache/sis/internal/netcdf/Convention.java
+++ b/storage/sis-netcdf/src/main/java/org/apache/sis/internal/netcdf/Convention.java
@@ -24,12 +24,17 @@ import java.util.Collections;
 import java.util.LinkedHashMap;
 import java.util.Locale;
 import java.awt.image.DataBuffer;
+import org.opengis.referencing.crs.ProjectedCRS;
 import org.opengis.referencing.crs.GeographicCRS;
 import org.opengis.referencing.operation.MathTransform;
-import org.apache.sis.internal.referencing.LazySet;
+import org.opengis.referencing.operation.TransformException;
 import org.apache.sis.referencing.operation.transform.TransferFunction;
+import org.apache.sis.referencing.datum.BursaWolfParameters;
+import org.apache.sis.referencing.CommonCRS;
+import org.apache.sis.internal.referencing.LazySet;
 import org.apache.sis.measure.MeasurementRange;
 import org.apache.sis.measure.NumberRange;
+import org.apache.sis.math.Vector;
 import org.apache.sis.util.Numbers;
 import org.apache.sis.util.resources.Errors;
 import ucar.nc2.constants.CDM;
@@ -474,10 +479,12 @@ public class Convention {
     /**
      * Returns the name of nodes (variables or groups) that may define the map projection
parameters.
      * The variables or groups will be inspected in the order they are declared in the returned
set.
-     * For each string in the set, {@link Decoder#findNode(String)} is invoked.
+     * For each string in the set, {@link Decoder#findNode(String)} is invoked and the return
value
+     * (if non-null) is given to {@link #projection(Node)}.
      *
      * <p>The default implementation returns the value of {@link CF#GRID_MAPPING},
or an empty set
-     * if the given variable does not contain that attribute.</p>
+     * if the given variable does not contain that attribute. Subclasses may override for
example
+     * if grid mapping information are hard-coded in a particular node for a specific product.</p>
      *
      * @param  data  the variable for which to get the grid mapping node.
      * @return name of nodes that may contain the grid mapping, or an empty set if none.
@@ -493,15 +500,10 @@ public class Convention {
     protected static final String BASE_CRS = "base_crs";
 
     /**
-     * Key associated to {@link MathTransform} value in the map returned by {@link #projection(Node)}.
-     */
-    protected static final String GRID_TO_CRS = "grid_to_crs";
-
-    /**
-     * Returns the map projection defined by the given node.  The given {@code node} is typically
the variable named by a
-     * {@value CF#GRID_MAPPING} attribute (this is the preferred, CF-compliant approach),
or if no grid mapping attribute
-     * is found {@code node} may be directly the data variable (non CF-compliant, but found
in practice).
-     * If non-null, the returned map contains the following information:
+     * Returns the map projection defined by the given node. The given {@code node} argument
is one of the nodes
+     * named by {@link #gridMapping(Variable)} (typically a variable named by {@value CF#GRID_MAPPING}
attribute),
+     * or if no grid mapping attribute is found {@code node} may be directly the data variable
(non CF-compliant,
+     * but found in practice). If non-null, the returned map contains the following information:
      *
      * <table class="sis">
      *   <caption>Content of the returned map</caption>
@@ -531,10 +533,10 @@ public class Convention {
      *     <td>Map projection parameter values</td>
      *     <td>Attributes found on grid mapping variable.</td>
      *   </tr><tr>
-     *     <td>{@value #GRID_TO_CRS}</td>
-     *     <td>{@link MathTransform}</td>
-     *     <td>Conversion from pixel indices to CRS.</td>
-     *     <td>None (not from CF-convention).</td>
+     *     <td>{@code "towgs84"}</td>
+     *     <td>{@link BursaWolfParameters}</td>
+     *     <td>Datum shift information.</td>
+     *     <td>Attributes found on grid mapping variable.</td>
      *   </tr>
      * </table>
      *
@@ -543,6 +545,8 @@ public class Convention {
      *
      * @param  node  the {@value CF#GRID_MAPPING} variable (preferred) or the data variable
(as a fallback) from which to read attributes.
      * @return the map projection definition, or {@code null} if none.
+     *
+     * @see <a href="http://cfconventions.org/cf-conventions/cf-conventions.html#grid-mappings-and-projections">CF-conventions</a>
      */
     public Map<String,Object> projection(final Node node) {
         final String method = node.getAttributeAsString(CF.GRID_MAPPING_NAME);
@@ -561,8 +565,27 @@ public class Convention {
                 break;
             } else switch (ln) {
                 case CF.GRID_MAPPING_NAME: continue;        // Already stored.
-                case "towgs84":            continue;        // TODO: process in a future
version.
-                case "crs_wkt":            continue;        // TODO: process in a future
version.
+                case "towgs84": {
+                    /*
+                     * Conversion to WGS 84 datum may be specified as Bursa-Wolf parameters.
Encoding this information
+                     * with the CRS is deprecated (the hard-coded WGS84 target datum is not
always suitable) but still
+                     * a common practice as of 2019. We require at least the 3 translation
parameters.
+                     */
+                    final Object[] values = node.getAttributeValues(name, true);
+                    if (values.length < 3) continue;
+                    final BursaWolfParameters bp = new BursaWolfParameters(CommonCRS.WGS84.datum(),
null);
+                    bp.setValues(Vector.create(values, false).doubleValues());
+                    value = bp;
+                    break;
+                }
+                case "crs_wkt": {
+                    /*
+                     * CF-Convention said that even if a WKT definition is provided, other
attributes shall be present
+                     * and have precedence over the WKT definition. Consequently purpose
of WKT in netCDF files is not
+                     * obvious (except for CompoundCRS). We ignore them for now.
+                     */
+                    continue;
+                }
                 default: {
                     final double n = node.getAttributeAsNumber(name);
                     if (Double.isNaN(n)) continue;
@@ -576,4 +599,22 @@ public class Convention {
         }
         return definition;
     }
+
+    /**
+     * Returns the <cite>grid to CRS</cite> transform for the given node. This
method is invoked after call
+     * to {@link #projection(Node)} method resulted in creation of a projected coordinate
reference system.
+     * The {@linkplain ProjectedCRS#getBaseCRS() base CRS} shall have (latitude, longitude)
axes in degrees.
+     * The projected CRS axes may have any order and units.
+     * The returned transform, if non-null, shall map cell corners.
+     *
+     * <p>The default implementation returns {@code null}.</p>
+     *
+     * @param  node  the same node than the one given to {@link #projection(Node)}.
+     * @param  crs   the projected coordinate reference system created from the information
given by {@code node}.
+     * @return the "grid corner to CRS" transform, or {@code null} if none or unknown.
+     * @throws TransformException if a coordinate operation was required but failed.
+     */
+    public MathTransform gridToCRS(final Node node, final ProjectedCRS crs) throws TransformException
{
+        return null;
+    }
 }
diff --git a/storage/sis-netcdf/src/main/java/org/apache/sis/internal/netcdf/GridMapping.java
b/storage/sis-netcdf/src/main/java/org/apache/sis/internal/netcdf/GridMapping.java
index eb9628e..7778b2a 100644
--- a/storage/sis-netcdf/src/main/java/org/apache/sis/internal/netcdf/GridMapping.java
+++ b/storage/sis-netcdf/src/main/java/org/apache/sis/internal/netcdf/GridMapping.java
@@ -29,6 +29,7 @@ import org.opengis.referencing.cs.CoordinateSystem;
 import org.opengis.referencing.crs.ProjectedCRS;
 import org.opengis.referencing.crs.GeographicCRS;
 import org.opengis.referencing.crs.CoordinateReferenceSystem;
+import org.opengis.referencing.operation.TransformException;
 import org.opengis.referencing.operation.CoordinateOperationFactory;
 import org.opengis.referencing.operation.OperationMethod;
 import org.opengis.referencing.operation.MathTransform;
@@ -161,8 +162,10 @@ final class GridMapping {
             final GeographicCRS baseCRS = (GeographicCRS) definition.get(Convention.BASE_CRS);
             final CartesianCS cs = CommonCRS.WGS84.universal(0,0).getCoordinateSystem();
                    // TODO
             final ProjectedCRS crs = node.decoder.getCRSFactory().createProjectedCRS(name,
baseCRS, conversion, cs);
-            return new GridMapping(crs, null, false);
-        } catch (ClassCastException | IllegalArgumentException | FactoryException e) {
+
+            final MathTransform gridToCRS = node.decoder.convention().gridToCRS(node, crs);
+            return new GridMapping(crs, gridToCRS, false);
+        } catch (ClassCastException | IllegalArgumentException | FactoryException | TransformException
e) {
             canNotCreate(node, Resources.Keys.CanNotCreateCRS_3, e);
         }
         return null;


Mime
View raw message