ucar.nc2.ft2.coverage.adapter.GeoGridCoordinate2D 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.adapter;
import ucar.nc2.dataset.CoordinateAxis2D;
import ucar.ma2.ArrayDouble;
import ucar.ma2.MAMath;
/**
* fork ucar.nc2.dt.grid.GridCoordinate2D for adaption of GridCoverage
*
* 2D Coordinate System has lat(x,y) and lon(x,y).
* This class implements finding the index (i,j) from (lat, lon) coord.
* This is for "one-off" computation, not a systematic lookup table for all points in a pixel array.
* Hueristic search of the 2D space for the cell that contains the point.
*
* @author caron
* @since Jul 10, 2009
*/
public class GeoGridCoordinate2D {
private static boolean debug;
private static final org.slf4j.Logger log = org.slf4j.LoggerFactory.getLogger(GeoGridCoordinate2D.class);
private final CoordinateAxis2D latCoord, lonCoord;
private final int nrows, ncols;
private ArrayDouble.D2 latEdge, lonEdge;
private MAMath.MinMax latMinMax, lonMinMax;
GeoGridCoordinate2D(CoordinateAxis2D latCoord, CoordinateAxis2D lonCoord) {
this.latCoord = latCoord;
this.lonCoord = lonCoord;
assert latCoord.getRank() == 2;
assert lonCoord.getRank() == 2;
int[] shape = latCoord.getShape();
nrows = shape[0];
ncols = shape[1];
}
private void findBounds() {
if (lonMinMax != null)
return;
lonEdge = lonCoord.getEdges();
latEdge = latCoord.getEdges();
// assume missing values have been converted to NaNs
latMinMax = MAMath.getMinMax(latEdge);
lonMinMax = MAMath.getMinMax(lonEdge);
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);
}
// brute force
public boolean findCoordElementForce(double wantLat, double wantLon, int[] rectIndex) {
findBounds();
if (wantLat < latMinMax.min)
return false;
if (wantLat > latMinMax.max)
return false;
if (wantLon < lonMinMax.min)
return false;
if (wantLon > lonMinMax.max)
return false;
boolean saveDebug = debug;
debug = 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)) {
debug = saveDebug;
return true;
}
}
}
// debug = saveDebug;
return false;
}
public boolean findCoordElement(double wantLat, double wantLon, int[] rectIndex) {
return findCoordElementNoForce(wantLat, wantLon, rectIndex);
}
/**
* 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 findCoordElementNoForce(double wantLat, double wantLon, int[] rectIndex) {
findBounds();
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);
}
}
}
/**
* 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;
return !(wantLon > lonEdge.get(row, col + 1));
}
/**
* 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 oo 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;
return sign == detIsPositive(x4, y4, x1, y1, wantLon, wantLat);
}
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=" + latCoord.getFullName() + ", " + latCoord.getFullName()); // 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]);
return (row != rectIndex[0]) || (col != rectIndex[1]); // nothing has changed
}
}
/**
* 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]);
return (row != rectIndex[0]) || (col != rectIndex[1]); // nothing has changed
}
}
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;
}
}
© 2015 - 2025 Weber Informatics LLC | Privacy Policy