All Downloads are FREE. Search and download functionalities are using the official Maven repository.

ucar.nc2.ft2.coverage.HorizCoordSys2D Maven / Gradle / Ivy

The newest version!
/*
 * Copyright (c) 1998-2018 John Caron and University Corporation for Atmospheric Research/Unidata
 * See LICENSE for license information.
 */
package ucar.nc2.ft2.coverage;

import ucar.ma2.*;
import ucar.nc2.util.Optional;
import ucar.unidata.geoloc.LatLonPoint;
import ucar.unidata.geoloc.LatLonPointImpl;
import ucar.unidata.geoloc.LatLonRect;

import java.util.ArrayList;
import java.util.Formatter;
import java.util.List;

/**
 * HorizCoordSys with LatLonAxis2D latAxis, lonAxis
 * some code forked from ucar.nc2.dataset.GridCoordinate2D
 */
public class HorizCoordSys2D extends HorizCoordSys {
  static private boolean debug = false;
  static private final org.slf4j.Logger log = org.slf4j.LoggerFactory.getLogger(HorizCoordSys2D.class);

  // private final LatLonAxis2D latAxis, lonAxis;
  private final int nrows, ncols;
  private Edges edges;

   HorizCoordSys2D(LatLonAxis2D latCoord, LatLonAxis2D lonCoord) {
    super(null, null, latCoord, lonCoord, null);
    int[] shape = latCoord.getShape(); // y, x

    nrows = shape[0];
    ncols = shape[1];
  }

  @Override
  public boolean isLatLon2D() {
    return true;
  }

    @Override
  public Optional subset(SubsetParams params) {

    LatLonRect llbb = (LatLonRect) params.get(SubsetParams.latlonBB);
    Integer horizStride = (Integer) params.get(SubsetParams.horizStride);
    if (horizStride == null || horizStride < 1) horizStride = 1;

    LatLonAxis2D lataxisSubset = null, lonaxisSubset = null;

    Formatter errMessages = new Formatter();
    if (llbb != null) {
      Optional> opt = computeBounds(llbb, horizStride);
      if (!opt.isPresent()) {
        errMessages.format("%s;%n", opt.getErrorMessage());
      } else {
        List ranges = opt.get();
        lataxisSubset = latAxis2D.subset(ranges.get(1), ranges.get(0));  // x, y
        lonaxisSubset = lonAxis2D.subset(ranges.get(1), ranges.get(0));
      }
    } else if (horizStride > 1) {
      try {
        lataxisSubset = latAxis2D.subset(new Range(0, ncols-1, horizStride), new Range(0, nrows-1, horizStride));
        lonaxisSubset = lonAxis2D.subset(new Range(0, ncols-1, horizStride), new Range(0, nrows-1, horizStride));
      } catch (InvalidRangeException e) {
        errMessages.format("%s;%n", e.getMessage());
      }
    }

    String errs = errMessages.toString();
    if (errs.length() > 0)
      return Optional.empty(errs);

    // makes a copy of the axis
    if (lataxisSubset == null) lataxisSubset = (LatLonAxis2D) latAxis2D.copy();
    if (lonaxisSubset == null) lonaxisSubset = (LatLonAxis2D) lonAxis2D.copy();

    return Optional.of(new HorizCoordSys2D(lataxisSubset, lonaxisSubset));
  }

  @Override
  public LatLonPoint getLatLon(int yindex, int xindex) {
    double lat = latAxis2D.getCoord(yindex, xindex);
    double lon = lonAxis2D.getCoord(yindex, xindex);
    return new LatLonPointImpl(lat, lon);
  }

  @Override
  public List getRanges() {
    return latAxis2D.getRanges(); // both are the same
  }

  @Override
  public Optional findXYindexFromCoord(double x, double y) {
    synchronized (this) {
      if (edges == null) edges = new Edges();
    }
    CoordReturn result = new CoordReturn();
    int[] index = new int[2];
    boolean ok = edges.findCoordElement(y, x, index);
    if (!ok)
      return Optional.empty("not in grid2D");

    result.x = index[1];
    result.y = index[0];
    result.xcoord = getLonAxis2D().getCoord(result.y, result.x);
    result.ycoord = getLatAxis2D().getCoord(result.y, result.x);
    return Optional.of(result);
  }

  @Override
  public List getCoordAxes() {
    List result = new ArrayList<>();
    result.add(latAxis2D);
    result.add(lonAxis2D);
    return result;
  }

  // return y, x ranges
  private Optional> computeBounds(LatLonRect llbb, int horizStride) {
    synchronized (this) {
      if (edges == null) edges = new Edges();
    }
    return edges.computeBoundsExhaustive(llbb, horizStride);
  }

  // assume this class is instantiated when these edges are needed
  private class Edges {
    private ArrayDouble.D2 latEdge, lonEdge;
    private MAMath.MinMax latMinMax, lonMinMax;

    Edges() {
      latEdge = (ArrayDouble.D2) latAxis2D.getCoordBoundsAsArray();
      lonEdge = (ArrayDouble.D2) lonAxis2D.getCoordBoundsAsArray();

      // assume missing values have been converted to NaNs
      latMinMax = MAMath.getMinMax(latEdge);
      lonMinMax = MAMath.getMinMax(lonEdge);

      // normalize to [minLon,minLon+360]
      int nlons = (int) lonEdge.getSize();
      for (int i = 0; i < nlons; i++) {
        double nonVal = lonEdge.getDouble(i);
        lonEdge.setDouble(i, LatLonPointImpl.lonNormalFrom(nonVal, lonMinMax.min));
      }

      if (debug)
        System.out.printf("Bounds (%d %d): lat= (%f,%f) lon = (%f,%f) %n",
                nrows, ncols, latMinMax.min, latMinMax.max, lonMinMax.min, lonMinMax.max);
    }

    /**
     * Find the best index for the given lat,lon point.
     *
     * @param wantLat   lat of point
     * @param wantLon   lon of point
     * @param rectIndex return (row,col) index, or best guess here. may not be null
     * @return false if not in the grid.
     */
    public boolean findCoordElement(double wantLat, double wantLon, int[] rectIndex) {
      double wantLonNormal = LatLonPointImpl.lonNormalFrom(wantLon, lonMinMax.min);
      return findCoordElementNoForce(wantLat, wantLonNormal, rectIndex);
    }

    private boolean findCoordElementNoForce(double wantLat, double wantLon, int[] rectIndex) {
      if (wantLat < latMinMax.min) return false;
      if (wantLat > latMinMax.max) return false;
      if (wantLon < lonMinMax.min) return false;
      if (wantLon > lonMinMax.max) return false;

      double gradientLat = (latMinMax.max - latMinMax.min) / nrows;
      double gradientLon = (lonMinMax.max - lonMinMax.min) / ncols;

      double diffLat = wantLat - latMinMax.min;
      double diffLon = wantLon - lonMinMax.min;

      // initial guess
      rectIndex[0] = (int) Math.round(diffLat / gradientLat);  // row
      rectIndex[1] = (int) Math.round(diffLon / gradientLon);  // col

      int count = 0;
      while (true) {
        count++;
        if (debug) System.out.printf("%nIteration %d %n", count);
        if (contains(wantLat, wantLon, rectIndex))
          return true;

        if (!jump2(wantLat, wantLon, rectIndex))
          return false;

        // bouncing around
        if (count > 10) {
          // last ditch attempt
          return incr(wantLat, wantLon, rectIndex);
          //if (!ok)
          //  log.error("findCoordElement didnt converge lat,lon = "+wantLat+" "+ wantLon);
          //return ok;
        }
      }
    }

    /**
     * Is the point (lat,lon) contained in the (row, col) rectangle ?
     *
     * @param wantLat   lat of point
     * @param wantLon   lon of point
     * @param rectIndex rectangle row, col, will be clipped to [0, nrows), [0, ncols)
     * @return true if contained
     */
    private boolean containsOld(double wantLat, double wantLon, int[] rectIndex) {
      rectIndex[0] = Math.max(Math.min(rectIndex[0], nrows - 1), 0);
      rectIndex[1] = Math.max(Math.min(rectIndex[1], ncols - 1), 0);

      int row = rectIndex[0];
      int col = rectIndex[1];

      if (debug) System.out.printf(" (%d,%d) contains (%f,%f) in (lat=%f %f) (lon=%f %f) ?%n",
              rectIndex[0], rectIndex[1], wantLat, wantLon,
              latEdge.get(row, col), latEdge.get(row + 1, col), lonEdge.get(row, col), lonEdge.get(row, col + 1));

      if (wantLat < latEdge.get(row, col)) return false;
      if (wantLat > latEdge.get(row + 1, col)) return false;
      if (wantLon < lonEdge.get(row, col)) return false;
      if (wantLon > lonEdge.get(row, col + 1)) return false;
      return true;
    }

    /**
     * Is the point (lat,lon) contained in the (row, col) rectangle ?
     *
     * @param wantLat   lat of point
     * @param wantLon   lon of point
     * @param rectIndex rectangle row, col, will be clipped to [0, nrows), [0, ncols)
     * @return true if contained
     */

/*
  http://mathforum.org/library/drmath/view/54386.html

      Given any three points on the plane (x0,y0), (x1,y1), and
  (x2,y2), the area of the triangle determined by them is
  given by the following formula:

        1 | x0 y0 1 |
    A = - | x1 y1 1 |,
        2 | x2 y2 1 |

  where the vertical bars represent the determinant.
  the value of the expression above is:

       (.5)(x1*y2 - y1*x2 -x0*y2 + y0*x2 + x0*y1 - y0*x1)

  The amazing thing is that A is positive if the three points are
  taken in a counter-clockwise orientation, and negative otherwise.

  To be inside a rectangle (or any convex body), as you trace
  around in a clockwise direction from p1 to p2 to p3 to p4 and
  back to p1, the "areas" of triangles p1 p2 p, p2 p3 p, p3 p4 p,
  and p4 p1 p must all be positive.  If you don't know the
  orientation of your rectangle, then they must either all be
  positive or all be negative.

  this method works for arbitrary convex regions on the plane.
*/
    private boolean contains(double wantLat, double wantLon, int[] rectIndex) {
      rectIndex[0] = Math.max(Math.min(rectIndex[0], nrows - 1), 0);
      rectIndex[1] = Math.max(Math.min(rectIndex[1], ncols - 1), 0);

      int row = rectIndex[0];
      int col = rectIndex[1];

      double x1 = lonEdge.get(row, col);
      double y1 = latEdge.get(row, col);

      double x2 = lonEdge.get(row, col + 1);
      double y2 = latEdge.get(row, col + 1);

      double x3 = lonEdge.get(row + 1, col + 1);
      double y3 = latEdge.get(row + 1, col + 1);

      double x4 = lonEdge.get(row + 1, col);
      double y4 = latEdge.get(row + 1, col);

      // must all have same determinate sign
      boolean sign = detIsPositive(x1, y1, x2, y2, wantLon, wantLat);
      if (sign != detIsPositive(x2, y2, x3, y3, wantLon, wantLat)) return false;
      if (sign != detIsPositive(x3, y3, x4, y4, wantLon, wantLat)) return false;
      if (sign != detIsPositive(x4, y4, x1, y1, wantLon, wantLat)) return false;

      return true;
    }

    private boolean detIsPositive(double x0, double y0, double x1, double y1, double x2, double y2) {
      double det = (x1 * y2 - y1 * x2 - x0 * y2 + y0 * x2 + x0 * y1 - y0 * x1);
      if (det == 0)
        log.warn("determinate = 0 on lat/lon=" + latAxis2D.getName() + ", " + lonAxis2D.getName()); // LOOK needed?
      return det > 0;
    }

    private boolean jumpOld(double wantLat, double wantLon, int[] rectIndex) {
      int row = Math.max(Math.min(rectIndex[0], nrows - 1), 0);
      int col = Math.max(Math.min(rectIndex[1], ncols - 1), 0);
      double gradientLat = latEdge.get(row + 1, col) - latEdge.get(row, col);
      double gradientLon = lonEdge.get(row, col + 1) - lonEdge.get(row, col);

      double diffLat = wantLat - latEdge.get(row, col);
      double diffLon = wantLon - lonEdge.get(row, col);

      int drow = (int) Math.round(diffLat / gradientLat);
      int dcol = (int) Math.round(diffLon / gradientLon);

      if (debug) System.out.printf("   jump from %d %d (grad=%f %f) (diff=%f %f) (delta=%d %d)",
              row, col, gradientLat, gradientLon,
              diffLat, diffLon, drow, dcol);

      if ((drow == 0) && (dcol == 0)) {
        if (debug) System.out.printf("%n   incr:");
        return incr(wantLat, wantLon, rectIndex);
      } else {
        rectIndex[0] = Math.max(Math.min(row + drow, nrows - 1), 0);
        rectIndex[1] = Math.max(Math.min(col + dcol, ncols - 1), 0);
        if (debug) System.out.printf(" to (%d %d)%n", rectIndex[0], rectIndex[1]);
        if ((row == rectIndex[0]) && (col == rectIndex[1])) return false; // nothing has changed
      }

      return true;
    }

    /**
     * jump to a new row,col
     *
     * @param wantLat   want this lat
     * @param wantLon   want this lon
     * @param rectIndex starting row, col and replaced by new guess
     * @return true if new guess, false if failure
     */
  /*
    choose x, y such that (matrix multiply) :

    (wantx) = (fxx fxy)  (x)
    (wanty)   (fyx fyy)  (y)

     where fxx = d(fx)/dx  ~= delta lon in lon direction
     where fxy = d(fx)/dy  ~= delta lon in lat direction
     where fyx = d(fy)/dx  ~= delta lat in lon direction
     where fyy = d(fy)/dy  ~= delta lat in lat direction

    2 linear equations in 2 unknowns, solve in usual way
   */
    private boolean jump2(double wantLat, double wantLon, int[] rectIndex) {
      int row = Math.max(Math.min(rectIndex[0], nrows - 1), 0);
      int col = Math.max(Math.min(rectIndex[1], ncols - 1), 0);
      double lat = latEdge.get(row, col);
      double lon = lonEdge.get(row, col);
      double diffLat = wantLat - lat;
      double diffLon = wantLon - lon;

      double dlatdy = latEdge.get(row + 1, col) - lat;
      double dlatdx = latEdge.get(row, col + 1) - lat;
      double dlondx = lonEdge.get(row, col + 1) - lon;
      double dlondy = lonEdge.get(row + 1, col) - lon;

      // solve for dlon

      double dx = (diffLon - dlondy * diffLat / dlatdy) / (dlondx - dlatdx * dlondy / dlatdy);
      // double dy =  (diffLat - dlatdx * diffLon / dlondx) / (dlatdy - dlatdx * dlondy / dlondx);
      double dy = (diffLat - dlatdx * dx) / dlatdy;

      if (debug)
        System.out.printf("   jump from %d %d (dlondx=%f dlondy=%f dlatdx=%f dlatdy=%f) (diffLat,Lon=%f %f) (deltalat,Lon=%f %f)",
                row, col, dlondx, dlondy, dlatdx, dlatdy,
                diffLat, diffLon, dy, dx);

      int drow = (int) Math.round(dy);
      int dcol = (int) Math.round(dx);

      if ((drow == 0) && (dcol == 0)) {
        if (debug) System.out.printf("%n   incr:");
        return incr(wantLat, wantLon, rectIndex);
      } else {
        rectIndex[0] = Math.max(Math.min(row + drow, nrows - 1), 0);
        rectIndex[1] = Math.max(Math.min(col + dcol, ncols - 1), 0);
        if (debug) System.out.printf(" to (%d %d)%n", rectIndex[0], rectIndex[1]);
        if ((row == rectIndex[0]) && (col == rectIndex[1])) return false; // nothing has changed
      }

      return true;
    }

    private boolean incr(double wantLat, double wantLon, int[] rectIndex) {
      int row = rectIndex[0];
      int col = rectIndex[1];
      double diffLat = wantLat - latEdge.get(row, col);
      double diffLon = wantLon - lonEdge.get(row, col);

      if (Math.abs(diffLat) > Math.abs(diffLon)) { // try lat first
        rectIndex[0] = row + ((diffLat > 0) ? 1 : -1);
        if (contains(wantLat, wantLon, rectIndex)) return true;
        rectIndex[1] = col + ((diffLon > 0) ? 1 : -1);
        if (contains(wantLat, wantLon, rectIndex)) return true;
      } else {
        rectIndex[1] = col + ((diffLon > 0) ? 1 : -1);
        if (contains(wantLat, wantLon, rectIndex)) return true;
        rectIndex[0] = row + ((diffLat > 0) ? 1 : -1);
        if (contains(wantLat, wantLon, rectIndex)) return true;
      }

      // back to original, do box search
      rectIndex[0] = row;
      rectIndex[1] = col;
      return box9(wantLat, wantLon, rectIndex);
    }

    // we think its got to be in one of the 9 boxes around rectIndex
    private boolean box9(double wantLat, double wantLon, int[] rectIndex) {
      int row = rectIndex[0];
      int minrow = Math.max(row - 1, 0);
      int maxrow = Math.min(row + 1, nrows);

      int col = rectIndex[1];
      int mincol = Math.max(col - 1, 0);
      int maxcol = Math.min(col + 1, ncols);

      if (debug) System.out.printf("%n   box9:");
      for (int i = minrow; i <= maxrow; i++)
        for (int j = mincol; j <= maxcol; j++) {
          rectIndex[0] = i;
          rectIndex[1] = j;
          if (contains(wantLat, wantLon, rectIndex)) return true;
        }

      return false;
    }

    // return y, x ranges
    Optional> computeBoundsExhaustive(LatLonRect rect, int horizStride) {
      LatLonPointImpl llpt = rect.getLowerLeftPoint();
      LatLonPointImpl urpt = rect.getUpperRightPoint();

      double miny = llpt.getLatitude();
      double maxy = urpt.getLatitude();

      // normalize to [minLon,minLon+360], edge already normalized to this
      double minx = LatLonPointImpl.lonNormalFrom(llpt.getLongitude(), lonMinMax.min);
      double maxx = LatLonPointImpl.lonNormalFrom(urpt.getLongitude(), lonMinMax.min);

      int shape[] = lonAxis2D.getShape();
      int ny = shape[0];
      int nx = shape[1];
      int minCol = Integer.MAX_VALUE, minRow = Integer.MAX_VALUE;
      int maxCol = -1, maxRow = -1;

      boolean allX = (minx > lonMinMax.max && maxx > lonMinMax.max && minx > maxx);
      boolean allY =  (miny <= latMinMax.min && maxy >= latMinMax.max);
      if (allX && allY) {
        // return full set
        return Optional.of(getRanges());
      }

      if (minx > lonMinMax.max && maxx > lonMinMax.max && minx < maxx) { // otherwise ignoring minx > maxx
        return Optional.empty("no intersection");
      } else if (minx > lonMinMax.max && maxx > lonMinMax.max && minx > maxx) {
        minCol = 0;   // all of x
        maxCol = nx;
        minx = lonMinMax.min;
      } else if (minx > lonMinMax.min && maxx > lonMinMax.max) {
        maxCol = nx;
      } else if (minx > lonMinMax.max && maxx < lonMinMax.max) {
        minCol = 0;
        minx = lonMinMax.min;
      }

      // probably not needed
      if (miny <= latMinMax.min) {
        minRow = 0;
      } else if (maxy >= latMinMax.max) {
        maxRow = ny;
      }

      // brute force, examine every point LOOK BAD
      for (int row = 0; row <= ny; row++) {
        for (int col = 0; col <= nx; col++) {
          double lat = latEdge.get(row, col);
          double lon = lonEdge.get(row, col);

          if ((lat >= miny) && (lat <= maxy) && (lon >= minx) && (lon <= maxx)) {
            if (col > maxCol) maxCol = col;
            if (col < minCol) minCol = col;
            if (row > maxRow) maxRow = row;
            if (row < minRow) minRow = row;
          }
        }
      }

      try {
        List list = new ArrayList<>();
        list.add(new Range(minRow, maxRow-1, horizStride));
        list.add(new Range(minCol, maxCol-1, horizStride));
        return Optional.of(list);

      } catch (InvalidRangeException e) {
        throw new RuntimeException(e);
      }

    }

    private double getMinOrMaxLon(double lon1, double lon2, boolean wantMin) {
      double midpoint = (lon1 + lon2) / 2;
      lon1 = LatLonPointImpl.lonNormal(lon1, midpoint);
      lon2 = LatLonPointImpl.lonNormal(lon2, midpoint);

      return wantMin ? Math.min(lon1, lon2) : Math.max(lon1, lon2);
    }

    // brute force
    public boolean findCoordElementExhaustive(double wantLat, double wantLon, int[] rectIndex) {
      if (wantLat < latMinMax.min) return false;
      if (wantLat > latMinMax.max) return false;
      if (wantLon < lonMinMax.min) return false;
      if (wantLon > lonMinMax.max) return false;

      for (int row = 0; row < nrows; row++) {
        for (int col = 0; col < ncols; col++) {

          rectIndex[0] = row;
          rectIndex[1] = col;

          if (contains(wantLat, wantLon, rectIndex)) {
            return true;
          }
        }
      }
      //debug = saveDebug;
      return false;
    }
  }


  /*
  public int[] findXYindexFromCoord(double x_coord, double y_coord) {
      int[] result = new int[2];

      int[] result2 = new int[2];
      boolean found = g2d.findCoordElement(y_coord, x_coord, result2);
      if (found) {
        result[0] = result2[1];
        result[1] = result2[0];
      } else {
        result[0] = -1;
        result[1] = -1;
      }
      return result;
  } */

}





© 2015 - 2024 Weber Informatics LLC | Privacy Policy