ucar.nc2.dataset.conv.ADASConvention Maven / Gradle / Ivy

Go to download
Show more of this group Show more artifacts with this name
Show all versions of cdm Show documentation
The NetCDF-Java Library is a Java interface to NetCDF files, as well as to many other types of scientific data formats.
The newest version!
/*
 * Copyright 1998-2014 University Corporation for Atmospheric Research/Unidata
 *
 *   Portions of this software were developed by the Unidata Program at the
 *   University Corporation for Atmospheric Research.
 *
 *   Access and use of this software shall impose the following obligations
 *   and understandings on the user. The user is granted the right, without
 *   any fee or cost, to use, copy, modify, alter, enhance and distribute
 *   this software, and any derivative works thereof, and its supporting
 *   documentation for any purpose whatsoever, provided that this entire
 *   notice appears in all copies of the software, derivative works and
 *   supporting documentation.  Further, UCAR requests that the user credit
 *   UCAR/Unidata in any publications that result from the use of this
 *   software or in any product that includes this software. The names UCAR
 *   and/or Unidata, however, may not be used in any advertising or publicity
 *   to endorse or promote any products or commercial entity unless specific
 *   written permission is obtained from UCAR/Unidata. The user also
 *   understands that UCAR/Unidata is not obligated to provide the user with
 *   any support, consulting, training or assistance of any kind with regard
 *   to the use, operation and performance of this software nor to provide
 *   the user with any updates, revisions, new versions or "bug fixes."
 *
 *   THIS SOFTWARE IS PROVIDED BY UCAR/UNIDATA "AS IS" AND ANY EXPRESS OR
 *   IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
 *   WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
 *   DISCLAIMED. IN NO EVENT SHALL UCAR/UNIDATA BE LIABLE FOR ANY SPECIAL,
 *   INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES WHATSOEVER RESULTING
 *   FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN ACTION OF CONTRACT,
 *   NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF OR IN CONNECTION
 *   WITH THE ACCESS, USE OR PERFORMANCE OF THIS SOFTWARE.
 */

package ucar.nc2.dataset.conv;

import ucar.ma2.*;
import ucar.nc2.*;
import ucar.nc2.constants.CDM;
import ucar.nc2.constants.CF;
import ucar.nc2.constants._Coordinate;
import ucar.nc2.constants.AxisType;
import ucar.nc2.units.SimpleUnit;
import ucar.nc2.util.CancelTask;
import ucar.nc2.dataset.*;
import ucar.nc2.dataset.transform.WRFEtaTransformBuilder;

import ucar.unidata.geoloc.*;
import ucar.unidata.geoloc.projection.*;
//import ucar.units.ConversionException;

import java.io.IOException;
import java.util.*;

/**
 * ADAS netcdf files.
 *
 * @author caron
 */

public class ADASConvention extends CoordSysBuilder {

  public ADASConvention() {
    this.conventionName = "ARPS/ADAS";
  }

  // private double originX = 0.0, originY = 0.0;
  private ProjectionCT projCT = null;
  private boolean debugProj = false;

  /**
   * create a NetcdfDataset out of this NetcdfFile, adding coordinates etc.
   */
  public void augmentDataset(NetcdfDataset ds, CancelTask cancelTask) throws IOException {
    if (null != ds.findVariable("x")) return; // check if its already been done - aggregating enhanced datasets.

    // old way
    Attribute att = ds.findGlobalAttribute("MAPPROJ");
    int projType = att.getNumericValue().intValue();

    double lat1 = findAttributeDouble(ds, "TRUELAT1", Double.NaN);
    double lat2 = findAttributeDouble(ds, "TRUELAT2", Double.NaN);
    double lat_origin = lat1;
    double lon_origin = findAttributeDouble(ds, "TRUELON", Double.NaN);
    double false_easting = 0.0;
    double false_northing = 0.0;

    // new way
    String projName = ds.findAttValueIgnoreCase(null, CF.GRID_MAPPING_NAME, null);
    if (projName != null) {
      projName = projName.trim();
      lat_origin = findAttributeDouble(ds, "latitude_of_projection_origin", Double.NaN);
      lon_origin = findAttributeDouble(ds, "longitude_of_central_meridian", Double.NaN);
      false_easting = findAttributeDouble(ds, "false_easting", 0.0);
      false_northing = findAttributeDouble(ds, "false_northing", 0.0);

      Attribute att2 = ds.findGlobalAttributeIgnoreCase("standard_parallel");
      if (att2 != null) {
        lat1 = att2.getNumericValue().doubleValue();
        lat2 = (att2.getLength() > 1) ? att2.getNumericValue(1).doubleValue() : lat1;
      }
    } else {
      if (projType == 2) projName = "lambert_conformal_conic";
    }

    Variable coord_var = ds.findVariable("x_stag");
    if (!Double.isNaN(false_easting) || !Double.isNaN(false_northing)) {
      String units = ds.findAttValueIgnoreCase(coord_var, CDM.UNITS, null);
      double scalef = 1.0;
      try {
        scalef = SimpleUnit.getConversionFactor(units, "km");
      } catch (IllegalArgumentException e) {
        log.error(units + " not convertible to km");
      }
      false_easting *= scalef;
      false_northing *= scalef;
    }

    ProjectionImpl proj = null;
    if ((projName != null) && projName.equalsIgnoreCase("lambert_conformal_conic")) {
      proj = new LambertConformal(lat_origin, lon_origin, lat1, lat2, false_easting, false_northing);
      projCT = new ProjectionCT("Projection", "FGDC", proj);
      if (false_easting == 0.0) calcCenterPoints(ds, proj); // old way
    } else {
      parseInfo.format("ERROR: unknown projection type = %s%n", projName);
    }

    if (debugProj && (proj != null)) {
      System.out.println(" using LC " + proj.paramsToString());
      double lat_check = findAttributeDouble(ds, "CTRLAT", Double.NaN);
      double lon_check = findAttributeDouble(ds, "CTRLON", Double.NaN);

      LatLonPointImpl lpt0 = new LatLonPointImpl(lat_check, lon_check);
      ProjectionPoint ppt0 = proj.latLonToProj(lpt0, new ProjectionPointImpl());
      System.out.println("CTR lpt0= " + lpt0 + " ppt0=" + ppt0);

      Variable xstag = ds.findVariable("x_stag");
      ArrayFloat.D1 xstagData = (ArrayFloat.D1) xstag.read();
      float center_x = xstagData.get((int) xstag.getSize() - 1);
      Variable ystag = ds.findVariable("y_stag");
      ArrayFloat.D1 ystagData = (ArrayFloat.D1) ystag.read();
      float center_y = ystagData.get((int) ystag.getSize() - 1);
      System.out.println("CTR should be x,y= " + center_x / 2000 + ", " + center_y / 2000);

      lpt0 = new LatLonPointImpl(lat_origin, lon_origin);
      ppt0 = proj.latLonToProj(lpt0, new ProjectionPointImpl());
      System.out.println("ORIGIN lpt0= " + lpt0 + " ppt0=" + ppt0);

      lpt0 = new LatLonPointImpl(lat_origin, lon_origin);
      ppt0 = proj.latLonToProj(lpt0, new ProjectionPointImpl());
      System.out.println("TRUE ORIGIN lpt0= " + lpt0 + " ppt0=" + ppt0);
    }

    if (projCT != null) {
      VariableDS v = makeCoordinateTransformVariable(ds, projCT);
      v.addAttribute(new Attribute(_Coordinate.AxisTypes, "GeoX GeoY"));
      ds.addVariable(null, v);
    }

    if (ds.findVariable("x_stag") != null)
      ds.addCoordinateAxis(makeCoordAxis(ds, "x"));
    if (ds.findVariable("y_stag") != null)
      ds.addCoordinateAxis(makeCoordAxis(ds, "y"));
    if (ds.findVariable("z_stag") != null)
      ds.addCoordinateAxis(makeCoordAxis(ds, "z"));

    Variable zsoil = ds.findVariable("ZPSOIL");
    if (zsoil != null)
      zsoil.addAttribute(new Attribute(_Coordinate.AxisType, AxisType.GeoZ.toString()));

    ds.finish();
  }

  // old
  private void calcCenterPoints(NetcdfDataset ds, Projection proj) throws IOException {
    double lat_check = findAttributeDouble(ds, "CTRLAT", Double.NaN);
    double lon_check = findAttributeDouble(ds, "CTRLON", Double.NaN);

    LatLonPointImpl lpt0 = new LatLonPointImpl(lat_check, lon_check);
    ProjectionPoint ppt0 = proj.latLonToProj(lpt0, new ProjectionPointImpl());
    System.out.println("CTR lpt0= " + lpt0 + " ppt0=" + ppt0);

    Variable xstag = ds.findVariable("x_stag");
    int nxpts = (int) xstag.getSize();
    ArrayFloat.D1 xstagData = (ArrayFloat.D1) xstag.read();
    float center_x = xstagData.get(nxpts - 1);
    double false_easting = center_x / 2000 - ppt0.getX() * 1000.0;
    System.out.println("false_easting= " + false_easting);

    Variable ystag = ds.findVariable("y_stag");
    int nypts = (int) ystag.getSize();
    ArrayFloat.D1 ystagData = (ArrayFloat.D1) ystag.read();
    float center_y = ystagData.get(nypts - 1);
    double false_northing = center_y / 2000 - ppt0.getY() * 1000.0;
    System.out.println("false_northing= " + false_northing);

    double dx = findAttributeDouble(ds, "DX", Double.NaN);
    double dy = findAttributeDouble(ds, "DY", Double.NaN);

    double w = dx * (nxpts - 1);
    double h = dy * (nypts - 1);
    double startx = ppt0.getX() * 1000.0 - w / 2;
    double starty = ppt0.getY() * 1000.0 - h / 2;

    xstag.setValues(nxpts, startx, dx);
    ystag.setValues(nypts, starty, dy);
  }

  /////////////////////////////////////////////////////////////////////////

  protected void makeCoordinateTransforms(NetcdfDataset ds) {
    if (projCT != null) {
      VarProcess vp = findVarProcess(projCT.getName(), null);
      vp.isCoordinateTransform = true;
      vp.ct = projCT;
    }
    super.makeCoordinateTransforms(ds);
  }

  protected AxisType getAxisType(NetcdfDataset ds, VariableEnhanced ve) {
    Variable v = (Variable) ve;
    String vname = v.getShortName();

    if (vname.equalsIgnoreCase("x") || vname.equalsIgnoreCase("x_stag"))
      return AxisType.GeoX;

    if (vname.equalsIgnoreCase("lon"))
      return AxisType.Lon;

    if (vname.equalsIgnoreCase("y") || vname.equalsIgnoreCase("y_stag"))
      return AxisType.GeoY;

    if (vname.equalsIgnoreCase("lat"))
      return AxisType.Lat;

    if (vname.equalsIgnoreCase("z") || vname.equalsIgnoreCase("z_stag"))
      return AxisType.GeoZ;

    if (vname.equalsIgnoreCase("Z"))
      return AxisType.Height;

    if (vname.equalsIgnoreCase("time"))
      return AxisType.Time;

    String unit = ve.getUnitsString();
    if (unit != null) {
      if (SimpleUnit.isCompatible("millibar", unit))
        return AxisType.Pressure;

      if (SimpleUnit.isCompatible("m", unit))
        return AxisType.Height;
    }


    return null;
  }

  /**
   * Does increasing values of Z go vertical  up?
   *
   * @param v for this axis
   * @return "up" if this is a Vertical (z) coordinate axis which goes up as coords get bigger,
   *         else return "down"
   */
  public String getZisPositive(CoordinateAxis v) {
    return "down"; //eta coords decrease upward
  }

  //////////////////////////////////////////////////////////////////////////////////////////////

  private CoordinateAxis makeCoordAxis(NetcdfDataset ds, String axisName) throws IOException {
    Variable stagV = ds.findVariable(axisName + "_stag");
    Array data_stag = stagV.read();
    int n = (int) data_stag.getSize() - 1;
    Array data = Array.factory(data_stag.getElementType(), new int[]{n});
    Index stagIndex = data_stag.getIndex();
    Index dataIndex = data.getIndex();
    for (int i = 0; i < n; i++) {
      double val = data_stag.getDouble(stagIndex.set(i)) + data_stag.getDouble(stagIndex.set(i + 1));
      data.setDouble(dataIndex.set(i), 0.5 * val);
    }

    DataType dtype = DataType.getType(data.getElementType());
    String units = ds.findAttValueIgnoreCase(stagV, CDM.UNITS, "m");
    CoordinateAxis v = new CoordinateAxis1D(ds, null, axisName, dtype, axisName, units, "synthesized non-staggered " + axisName + " coordinate");
    v.setCachedData(data, true);
    return v;
  }

  private double findAttributeDouble(NetcdfDataset ds, String attname, double defValue) {
    Attribute att = ds.findGlobalAttributeIgnoreCase(attname);
    if (att == null) return defValue;
    return att.getNumericValue().doubleValue();
  }

  /* private VerticalCT makeWRFEtaVerticalCoordinateTransform(NetcdfDataset ds, CoordinateSystem cs) {
    if ((null == ds.findVariable("PH")) || (null == ds.findVariable("PHB")) ||
        (null == ds.findVariable("P")) || (null == ds.findVariable("PB")))
      return null;

  	VerticalCT.Type type = VerticalCT.Type.WRFEta;
  	VerticalCT ct = new VerticalCT(type.toString(), conventionName, type);

  	ct.addParameter(new Parameter("height formula", "height(x,y,z) = (PH(x,y,z) + PHB(x,y,z)) / 9.81"));
  	ct.addParameter(new Parameter("perturbation geopotential variable name", "PH"));
  	ct.addParameter(new Parameter("base state geopotential variable name", "PHB"));
  	ct.addParameter(new Parameter("pressure formula", "pressure(x,y,z) = P(x,y,z) + PB(x,y,z)"));
  	ct.addParameter(new Parameter("perturbation pressure variable name", "P"));
  	ct.addParameter(new Parameter("base state pressure variable name", "PB"));
  	ct.addParameter(new Parameter("staggered x", ""+isStaggered(cs.getXaxis())));
  	ct.addParameter(new Parameter("staggered y", ""+isStaggered(cs.getYaxis())));
  	ct.addParameter(new Parameter("staggered z", ""+isStaggered(cs.getZaxis())));
  	ct.addParameter(new Parameter("eta", ""+cs.getZaxis().getName()));

    parseInfo.append(" added vertical coordinate transform = "+type+"\n");
  	return ct;
  }

  private boolean isStaggered(CoordinateAxis axis) {
  	if (axis == null) return false;
  	String name = axis.getName();
  	if (name == null) return false;
  	if (name.endsWith("stag")) return true;
  	return false;
  } */

  /**
   * Assign CoordinateTransform objects to Coordinate Systems.
   */
  protected void assignCoordinateTransforms(NetcdfDataset ncDataset) {
    super.assignCoordinateTransforms(ncDataset);

    // any cs whose got a vertical coordinate with no units
    List csys = ncDataset.getCoordinateSystems();
    for (CoordinateSystem cs : csys) {
      if (cs.getZaxis() != null) {
        String units = cs.getZaxis().getUnitsString();
        if ((units == null) || (units.trim().length() == 0)) {
          VerticalCT vct = makeWRFEtaVerticalCoordinateTransform(ncDataset, cs);
          if (vct != null) {
            cs.addCoordinateTransform(vct);
            parseInfo.format("***Added WRFEta verticalCoordinateTransform to %s%n", cs.getName());
          }
        }
      }
    }
  }

  private VerticalCT makeWRFEtaVerticalCoordinateTransform(NetcdfDataset ds, CoordinateSystem cs) {
    if ((null == ds.findVariable("PH")) || (null == ds.findVariable("PHB")) ||
        (null == ds.findVariable("P")) || (null == ds.findVariable("PB")))
      return null;

    WRFEtaTransformBuilder builder = new WRFEtaTransformBuilder(cs);
    return (VerticalCT) builder.makeCoordinateTransform(ds, null);
  }

}