ucar.nc2.dt.grid.GridCoordinate2D 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.dt.grid;
import ucar.nc2.dataset.CoordinateAxis2D;
import ucar.nc2.dataset.CoordinateAxis;
import ucar.nc2.dt.GridCoordSystem;
import ucar.ma2.ArrayDouble;
import ucar.ma2.MAMath;
import java.io.IOException;
/**
* 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 GridCoordinate2D {
static private boolean debug = false;
static private final org.slf4j.Logger log = org.slf4j.LoggerFactory.getLogger(GridCoordinate2D.class);
static final String testfilename = "Q:/cdmUnitTest/ft/fmrc/rtofs/ofs.20091122/ofs_atl.t00z.F024.grb.grib2";
private final CoordinateAxis2D latCoord, lonCoord;
private final int nrows, ncols;
private ArrayDouble.D2 latEdge, lonEdge;
private MAMath.MinMax latMinMax, lonMinMax;
GridCoordinate2D(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 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 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;
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=" + 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]);
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;
}
private static void doOne(GridCoordinate2D g2d, double wantLat, double wantLon) {
int[] result = new int[2];
if (g2d.findCoordElementForce(wantLat, wantLon, result))
System.out.printf("Brute (%f %f) == (%d %d) %n", wantLat, wantLon, result[0], result[1]);
else {
System.out.printf("Brute (%f %f) FAIL", wantLat, wantLon);
return;
}
if (g2d.findCoordElement(wantLat, wantLon, result))
System.out.printf("(%f %f) == (%d %d) %n", wantLat, wantLon, result[0], result[1]);
else
System.out.printf("(%f %f) FAIL %n", wantLat, wantLon);
System.out.printf("----------------------------------------%n");
}
public static void test1() throws IOException {
String filename = "D:/work/asaScience/EGM200_3.ncml";
GridDataset gds = GridDataset.open(filename);
GeoGrid grid = gds.findGridByName("u_wind");
GridCoordSystem gcs = grid.getCoordinateSystem();
CoordinateAxis lonAxis = gcs.getXHorizAxis();
assert lonAxis instanceof CoordinateAxis2D;
CoordinateAxis latAxis = gcs.getYHorizAxis();
assert latAxis instanceof CoordinateAxis2D;
GridCoordinate2D g2d = new GridCoordinate2D((CoordinateAxis2D) latAxis, (CoordinateAxis2D) lonAxis);
doOne(g2d, 35.0, -6.0);
doOne(g2d, 34.667302, -5.008376); // FAIL
doOne(g2d, 34.667303, -6.394240);
doOne(g2d, 36.6346, -5.0084);
doOne(g2d, 36.6346, -6.394240);
gds.close();
}
public static void test2() throws IOException {
String filename = "C:/data/fmrc/apex_fmrc/Run_20091025_0000.nc";
GridDataset gds = GridDataset.open(filename);
GeoGrid grid = gds.findGridByName("temp");
GridCoordSystem gcs = grid.getCoordinateSystem();
CoordinateAxis lonAxis = gcs.getXHorizAxis();
assert lonAxis instanceof CoordinateAxis2D;
CoordinateAxis latAxis = gcs.getYHorizAxis();
assert latAxis instanceof CoordinateAxis2D;
GridCoordinate2D g2d = new GridCoordinate2D((CoordinateAxis2D) latAxis, (CoordinateAxis2D) lonAxis);
doOne(g2d, 40.166959,-73.954234);
gds.close();
}
public static void test3() throws IOException {
String filename = testfilename;
GridDataset gds = GridDataset.open(filename);
GeoGrid grid = gds.findGridByName("Sea_Surface_Height_Relative_to_Geoid_surface");
GridCoordSystem gcs = grid.getCoordinateSystem();
CoordinateAxis lonAxis = gcs.getXHorizAxis();
assert lonAxis instanceof CoordinateAxis2D;
CoordinateAxis latAxis = gcs.getYHorizAxis();
assert latAxis instanceof CoordinateAxis2D;
GridCoordinate2D g2d = new GridCoordinate2D((CoordinateAxis2D) latAxis, (CoordinateAxis2D) lonAxis);
doOne(g2d, -15.554099426977835, -0.7742870290336263);
gds.close();
}
public static void main(String[] args) throws IOException {
debug = true;
test3();
}
}