Writing a Coordinate Transform: Projections and Vertical Transforms

Overview

A CoordinateTransform represents a mathematical function that transforms a dataset’s coordinates into coordinates in a reference CoordinateSystem. Currently, the CDM has two kinds of transforms: Projections and VerticalTransforms. A Projection maps between cartesian x and y coordinates (called GeoX and GeoY) and latitude, longitude coordinates, by implementing the ucar.unidata.geoloc.Projection interface. A VerticalTransform takes a GeoZ coordinate and usually other data fields such as surface pressure, and produces a 3D height or pressure vertical coordinate field.

A CoordinateSystem may have 0 or more CoordinateTransforms, each of which is either a ProjectionCT containing a ucar.unidata.geoloc.Projection or a VerticalCT containing a ucar.unidata.geoloc.vertical.VerticalTransform:

Coord Sys Object Model

The netCDF-Java library implements a standard set of ucar.unidata.geoloc.Projection and ucar.unidata.geoloc.vertical.VerticalTransform classes, following the specifications of the CF-1.0 Conventions.

Implementing a Coordinate Transform

The steps to using your own CoordinateTransform in the netCDF-Java library:

Write a class that implements ucar.nc2.dataset.CoordTransBuilderIF, by subclassing ucar.nc2.dataset.transform.AbstractTransformBuilder. Add these classes to your classpath. From your application, call ucar.nc2.dataset.CoordTransBuilder.registerTransform( String transformName, Class c). The Coordinate System Builder for your dataset must recognize the transform and add it to the coordinate system. If you use the CF-1.0 or the _Coordinate conventions, this means that the dataset must contain a CoordinateTransform variable that contains the parameters of the transform.

The classes that you will use are shown in the following diagram, which has an example of both a Projection ( LambertConformal) and a VerticalTransform (CFOceanSigma).

Coordinate Transforms

Implementing a Projection

You should implement the ucar.unidata.geoloc.Projection interface by subclassing the abstract class ucar.unidata.geoloc.projection.ProjectionImpl. The methods you need to implement are:

  public ProjectionPoint latLonToProj(LatLonPoint latlon, ProjectionPointImpl destPoint);
  public LatLonPoint projToLatLon(ProjectionPoint ppt, LatLonPointImpl destPoint);
  public boolean crossSeam(ProjectionPoint pt1, ProjectionPoint pt2);
  public ProjectionImpl constructCopy();
  public boolean equals(Object proj);
  public int hashCode()
  public String paramsToString();

The latLonToProj() and projToLatLon() methods are inverses of each other, and map between lat, lon (in units of decimal degrees) to cartesian x,y, the coordinates that your dataset uses, usually in units of “km on the projection plane”. The crossSeam() method returns true when a line between two points in projection coordinates would cross a seam in the projection plane, such as for a cylindrical or conic projections. This helps drawing routines to eliminate spurious lines. The constructCopy() method constructs a new, equivalent Projection object, which avoids the problems with clone (see Bloch, Effective Java, item 10). The equals() method should be overridden to make Projections equal that have the same parameters. You should also override hashCode() to make it consistent with equals (see Bloch, Effective Java, item 8). The paramsToString() returns a String representation of the Projection parameters. Examine the classes in ucar.unidata.geoloc.projection for implementation examples.

Implementing a VerticalTransform

You should implement the ucar.unidata.geoloc.vertical.VerticalTransform interface by subclassing the abstract class ucar.unidata.geoloc.vertical.VerticalTransformImpl. The methods you need to implement are:

  public ucar.ma2.ArrayDouble.D3 getCoordinateArray(int timeIndex);
  public String getUnitString();

The getCoordinateArray() method returns a 3D vertical coordinate for the given time step (the time step is ignored if isTimeDependent() is false). The returned array must use dimensions in the order of z, y, x. The getUnitString() method returns the unit of the transformed vertical coordinate, which should be udunits compatible with height or pressure. Examine the classes in ucar.unidata.geoloc.vertical for implementation examples.

Implementing and registering CoordTransBuilderIF

The Projection and VerticalTransform implement the mathematical transformation itself. Now we need to add the glue classes that allow runtime discovery and object instantiation. To do so, you must add a class that implements the ucar.nc2.dataset.CoordTransBuilderIF interface. You should subclass the abstract class ucar.nc2.dataset.transform.AbstractTransformBuilder, and implement the following methods:

  public String getTransformName();
  public TransformType getTransformType();
  public CoordinateTransform makeCoordinateTransform (NetcdfDataset ds, Variable ctv);
  public ucar.unidata.geoloc.vertical.VerticalTransform makeMathTransform(NetcdfDataset ds, Dimension timeDim, VerticalCT vCT);

Give your transform a unique name, which is returned by the getTransformName() method. The getTransformType() method should return either ucar.nc2.dataset.TransformType.Projection or TransformType.Vertical. The makeCoordinateTransform() method is the guts of the class, it takes as parameters the NetcdfDataset and the CoordinateTransform variable that contains the transformation parameters. The makeMathTransform() is used only for VerticalTransforms to defer the creation of the VerticalTransform until the CoordinateSystem has been fully constructed and, for example, the time dimension has been identified.

You then need to tell the netCDF-Java library about your transform class :

CoordTransBuilder.registerTransform("YourTransformName", YourTransformClass.class);

The name is the same as getTransformType() returns, and must be referenced in your dataset by the CoordinateTransform variable.

Projection Example

Following is an example from the standard implementation classes in ucar.nc2.internal.dataset.transform.

class LambertConformalConic extends AbstractTransformBuilder
    implements HorizTransformBuilderIF {

  public ProjectionCT makeCoordinateTransform(AttributeContainer ctv,
      String geoCoordinateUnits) {
    // IMPLEMENTED BELOW
  }

  public String getTransformName() {
    // The name of the transformation. This is referenced in your dataset
    return "lambert_conformal_conic";
  }
}

Implementation of makeCoordinateTransform()

Below is an example implementation of the makeCoordinateTransform() method.

double[] pars = readAttributeDouble2(ctv.findAttribute(CF.STANDARD_PARALLEL));
if (pars == null)
  return null;

double lon0 = readAttributeDouble(ctv, CF.LONGITUDE_OF_CENTRAL_MERIDIAN, Double.NaN);
double lat0 = readAttributeDouble(ctv, CF.LATITUDE_OF_PROJECTION_ORIGIN, Double.NaN);
double false_easting = readAttributeDouble(ctv, CF.FALSE_EASTING, 0.0);
double false_northing = readAttributeDouble(ctv, CF.FALSE_NORTHING, 0.0);

if ((false_easting != 0.0) || (false_northing != 0.0)) {
  double scalef = getFalseEastingScaleFactor(geoCoordinateUnits);
  false_easting *= scalef;
  false_northing *= scalef;
}

double earth_radius = getEarthRadiusInKm(ctv);
double semi_major_axis = readAttributeDouble(ctv, CF.SEMI_MAJOR_AXIS, Double.NaN);
double semi_minor_axis = readAttributeDouble(ctv, CF.SEMI_MINOR_AXIS, Double.NaN);
double inverse_flattening = readAttributeDouble(ctv, CF.INVERSE_FLATTENING, 0.0);

ucar.unidata.geoloc.ProjectionImpl proj;

// check for ellipsoidal earth
if (!Double.isNaN(semi_major_axis) && (!Double.isNaN(semi_minor_axis) || inverse_flattening != 0.0)) {
  Earth earth = new Earth(semi_major_axis, semi_minor_axis, inverse_flattening);
  proj = new ucar.unidata.geoloc.projection.proj4.LambertConformalConicEllipse(lat0, lon0, pars[0], pars[1],
      false_easting, false_northing, earth);

} else {
  proj = new ucar.unidata.geoloc.projection.LambertConformal(lat0, lon0, pars[0], pars[1], false_easting,
      false_northing, earth_radius);
}

return new ProjectionCT(ctv.getName(), "FGDC", proj);

Explanation of makeCoordinateTransform()

  • Various parameters are read from the attributes of the Coordinate Transform Variable.
  • A Projection is created from the parameters.
  • A ProjectionCT wraps the Projection and is returned.

Vertical Transform Example

Following is an example from the standard implementation classes in ucar.nc2.internal.dataset.transform.vertical.

public class CFOceanSigma extends AbstractTransformBuilder {
  private String sigma, ps, ptop;

  protected String getFormula(NetcdfDataset ds, Variable ctv) {
    String formula = ctv.findAttributeString("formula_terms", null);
    if (null == formula) {
      if (null != errBuffer)
        errBuffer.format(
            "CoordTransBuilder %s: needs attribute 'formula_terms' on Variable %s%n",
            getTransformName(), ctv.attributes().getName());
      return null;
    }
    return formula;
  }

  public VerticalCT makeCoordinateTransform(NetcdfDataset ds, AttributeContainer ctv) {
    // IMPLEMENTED BELOW
  }

  public String getTransformName() {
    // The name of the transformation, this is referenced in your dataset
    return "atmosphere_sigma_coordinate";
  }

  public OceanSigma makeMathTransform(NetcdfDataset ds, Dimension timeDim, VerticalCT vCT) {
    // VerticalTransform is created when method called by the VerticalCT
    return new OceanSigma(ds, timeDim, vCT.getParameters());
  }
}

Implementation of makeCoordinateTransform()

Below is an example implementation of the makeCoordinateTransform() method.

String formula_terms = getFormula(ctv);
if (null == formula_terms)
  return null;

String[] values = parseFormula(formula_terms, "sigma eta depth");
if (values == null)
  return null;

sigma = values[0];
eta = values[1];
depth = values[2];

VerticalCT rs =
    new VerticalCT("OceanSigma_Transform_" + ctv.getName(), getTransformName(), VerticalCT.Type.OceanSigma, this);
rs.addParameter(new Parameter("standard_name", getTransformName()));
rs.addParameter(new Parameter("formula_terms", formula_terms));
rs.addParameter((new Parameter("height_formula", "height(x,y,z) = eta(n,j,i) + sigma(k)*(depth(j,i)+eta(n,j,i))")));

if (!addParameter(rs, OceanSigma.SIGMA, ds, sigma))
  return null;
if (!addParameter(rs, OceanSigma.ETA, ds, eta))
  return null;
if (!addParameter(rs, OceanSigma.DEPTH, ds, depth))
  return null;

return rs;

Corresponding Vertical Transform Example

Following is an example from the standard implementation classes in ucar.unidata.geoloc.vertical.

class AtmosSigma extends VerticalTransformImpl {
  // The parameter names as public constant Strings
  public static final String PTOP = "Pressure at top";
  public static final String PS = "surfacePressure variable name";
  public static final String SIGMA = "sigma variable name";
  private Variable psVar; // surface pressure
  private double[] sigma; // The sigma array, function of z
  private double ptop; // Top of the model

  public AtmosSigma(Dimension timeDim) {
    // IMPLEMENTED BELOW
  }

  public ArrayDouble.D3 getCoordinateArray(int timeIndex) {
    // IMPLEMENTED BELOW
  }

  public ArrayDouble.D1 getCoordinateArray1D(int timeIndex, int xIndex, int yIndex) {
    // TO BE IMPLEMENTED
  }
}

Implementation of the constructor

Below is an example implementation of the AtmosSigma() constructor.

super(timeDim);

String psName = getParameterStringValue(params, PS);
psVar = ds.findVariable(psName);

String ptopName = getParameterStringValue(params, PTOP);
Variable ptopVar = ds.findVariable(ptopName);
try {
  this.ptop = ptopVar.readScalarDouble();
} catch (IOException e) {
  throw new IllegalArgumentException("AtmosSigma failed to read " + ptopVar + " err= " + e.getMessage());
}

String sigmaName = getParameterStringValue(params, SIGMA);
Variable sigmaVar = ds.findVariable(sigmaName);

try {
  Array data = sigmaVar.read();
  sigma = (double[]) data.get1DJavaArray(double.class);
} catch (IOException e) {
  throw new IllegalArgumentException("AtmosSigma failed to read " + sigmaName + " err= " + e.getMessage());
}

units = psVar.findAttributeString(CDM.UNITS, "none");

String ptopUnitStr = ptopVar.findAttributeString(CDM.UNITS, "none");
if (!units.equalsIgnoreCase(ptopUnitStr)) {
  // Convert ptopVar to units of psVar
  SimpleUnit psUnit = SimpleUnit.factory(units);
  SimpleUnit ptopUnit = SimpleUnit.factory(ptopUnitStr);
  double factor = ptopUnit.convertTo(1.0, psUnit);
  this.ptop = this.ptop * factor;
}

Explanation of AtmosSigma class constructor:

  • The psVar variable holding the surface pressure
  • The ptopVar variable is the value of the ptop scalar variable
  • The sigmaVar variable is the value of the sigma[z] coordinate
  • The returned converted coordinates will be in the units of the surface pressure

    Implementation of getCoordinateArray()

Below is an example implementation of the getCoordinateArray() method.

Array ps = readArray(psVar, timeIndex);
Index psIndex = ps.getIndex();

int nz = sigma.length;
int[] shape2D = ps.getShape();
int ny = shape2D[0];
int nx = shape2D[1];

ArrayDouble.D3 result = new ArrayDouble.D3(nz, ny, nx);

for (int y = 0; y < ny; y++) {
  for (int x = 0; x < nx; x++) {
    double psVal = ps.getDouble(psIndex.set(y, x));
    for (int z = 0; z < nz; z++) {
      result.set(z, y, x, ptop + sigma[z] * (psVal - ptop));
    }
  }
}

return result;

Explanation of getCoordinateArray():

  • Reads the surface pressure variable at the given time step through a utility method in the superclass
  • Creates the result array
  • Extracts the surface pressure at the given x,y point
  • Loops over z, the converted coordinate = ptop + sigma(z)*(surfacePressure(x,y)-ptop)