Skip to content

Commit

Permalink
Merge branch 'maint-5.x' into jmasson-netcdf-java-1235
Browse files Browse the repository at this point in the history
  • Loading branch information
@mnlerman
mnlerman authored Sep 19, 2023
2 parents 6282b64 + 8935ad7 commit 0f8ad82
Show file tree
Hide file tree
Showing 6 changed files with 119 additions and 70 deletions.
22 changes: 1 addition & 21 deletions cdm/core/src/main/java/ucar/nc2/dataset/transform/Geostationary.java
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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);

Expand Down Expand Up @@ -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);
}
Expand Down
61 changes: 28 additions & 33 deletions cdm/core/src/main/java/ucar/nc2/ft2/coverage/HorizCoordSys.java
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand Down Expand Up @@ -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) {
Expand All @@ -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;
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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();
Expand Down Expand Up @@ -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);
Expand Down Expand Up @@ -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())) {
Expand All @@ -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) {
Expand Down
20 changes: 20 additions & 0 deletions cdm/core/src/main/java/ucar/unidata/geoloc/ProjectionImpl.java
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
*/
Expand All @@ -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;
Expand Down Expand Up @@ -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
*
Expand Down
5 changes: 4 additions & 1 deletion cdm/core/src/main/java/ucar/unidata/geoloc/projection/LambertConformal.java
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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;
Expand Down
30 changes: 24 additions & 6 deletions cdm/core/src/main/java/ucar/unidata/geoloc/projection/sat/Geostationary.java
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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;

Expand All @@ -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) {
Expand All @@ -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) {
Expand All @@ -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);

Expand All @@ -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;
}
Expand Down
51 changes: 42 additions & 9 deletions cdm/core/src/test/java/ucar/nc2/ft2/coverage/TestHorizCoordSys.java
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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() {
Expand All @@ -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);
Expand All @@ -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"));
Expand All @@ -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);
Expand Down

0 comments on commit 0f8ad82

Please sign in to comment.