cdm/src/main/java/ucar/nc2/ft2/coverage/adapter/GeoGridCoordinate2D.java
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 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;
  }
/*
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.