Fix projection units

cdm/core/src/main/java/ucar/nc2/dataset/transform/Geostationary.java

@@ -78,21 +78,6 @@ public TransformType getTransformType() {
     return TransformType.Projection;
   }
 
-  private double getScaleFactor(String geoCoordinateUnits) {
-    // default value of -1.0 interpreted as no scaling in the class
-    // ucar.unidata.geoloc.projection.sat.Geostationary
-    double scaleFactor = defaultScaleFactor;
-    String neededMapCoordinateUnit = "radian";
-
-    if (SimpleUnit.isCompatible(geoCoordinateUnits, neededMapCoordinateUnit)) {
-      scaleFactor = SimpleUnit.getConversionFactor(geoCoordinateUnits, neededMapCoordinateUnit);
-    }
-
-    logger.debug("geoCoordinateUnits {}, scaleFactor {}", geoCoordinateUnits, scaleFactor);
-
-    return scaleFactor;
-  }
-
   public ProjectionCT makeCoordinateTransform(AttributeContainer ctv, String geoCoordinateUnits) {
     readStandardParams(ctv, geoCoordinateUnits);
 
@@ -145,13 +130,8 @@ public ProjectionCT makeCoordinateTransform(AttributeContainer ctv, String geoCo
       isSweepX = fixed_angle.equals("y");
     }
 
-    // scales less than zero indicate no scaling of axis (i.e. map coords have units of radians)
-    double geoCoordinateScaleFactor;
-
-    geoCoordinateScaleFactor = getScaleFactor(geoCoordinateUnits);
-
     ProjectionImpl proj = new ucar.unidata.geoloc.projection.sat.Geostationary(subLonDegrees, perspective_point_height,
-        semi_minor_axis, semi_major_axis, inv_flattening, isSweepX, geoCoordinateScaleFactor);
+        semi_minor_axis, semi_major_axis, inv_flattening, isSweepX);
 
     return new ProjectionCT(ctv.getName(), "FGDC", proj);
   }

cdm/core/src/main/java/ucar/nc2/ft2/coverage/HorizCoordSys.java

@@ -10,6 +10,7 @@
 import org.slf4j.Logger;
 import org.slf4j.LoggerFactory;
 import ucar.ma2.*;
+import ucar.nc2.units.SimpleUnit;
 import ucar.nc2.util.Optional;
 import ucar.unidata.geoloc.*;
 import javax.annotation.concurrent.Immutable;
@@ -58,6 +59,8 @@ public static HorizCoordSys factory(CoverageCoordAxis1D xAxis, CoverageCoordAxis
   private final boolean isProjection;
   private final boolean isLatLon1D;
   private boolean isLatLon2D; // isProjection and isLatLon2D may both be "true".
+  // scale factor for x, y axis if they have different units than the projection's default units
+  private final double coordinateConversionFactor;
 
   protected HorizCoordSys(CoverageCoordAxis1D xAxis, CoverageCoordAxis1D yAxis, CoverageCoordAxis latAxis,
       CoverageCoordAxis lonAxis, CoverageTransform transform) {
@@ -68,6 +71,7 @@ protected HorizCoordSys(CoverageCoordAxis1D xAxis, CoverageCoordAxis1D yAxis, Co
     this.isLatLon1D = latAxis instanceof CoverageCoordAxis1D && lonAxis instanceof CoverageCoordAxis1D;
     this.isLatLon2D = latAxis instanceof LatLonAxis2D && lonAxis instanceof LatLonAxis2D;
     assert isProjection || isLatLon1D || isLatLon2D : "missing horiz coordinates (x,y,projection or lat,lon)";
+    coordinateConversionFactor = getCoordinateConversionFactor();
 
     if (isProjection && isLatLon2D) {
       boolean ok = true;
@@ -163,16 +167,14 @@ public Optional<HorizCoordSys> subset(SubsetParams params) {
 
           // we have to transform latlon to projection coordinates
           ProjectionImpl proj = transform.getProjection();
-          final ProjectionPoint projectionPointInKm = proj.latLonToProj(latlon);
-          final double xInCorrectUnits = convertFromKm(projectionPointInKm.getX(), xAxis.units, xAxis.name);
-          optb = xhelper.subsetContaining(xInCorrectUnits);
+          final ProjectionPoint projectionRectInDefaultUnits = proj.latLonToProj(latlon);
+          optb = xhelper.subsetContaining(convertFromDefaultUnits(projectionRectInDefaultUnits.getX()));
           if (optb.isPresent())
             xaxisSubset = new CoverageCoordAxis1D(optb.get());
           else
             errMessages.format("xaxis: %s;%n", optb.getErrorMessage());
 
-          final double yInCorrectUnits = convertFromKm(projectionPointInKm.getY(), yAxis.units, yAxis.name);
-          optb = yhelper.subsetContaining(yInCorrectUnits);
+          optb = yhelper.subsetContaining(convertFromDefaultUnits(projectionRectInDefaultUnits.getY()));
           if (optb.isPresent())
             yaxisSubset = new CoverageCoordAxis1D(optb.get());
           else
@@ -233,17 +235,18 @@ public Optional<HorizCoordSys> subset(SubsetParams params) {
           if (isProjection) {
             // we have to transform latlon to projection coordinates
             ProjectionImpl proj = transform.getProjection();
-            final ProjectionRect projectionRectInKm = proj.latLonToProjBB(llbb); // allow projection to override
-            final double xMinInCorrectUnits = convertFromKm(projectionRectInKm.getMinX(), xAxis.units, xAxis.name);
-            final double xMaxInCorrectUnits = convertFromKm(projectionRectInKm.getMaxX(), xAxis.units, xAxis.name);
+            final ProjectionRect projectionRectInDefaultUnits = proj.latLonToProjBB(llbb); // allow projection to
+                                                                                           // override
+            final double xMinInCorrectUnits = convertFromDefaultUnits(projectionRectInDefaultUnits.getMinX());
+            final double xMaxInCorrectUnits = convertFromDefaultUnits(projectionRectInDefaultUnits.getMaxX());
             opt = xAxis.subset(xMinInCorrectUnits, xMaxInCorrectUnits, horizStride);
             if (opt.isPresent())
               xaxisSubset = (CoverageCoordAxis1D) opt.get();
             else
               errMessages.format("xaxis: %s;%n", opt.getErrorMessage());
 
-            final double yMinInCorrectUnits = convertFromKm(projectionRectInKm.getMinY(), yAxis.units, yAxis.name);
-            final double yMaxInCorrectUnits = convertFromKm(projectionRectInKm.getMaxY(), yAxis.units, yAxis.name);
+            final double yMinInCorrectUnits = convertFromDefaultUnits(projectionRectInDefaultUnits.getMinY());
+            final double yMaxInCorrectUnits = convertFromDefaultUnits(projectionRectInDefaultUnits.getMaxY());
             opt = yAxis.subset(yMinInCorrectUnits, yMaxInCorrectUnits, horizStride);
             if (opt.isPresent())
               yaxisSubset = (CoverageCoordAxis1D) opt.get();
@@ -319,9 +322,7 @@ public LatLonPoint getLatLon(int yindex, int xindex) {
       double x = xAxis.getCoordMidpoint(xindex);
       double y = yAxis.getCoordMidpoint(yindex);
       ProjectionImpl proj = transform.getProjection();
-      final double xInKm = convertToKm(x, xAxis.units, xAxis.name);
-      final double yInKm = convertToKm(y, yAxis.units, yAxis.name);
-      return proj.projToLatLon(xInKm, yInKm);
+      return proj.projToLatLon(convertToDefaultUnits(x), convertToDefaultUnits(y));
     } else {
       double lat = latAxis.getCoordMidpoint(yindex);
       double lon = lonAxis.getCoordMidpoint(xindex);
@@ -662,8 +663,7 @@ private List<LatLonPoint> calcLatLonBoundaryPointsFromProjection(int maxPointsIn
 
     for (ProjectionPoint projPoint : projPoints) {
       final ProjectionPoint projPointInKm =
-          ProjectionPoint.create(convertToKm(projPoint.getX(), xAxis.units, xAxis.name),
-              convertToKm(projPoint.getY(), yAxis.units, yAxis.name));
+          ProjectionPoint.create(convertToDefaultUnits(projPoint.getX()), convertToDefaultUnits(projPoint.getY()));
 
       final LatLonPoint latLonPoint = transform.getProjection().projToLatLon(projPointInKm);
       if (!Double.isNaN(latLonPoint.getLatitude()) && !Double.isNaN(latLonPoint.getLongitude())) {
@@ -674,27 +674,22 @@ private List<LatLonPoint> calcLatLonBoundaryPointsFromProjection(int maxPointsIn
     return latLonPoints;
   }
 
-  // TODO is there a better place to handle units?
-  // Some projections are actually just rotations (RotatedPole)
-  // so the "projection" coordinates have units "degrees" and don't need to be converted
-  private static double convertToKm(double coordinate, String unit, String axisName) {
-    if (unit.equals("km") || unit.equals("kilometers")) {
-      return coordinate;
-    } else if (unit.equals("m") || unit.equals("meters")) {
-      return 0.001 * coordinate;
-    } else {
-      return coordinate;
+  private double getCoordinateConversionFactor() {
+    if (!isProjection) {
+      return 1.0;
     }
+
+    final String defaultUnits = transform.getProjection().getDefaultUnits();
+    final String unit = xAxis.getUnits();
+    return SimpleUnit.isCompatible(unit, defaultUnits) ? SimpleUnit.getConversionFactor(unit, defaultUnits) : 1.0;
   }
 
-  private static double convertFromKm(double coordinateInKm, String desiredUnit, String axisName) {
-    if (desiredUnit.equals("km") || desiredUnit.equals("kilometers")) {
-      return coordinateInKm;
-    } else if (desiredUnit.equals("m") || desiredUnit.equals("meters")) {
-      return 1000 * coordinateInKm;
-    } else {
-      return coordinateInKm;
-    }
+  private double convertToDefaultUnits(double coordinate) {
+    return coordinate * coordinateConversionFactor;
+  }
+
+  private double convertFromDefaultUnits(double coordinateInDefaultUnit) {
+    return coordinateInDefaultUnit / coordinateConversionFactor;
   }
 
   private List<LatLonPoint> calcLatLon2DBoundaryPoints(int maxPointsInYEdge, int maxPointsInXEdge) {

cdm/core/src/main/java/ucar/unidata/geoloc/ProjectionImpl.java

@@ -73,6 +73,11 @@ public abstract class ProjectionImpl implements Projection, java.io.Serializable
    */
   protected String name; // LOOK should be final, IDV needs setName()
 
+  /**
+   * name of the default units for this projection
+   */
+  protected String defaultUnits;
+
   /**
    * flag for latlon
    */
@@ -96,6 +101,12 @@ public abstract class ProjectionImpl implements Projection, java.io.Serializable
    */
   public abstract ProjectionImpl constructCopy();
 
+  protected ProjectionImpl(String name, String defaultUnits, boolean isLatLon) {
+    this.name = name;
+    this.defaultUnits = defaultUnits;
+    this.isLatLon = isLatLon;
+  }
+
   protected ProjectionImpl(String name, boolean isLatLon) {
     this.name = name;
     this.isLatLon = isLatLon;
@@ -209,6 +220,15 @@ public void setName(String name) {
     this.name = name;
   }
 
+  /**
+   * Get the name of the default units for this projection
+   *
+   * @return the name of the default units
+   */
+  public String getDefaultUnits() {
+    return defaultUnits;
+  }
+
   /**
    * Get parameters as list of ucar.unidata.util.Parameter
    *

cdm/core/src/main/java/ucar/unidata/geoloc/projection/LambertConformal.java

@@ -23,6 +23,9 @@
  */
 
 public class LambertConformal extends ProjectionImpl {
+  private static final String NAME = "LambertConformal";
+  private static final String DEFAULT_UNITS = "km";
+
   private final double earth_radius;
   private double lat0, lon0; // lat/lon in radians
   private double par1, par2; // standard parallel 1 and 2 degrees
@@ -100,7 +103,7 @@ public LambertConformal(double lat0, double lon0, double par1, double par2, doub
   public LambertConformal(double lat0, double lon0, double par1, double par2, double false_easting,
       double false_northing, double earth_radius) {
 
-    super("LambertConformal", false);
+    super(NAME, DEFAULT_UNITS, false);
 
     this._lat0 = lat0;
     this._lon0 = lon0;

cdm/core/src/main/java/ucar/unidata/geoloc/projection/sat/Geostationary.java

@@ -73,6 +73,9 @@ public class Geostationary extends ProjectionImpl {
   private static final org.slf4j.Logger logger = org.slf4j.LoggerFactory.getLogger(Geostationary.class);
 
   private static final String NAME = CF.GEOSTATIONARY;
+  private static final String DEFAULT_UNITS = "radians";
+
+  // Remove in v6.x
   private boolean scaleGeoCoordinate;
   private double geoCoordinateScaleFactor = Double.MIN_VALUE;
 
@@ -81,13 +84,18 @@ public class Geostationary extends ProjectionImpl {
   public Geostationary(double subLonDegrees, double perspective_point_height, double semi_minor_axis,
       double semi_major_axis, double inv_flattening, boolean isSweepX) {
 
-    // scale factors (last two doubles in the sig) less than zero indicate no scaling of map x, y coordinates
+    // scale factors less than zero indicate no scaling of map x, y coordinates
     this(subLonDegrees, perspective_point_height, semi_minor_axis, semi_major_axis, inv_flattening, isSweepX, -1.0);
   }
 
+  /**
+   * @deprecated Remove in v6.x.
+   *             Use constructor without geoCoordinateScaleFactor as units are handled outside of projection classes
+   */
+  @Deprecated
   public Geostationary(double subLonDegrees, double perspective_point_height, double semi_minor_axis,
       double semi_major_axis, double inv_flattening, boolean isSweepX, double geoCoordinateScaleFactor) {
-    super(NAME, false);
+    super(NAME, DEFAULT_UNITS, false);
 
     String sweepAngleAxis = "y";
     if (isSweepX) {
@@ -112,19 +120,19 @@ public Geostationary(double subLonDegrees, double perspective_point_height, doub
   }
 
   public Geostationary() {
-    super(NAME, false);
+    super(NAME, DEFAULT_UNITS, false);
     navigation = new GEOSTransform();
     makePP();
   }
 
   public Geostationary(double subLonDegrees) {
-    super(NAME, false);
+    super(NAME, DEFAULT_UNITS, false);
     navigation = new GEOSTransform(subLonDegrees, GEOSTransform.GOES);
     makePP();
   }
 
   public Geostationary(double subLonDegrees, boolean isSweepX) {
-    super(NAME, false);
+    super(NAME, DEFAULT_UNITS, false);
 
     String sweepAngleAxis = "y";
     if (isSweepX) {
@@ -137,8 +145,13 @@ public Geostationary(double subLonDegrees, boolean isSweepX) {
     makePP();
   }
 
+  /**
+   * @deprecated Remove in v6.x.
+   *             Use constructor without geoCoordinateScaleFactor as units are handled outside of projection classes
+   */
+  @Deprecated
   public Geostationary(double subLonDegrees, String sweepAngleAxis, double geoCoordinateScaleFactor) {
-    super(NAME, false);
+    super(NAME, DEFAULT_UNITS, false);
 
     String scanGeometry = GEOSTransform.sweepAngleAxisToScanGeom(sweepAngleAxis);
 
@@ -164,6 +177,11 @@ private void makePP() {
     addParameter(CF.SEMI_MINOR_AXIS, navigation.r_pol * 1000.0);
   }
 
+  /**
+   * @deprecated Remove in v6.x.
+   *             Units are handled outside of projection classes
+   */
+  @Deprecated
   private boolean isGeoCoordinateScaled() {
     return scaleGeoCoordinate && geoCoordinateScaleFactor > Double.MIN_VALUE;
   }

cdm/core/src/test/java/ucar/nc2/ft2/coverage/TestHorizCoordSys.java

@@ -17,12 +17,15 @@
 import ucar.nc2.constants.AxisType;
 import ucar.nc2.constants.CF;
 import ucar.nc2.ft2.coverage.CoverageCoordAxis.Spacing;
+import ucar.nc2.util.Optional;
+import ucar.unidata.geoloc.LatLonPoint;
 import ucar.unidata.geoloc.LatLonPointNoNormalize;
 import ucar.unidata.geoloc.ProjectionPoint;
 
 public class TestHorizCoordSys {
   private static final Logger logger = LoggerFactory.getLogger(MethodHandles.lookup().lookupClass());
   private static final Locale DEFAULT_LOCALE = Locale.getDefault();
+  private static final double TOLERANCE = 1.0e-8;
 
   @After
   public void resetLocale() {
@@ -38,14 +41,7 @@ public void shouldRemoveNansWhenComputingLatLon() {
     final CoverageCoordAxis1D xAxis = createCoverageCoordAxis1D(AxisType.GeoX, xValues);
     final CoverageCoordAxis1D yAxis = createCoverageCoordAxis1D(AxisType.GeoY, yValues);
 
-    final AttributeContainerMutable attributes = new AttributeContainerMutable("attributes");
-    attributes.addAttribute(CF.GRID_MAPPING_NAME, CF.GEOSTATIONARY);
-    attributes.addAttribute(CF.LONGITUDE_OF_PROJECTION_ORIGIN, -75.0);
-    attributes.addAttribute(CF.PERSPECTIVE_POINT_HEIGHT, 35786023.0);
-    attributes.addAttribute(CF.SEMI_MINOR_AXIS, 6356752.31414);
-    attributes.addAttribute(CF.SEMI_MAJOR_AXIS, 6378137.0);
-    attributes.addAttribute(CF.INVERSE_FLATTENING, 298.2572221);
-    attributes.addAttribute(CF.SWEEP_ANGLE_AXIS, "x");
+    final AttributeContainerMutable attributes = createGeostationaryAttributes();
 
     final CoverageTransform transform = new CoverageTransform("transform", attributes, true);
     final HorizCoordSys horizCoordSys = HorizCoordSys.factory(xAxis, yAxis, null, null, transform);
@@ -62,6 +58,27 @@ public void shouldRemoveNansWhenComputingLatLon() {
     }
   }
 
+  @Test
+  public void shouldHandleUnitsInProjectionWhenSubsetting() {
+    final double[] xValues = new double[] {0, -24444.000000022017};
+    final double[] yValues = new double[] {0, 98028.00000000928};
+
+    final CoverageCoordAxis1D xAxis = createCoverageCoordAxis1D(AxisType.GeoX, "microradians", xValues);
+    final CoverageCoordAxis1D yAxis = createCoverageCoordAxis1D(AxisType.GeoY, "microradians", yValues);
+
+    final AttributeContainerMutable attributes = createGeostationaryAttributes();
+
+    final CoverageTransform transform = new CoverageTransform("transform", attributes, true);
+    final HorizCoordSys horizCoordSys = HorizCoordSys.factory(xAxis, yAxis, null, null, transform);
+
+    final LatLonPoint latLon = LatLonPoint.create(35, -85);
+    final SubsetParams subsetParams = new SubsetParams().setLatLonPoint(latLon);
+    final Optional<HorizCoordSys> subsetHorizCoordSys = horizCoordSys.subset(subsetParams);
+    assertThat(subsetHorizCoordSys.isPresent()).isTrue();
+    assertThat(subsetHorizCoordSys.get().getXAxis().getCoord(0)).isWithin(TOLERANCE).of(xValues[1]);
+    assertThat(subsetHorizCoordSys.get().getYAxis().getCoord(0)).isWithin(TOLERANCE).of(yValues[1]);
+  }
+
   @Test
   public void shouldUsePeriodsAsDecimalSeparatorsInWKT() throws IOException, URISyntaxException {
     Locale.setDefault(new Locale("fr", "FR"));
@@ -81,8 +98,24 @@ public void shouldUsePeriodsAsDecimalSeparatorsInWKT() throws IOException, URISy
     }
   }
 
+  private AttributeContainerMutable createGeostationaryAttributes() {
+    final AttributeContainerMutable attributes = new AttributeContainerMutable("attributes");
+    attributes.addAttribute(CF.GRID_MAPPING_NAME, CF.GEOSTATIONARY);
+    attributes.addAttribute(CF.LONGITUDE_OF_PROJECTION_ORIGIN, -75.0);
+    attributes.addAttribute(CF.PERSPECTIVE_POINT_HEIGHT, 35786023.0);
+    attributes.addAttribute(CF.SEMI_MINOR_AXIS, 6356752.31414);
+    attributes.addAttribute(CF.SEMI_MAJOR_AXIS, 6378137.0);
+    attributes.addAttribute(CF.INVERSE_FLATTENING, 298.2572221);
+    attributes.addAttribute(CF.SWEEP_ANGLE_AXIS, "x");
+    return attributes;
+  }
+
   private CoverageCoordAxis1D createCoverageCoordAxis1D(AxisType type, double[] values) {
-    final CoverageCoordAxisBuilder coordAxisBuilder = new CoverageCoordAxisBuilder("name", "unit", "description",
+    return createCoverageCoordAxis1D(type, "units", values);
+  }
+
+  private CoverageCoordAxis1D createCoverageCoordAxis1D(AxisType type, String units, double[] values) {
+    final CoverageCoordAxisBuilder coordAxisBuilder = new CoverageCoordAxisBuilder("name", units, "description",
         DataType.DOUBLE, type, null, CoverageCoordAxis.DependenceType.independent, null, Spacing.irregularPoint,
         values.length, values[0], values[values.length - 1], values[1] - values[0], values, null);
     return new CoverageCoordAxis1D(coordAxisBuilder);