org.integratedmodelling.engine.geospace.extents.Grid Maven / Gradle / Ivy
The newest version!
/*******************************************************************************
* Copyright (C) 2007, 2015:
*
* - Ferdinando Villa
* - integratedmodelling.org
* - any other authors listed in @author annotations
*
* All rights reserved. This file is part of the k.LAB software suite,
* meant to enable modular, collaborative, integrated
* development of interoperable data and model components. For
* details, see http://integratedmodelling.org.
*
* This program is free software; you can redistribute it and/or
* modify it under the terms of the Affero General Public License
* Version 3 or any later version.
*
* This program is distributed in the hope that it will be useful,
* but without any warranty; without even the implied warranty of
* merchantability or fitness for a particular purpose. See the
* Affero General Public License for more details.
*
* You should have received a copy of the Affero General Public License
* along with this program; if not, write to the Free Software
* Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
* The license is also available at: https://www.gnu.org/licenses/agpl.html
*******************************************************************************/
package org.integratedmodelling.engine.geospace.extents;
import java.util.ArrayList;
import java.util.Collection;
import java.util.Iterator;
import org.geotools.coverage.grid.GeneralGridEnvelope;
import org.geotools.geometry.jts.ReferencedEnvelope;
import org.integratedmodelling.api.modelling.IExtent;
import org.integratedmodelling.api.space.IGrid;
import org.integratedmodelling.api.space.IGridMask;
import org.integratedmodelling.api.space.IShape;
import org.integratedmodelling.collections.Pair;
import org.integratedmodelling.common.space.IGeometricShape;
import org.integratedmodelling.common.space.SpaceLocator;
import org.integratedmodelling.common.utils.MiscUtilities;
import org.integratedmodelling.common.vocabulary.Unit;
import org.integratedmodelling.engine.geospace.Geospace;
import org.integratedmodelling.engine.geospace.coverage.raster.RasterActivationLayer;
import org.integratedmodelling.engine.geospace.gis.ThinklabRasterizer;
import org.integratedmodelling.engine.geospace.literals.ShapeValue;
import org.integratedmodelling.exceptions.KlabException;
import org.integratedmodelling.exceptions.KlabRuntimeException;
import org.integratedmodelling.exceptions.KlabUnsupportedOperationException;
import org.integratedmodelling.exceptions.KlabValidationException;
import org.opengis.coverage.grid.GridEnvelope;
import org.opengis.referencing.crs.CoordinateReferenceSystem;
import com.vividsolutions.jts.geom.Envelope;
import com.vividsolutions.jts.geom.Geometry;
import com.vividsolutions.jts.geom.GeometryFactory;
/**
* Cell extents are the extents generated by the Raster conceptual models. Vector
* CMs should be prepared to deal with them too.
*
* A GridExtent may have been defined from a ShapeExtent by rasterizing objects. For this
* reason, we implement ILineageTraceable so that the original objects can be reconstructed.
* This must be transferred to the conceptualized grid in some (preferably elegant) way (TODO).
*
* Provides storage metadata that we float to the models and contexts that contain it, so that
* they can be queried for coverage with another context as parameter.
*
* @author Ferdinando, Ioannis
*
*/
public class Grid extends Area implements IGrid {
/*
* only set in a subgrid of the grid having the parent ID.
*/
int[] offsetInSupergrid = null;
String superGridId = null;
GeometryFactory gFactory = null;
int xDivs = 0;
int yDivs = 0;
double cellLength = 0.0;
double cellHeight = 0.0;
double xOrigin = 0.0;
double yOrigin = 0.0;
IGridMask activationLayer = null;
private double cellHeightMeters;
private double cellWidthMeters;
private double cellAreaMeters = -1.0;
private CoordinateReferenceSystem metersCRS = null;
Geometry boundary = null;
public ShapeValue shape;
public class CellImpl implements Cell {
int x;
int y;
CellImpl(int x, int y) {
this.x = x;
this.y = y;
}
@Override
public int getX() {
return x;
}
@Override
public int getY() {
return y;
}
@Override
public Cell N() {
int trow = getYCells() - 1;
if (y < trow) {
return new CellImpl(x, y + 1);
}
return null;
}
@Override
public Cell S() {
if (y > 0) {
return new CellImpl(x, y - 1);
}
return null;
}
@Override
public Cell E() {
if (x > 0) {
return new CellImpl(x - 1, y);
}
return null;
}
@Override
public Cell W() {
int tcol = getXCells() - 1;
if (x < tcol) {
return new CellImpl(x + 1, y);
}
return null;
}
@Override
public Cell NW() {
int trow = getYCells() - 1, tcol = getXCells() - 1;
if (y < trow && x < tcol) {
return new CellImpl(x + 1, y + 1);
}
return null;
}
@Override
public Cell NE() {
int trow = getYCells() - 1;
if (y < trow && x > 0) {
return new CellImpl(x - 1, y + 1);
}
return null;
}
@Override
public Cell SE() {
if (y > 0 && x > 0) {
return new CellImpl(x - 1, y - 1);
}
return null;
}
@Override
public Cell SW() {
int tcol = getXCells() - 1;
if (y > 0 && x < tcol) {
return new CellImpl(x + 1, y - 1);
}
return null;
}
/*
* return all the neighboring cells that actually exist in
* the grid.
*/
@Override
public Collection getNeighbors() {
ArrayList ret = new ArrayList<>();
Cell c = N();
if (c != null)
ret.add(c);
c = S();
if (c != null)
ret.add(c);
c = W();
if (c != null)
ret.add(c);
c = E();
if (c != null)
ret.add(c);
c = NW();
if (c != null)
ret.add(c);
c = NE();
if (c != null)
ret.add(c);
c = SW();
if (c != null)
ret.add(c);
c = SE();
if (c != null)
ret.add(c);
return ret;
}
@Override
public boolean equals(Object o) {
return o instanceof CellImpl && ((CellImpl) o).x == x && ((CellImpl) o).y == y;
}
@Override
public int hashCode() {
return (x + ":" + y).hashCode();
}
@Override
public Cell move(int xOfs, int yOfs) {
int newX = x + xOfs;
int newY = y + yOfs;
if (newX >= 0 && newX < xDivs) {
if (newY >= 0 && newY < yDivs) {
// if (isActive(newX, newY)) {
return new CellImpl(newX, newY);
// }
}
}
return null;
}
@Override
public double getMinX() {
return Grid.this.getMinX() + (x * getCellWidth());
}
@Override
public double getMaxX() {
return this.getMinX() + getCellWidth();
}
@Override
public double getMinY() {
return Grid.this.getMinY() + ((yDivs - y - 1) * getCellHeight());
}
@Override
public double getMaxY() {
return this.getMinY() + getCellHeight();
}
@Override
public Integer getOffsetInGrid() {
return getOffset(x, y);
}
@Override
public double[] getCenter() {
return getCoordinates(getOffsetInGrid());
}
public Geometry getEnvelope() {
return ShapeValue.makeCell(getMinX(), getMinY(), getMaxX(), getMaxY());
}
public Geometry getGridEnvelope() {
return getBoundary();
}
public Grid getGrid() {
return Grid.this;
}
@Override
public IShape getGeometry() {
return new ShapeValue(getEnvelope(), crs);
}
}
public Grid() {
/*
* for conceptualization - make sure you use this properly
*/
}
@Override
public SpaceLocator getLocator(int x, int y) {
return new SpaceLocator(x, yDivs - y - 1);
}
public Grid(CoordinateReferenceSystem crs, double x1, // lonLowerBound
double y1, // latLowerBound
double x2, // lonUpperBound
double y2, // latUpperBound
int xDivs, int yDivs) {
// super(crs, x1, y1, x2, y2);
this.xOrigin = x1;
this.yOrigin = y1;
this.crs = crs;
this.envelope = new ReferencedEnvelope(x1, x2, y1, y2, crs);
setResolution(xDivs, yDivs);
}
public Grid(Grid gridExtent) {
super(gridExtent.getCRS(), gridExtent.getWest(), gridExtent.getSouth(), gridExtent
.getEast(), gridExtent.getNorth());
this.xOrigin = gridExtent.getEast();
this.yOrigin = gridExtent.getSouth();
this.setResolution(gridExtent.getXCells(), gridExtent.getYCells());
}
/**
* Create a grid from a shape in the CRS of the shape and using the given
* resolution for the larger extent.
*
* @param shape
* @param resolutionInMeters
* @throws KlabException
*/
public Grid(ShapeValue shape, String resolutionInMeters) throws KlabException {
super(shape);
ShapeValue bb = shape.getBoundingBox().convertToMeters();
double meters = parseResolution(resolutionInMeters);
int linearResolution = (int) (bb.getEnvelope().getWidth() / meters);
setAdjustedEnvelope(shape, linearResolution);
}
/**
* Parse a string like "1 km" and return the meters in it. Throw an exception if this
* cannot be parsed.
*
* @param string
* @return the resolution in meters
* @throws KlabValidationException
*/
public static double parseResolution(String string) throws KlabValidationException {
Pair pd = MiscUtilities.splitNumberFromString(string);
if (pd.getFirst() == null || pd.getSecond() == null)
throw new KlabValidationException("wrong resolution specification: " + string);
Unit uu = new Unit(pd.getSecond());
Unit mm = new Unit("m");
return mm.convert(pd.getFirst().doubleValue(), uu).doubleValue();
}
public Grid(ShapeValue shape, double resolutionInMeters) throws KlabException {
super(shape);
// ShapeValue sh = shape.transform(Geospace.get().getDefaultCRS());
// ReferencedEnvelope bb = sh.getBoundingBox().getEnvelope();
//
// double x1 = GeoNS.mercX(bb.getMinX());
// double x2 = GeoNS.mercX(bb.getMaxX());
//
// double width = (x2 > x1) ? (x2 - x1) : (x1 - x2);
//
// int linearResolution = (int) (width / resolutionInMeters);
setAdjustedEnvelope(shape, resolutionInMeters);
}
// public Grid(ShapeValue shape, double resolutionInMeters) throws ThinklabException {
// super(shape);
//
// ShapeValue sh = shape.transform(Geospace.get().getDefaultCRS());
// Envelope enve = sh.getEnvelope();
//
// // double[] bl = GeoNS.LatLonToUTMXY(enve.getMinY(), enve.getMinX());
// // double[] tr = GeoNS.LatLonToUTMXY(enve.getMaxY(), enve.getMaxX());
//
// double distX = GeoNS.getDistance(new Coordinate(enve.getMinY(), enve.getMaxX()), new
// Coordinate(enve.getMinY(), enve.getMinX())) * 1000;
// double distY = GeoNS.getDistance(new Coordinate(enve.getMaxY(), enve.getMinX()), new
// Coordinate(enve.getMinY(), enve.getMinX())) * 1000;
// // double distX = tr[0] - bl[0];
// // double distY = tr[1] - bl[1];
//
// int xcells = (int)Math.ceil(distX/resolutionInMeters);
// int ycells = (int)Math.ceil(distY/resolutionInMeters);
//
// /*
// * delta to add to max x,y
// */
// double dfacX = 1.0 - (xcells*resolutionInMeters)/distX;
// double dfacY = 1.0 - (ycells*resolutionInMeters)/distY;
//
// Envelope oEnv = shape.getEnvelope();
// double dx = (oEnv.getMaxX() - oEnv.getMinX()) * (dfacX/2);
// double dy = (oEnv.getMaxY() - oEnv.getMinY()) * (dfacY/2);
//
// this.shape = shape;
// this.envelope = new ReferencedEnvelope(oEnv.getMinX() - dx, oEnv.getMaxX() + dx,
// oEnv.getMinY() - dy, oEnv.getMaxY() + dy, shape.getCRS());
// setResolution(xcells, ycells);
// activationLayer = ThinklabRasterizer.createMask(shape, this);
// this.cellWidthMeters = resolutionInMeters;
// this.cellHeightMeters = resolutionInMeters;
// this.cellAreaMeters = resolutionInMeters*resolutionInMeters;
// }
public Grid transform(CoordinateReferenceSystem crs) throws KlabException {
if (shape.getCRS().equals(crs))
return this;
ShapeValue shape = getShape();
shape = shape.transform(crs);
return new Grid(shape, getXCells(), getYCells());
}
/**
* Set the envelope from shape.
* @param shape could be in any CRS.
* @param squareSize should be in meters
* @throws KlabException
*/
private void setAdjustedEnvelope(ShapeValue shape, double squareSize) throws KlabException {
int x = 0, y = 0;
double dx = 0, dy = 0;
ReferencedEnvelope env = shape.getEnvelope();
// Case one: both CRS and square size are in meters.
if (shape.getCRS().getCoordinateSystem().getAxis(0).getUnit().toString().equals("m")) {
double height = env.getHeight();
double width = env.getWidth();
x = (int) Math.ceil(width / squareSize);
y = (int) Math.ceil(height / squareSize);
dx = (x * squareSize) - width;
dy = (y * squareSize) - height;
ReferencedEnvelope envelop_ = new ReferencedEnvelope(env.getMinX() - dx, env.getMaxX()
+ dx, env.getMinY() - dy, env.getMaxY() + dy, shape.getCRS());
// After doing the calculations in the original CRS,
// shape and envelope is transformed to default.
this.shape = shape.transform(Geospace.get().getDefaultCRS());
try {
this.envelope = envelop_.transform(Geospace.get().getDefaultCRS(), true);
} catch (Exception e) {
// TODO: Shouldn't happen
e.printStackTrace();
}
this.crs = Geospace.get().getDefaultCRS();
}
// Case 2: CRS uses degrees
else {
// lets get height and width in meters
double height = Grid.haversine(env.getMinY(), env.getMinX(), env.getMaxY(), env.getMinX());
double width = Grid.haversine(env.getMinY(), env.getMinX(), env.getMinY(), env.getMaxX());
x = (int) Math.ceil(width / squareSize);
y = (int) Math.ceil(height / squareSize);
// Here I tried to adjust further dx and dy based on the size of the grid cell
// at the center of the Grid. Possibly not right either, so I left it out.
// double centralX = (env.getMaxX() + env.getMinX())/2;
// double stepX = (env.getMaxX() - env.getMinX())/x;
//
// double centralY = (env.getMaxY() + env.getMinY())/2;
// double stepY = (env.getMaxY() - env.getMinY())/x;
//
// double actualStepWidth = Grid.haversine(centralY,centralX,
// centralY,centralX+stepX);
//
// double actualStepHeight = Grid.haversine(centralY, centralX,
// centralY+stepY,centralX);
// dx = stepX * (squareSize/actualStepWidth) / 2;
// dy = stepY * (squareSize/actualStepHeight) / 2;
}
if (!this.crs.toString().equals(Geospace.get().getDefaultCRS().toString())) {
this.shape = this.shape.transform(Geospace.get().getDefaultCRS());
try {
this.envelope = this.envelope.transform(Geospace.get().getDefaultCRS(), true);
} catch (Exception e) {
// TODO: Shouldn't happen
e.printStackTrace();
}
this.crs = Geospace.get().getDefaultCRS();
}
this.setResolution(x, y);
activationLayer = ThinklabRasterizer.createMask(shape, this);
}
/**
* The haversine formula calculates great-circle distance between two points on a sphere
* from their longitudes and latitudes.
*
* From http://rosettacode.org/wiki/Haversine_formula#Java
*
* @param lat1 PointOne latitude
* @param lon1 PointOne longitude
* @param lat2 PointTwo latitude
* @param lon2 PointTwo longitude
* @return distance in meters
*/
public static double haversine(double lat1, double lon1, double lat2, double lon2) {
double R = 6372800; // in m
double dLat = Math.toRadians(lat2 - lat1);
double dLon = Math.toRadians(lon2 - lon1);
lat1 = Math.toRadians(lat1);
lat2 = Math.toRadians(lat2);
double a = Math.sin(dLat / 2) * Math.sin(dLat / 2) + Math.sin(dLon / 2) * Math.sin(dLon / 2)
* Math.cos(lat1) * Math.cos(lat2);
double c = 2 * Math.asin(Math.sqrt(a));
return R * c;
}
// /*
// * set envelope from shape's bounding box, adjusted to keep cells square when
// * measured in meters.
// */
// private void setAdjustedEnvelope(ShapeValue shape, int majorAxisResolution) throws ThinklabException {
//
// int x = 0, y = 0;
// double dx = 0, dy = 0;
// ReferencedEnvelope env = shape.getEnvelope();
//
// if (Geospace.get().squareCellsM()) {
// //
// // CoordinateReferenceSystem meters = Geospace.get().getMetersCRS();
// // if (Geospace.getCRSIdentifier(shape.getCRS(), false).equals("EPSG:4326")) {
// // meters = Geospace.get().getMetersCRS(env.getMinX() + env.getWidth() / 2,
// // env.getMinY() + env.getHeight() / 2);
// // }
// //
// // try {
// // env = env.transform(meters, true);
// // } catch (Exception e) {
// // throw new ThinklabValidationException(e);
// // }
// //
// // double width = env.getWidth();
// // double height = env.getHeight();
// //
// // if (width > height) {
// // double minorAxisResolution = Math.ceil(majorAxisResolution * height / width);
// // height = width * minorAxisResolution / majorAxisResolution;
// // x = majorAxisResolution;
// // y = (int) minorAxisResolution;
// // dy = (height - env.getHeight()) / 2.0;
// // } else {
// // double minorAxisResolution = Math.ceil(majorAxisResolution * width / height);
// // width = height * minorAxisResolution / majorAxisResolution;
// // x = (int) minorAxisResolution;
// // y = majorAxisResolution;
// // dx = (width - env.getWidth()) / 2.0;
// // }
// //
// // // System.out.println("ORIG-U: " + shape.getEnvelope());
// // // System.out.println("ORIG-M: " + env);
// //
// // this.shape = shape;
// // env = new ReferencedEnvelope(env.getMinX() - dx, env.getMaxX() + dx, env.getMinY() - dy,
// // env.getMaxY() + dy, meters);
// //
// // try {
// // this.envelope = env.transform(shape.getCRS(), true);
// // } catch (Exception e) {
// // throw new ThinklabValidationException(e);
// // }
//
// } else {
//
// double width = env.getWidth();
// double height = env.getHeight();
//
// if (width > height) {
// double minorAxisResolution = Math.ceil(majorAxisResolution * height / width);
// height = width * minorAxisResolution / majorAxisResolution;
// x = majorAxisResolution;
// y = (int) minorAxisResolution;
// dy = (height - env.getHeight()) / 2.0;
// } else {
// double minorAxisResolution = Math.ceil(majorAxisResolution * width / height);
// width = height * minorAxisResolution / majorAxisResolution;
// x = (int) minorAxisResolution;
// y = majorAxisResolution;
// dx = (width - env.getWidth()) / 2.0;
// }
//
// this.shape = shape;
// this.envelope = new ReferencedEnvelope(env.getMinX() - dx, env.getMaxX() + dx,
// env.getMinY() - dy, env.getMaxY() + dy, shape.getCRS());
// }
//
// // long xx = (long) x * (long) y;
// // if (xx > 2999999) {
// // throw new ThinklabValidationException(xx + " cells: please reduce resolution");
// // }
//
// this.setResolution(x, y);
//
// activationLayer = ThinklabRasterizer.createMask(shape, this);
//
// }
/**
* Create a grid from a shape in the CRS of the shape and using the given
* resolution for the larger extent.
*
* @param shape
* @param linearResolution
* @throws KlabException
*/
public Grid(ShapeValue shape, int linearResolution) throws KlabException {
super(shape);
setAdjustedEnvelope(shape, linearResolution);
}
/**
* Create a grid extent with a rasterized shape in the given grid. Will not
* clip the shape or anything: use ONLY when you already know the precise aspect factor for the extent
* resulting from the shape.
*
* @param shape
* @param x
* @param y
* @throws KlabException
*/
public Grid(ShapeValue shape, int x, int y) throws KlabException {
super(shape);
this.setResolution(x, y);
activationLayer = ThinklabRasterizer.createMask(shape, this);
}
/*
* TEMPORARY - does not create the activation layer; the parameter is a trick. Activation layers
* need to be small and smart - this is to avoid creating a monster for the time being.
*
* Create a grid extent with a rasterized shape in the given grid. Will not
* clip the shape or anything: use ONLY when you already know the precise aspect factor for the extent
* resulting from the shape.
*
* @param shape
* @param x
* @param y
* @throws ThinklabException
*/
@Deprecated
public Grid(ShapeValue shape, int x, int y, boolean dummy) throws KlabException {
super(shape);
this.setResolution(x, y);
}
// /**
// * Used to set the activation layer from a shape. This also create a ShapeExtent and
// * sets the lineage to it, so that our "object" can be reconstructed.
// *
// * @param extent
// * @param shape
// * @throws ThinklabException
// */
// public Grid(Grid extent, ShapeValue shape) throws ThinklabException {
// this(extent);
// activationLayer = ThinklabRasterizer.createMask(shape, this);
// }
public void setResolution(int xDivs, int yDivs) {
this.xDivs = xDivs;
this.yDivs = yDivs;
cellLength = getEnvelope().getWidth() / xDivs;
cellHeight = getEnvelope().getHeight() / yDivs;
}
@Override
public int[] getXYOffsets(int index) {
int xx = index % getXCells();
int yy = getYCells() - (index / getXCells()) - 1;
return new int[] { xx, yy };
}
public static int[] getXYCoordinates(int index, int width, int height) {
int xx = index % width;
int yy = /* height - (*/index / width/*) - 1*/;
return new int[] { xx, yy };
}
@Override
public int getOffset(int x, int y) {
return ((yDivs - y - 1) * getXCells()) + x;
}
/**
* Activate the cells that correspond to the passed shape. Note that this
* doesn't change the multiplicity.
*
* @param roi
* @throws KlabException
*/
public void createActivationLayer(ShapeValue roi) throws KlabException {
this.activationLayer = ThinklabRasterizer.createMask(roi, this);
this.shape = roi;
}
@Override
public boolean equals(Object obj) {
if (obj instanceof Grid) {
Grid ge = (Grid) obj;
/*
* TODO compare activation layers
*/
return xDivs == ge.xDivs && yDivs == ge.yDivs && envelope.equals(ge.envelope);
}
return false;
}
private void checkMeters() {
if (cellAreaMeters < -0.1) {
try {
/*
* TODO bullcrap
*/
metersCRS = Geospace.get().getMetersCRS();
if (Geospace.getCRSIdentifier(getCRS(), false).equals("EPSG:4326")) {
metersCRS = Geospace.get().getMetersCRS(envelope.getMinX()
+ envelope.getWidth() / 2, envelope.getMinY() + envelope.getHeight() / 2);
}
Envelope env = this.envelope.transform(metersCRS, true);
this.cellWidthMeters = env.getWidth() / getXCells();
this.cellHeightMeters = env.getHeight() / getYCells();
this.cellAreaMeters = this.cellWidthMeters * this.cellHeightMeters;
} catch (Exception e) {
throw new KlabRuntimeException(e);
}
}
}
@Override
public double getCellWidth() {
return cellLength;
}
@Override
public double getCellHeight() {
return cellHeight;
}
public double getCellWidthMeters() {
checkMeters();
return this.cellWidthMeters;
}
public double getCellHeightMeters() {
checkMeters();
return this.cellHeightMeters;
}
public double getCellAreaMeters() {
checkMeters();
return this.cellAreaMeters;
}
public double getTotalAreaSquareMeters() {
double ret = getCellAreaMeters();
if (activationLayer != null) {
ret *= activationLayer.totalActiveCells();
} else {
ret *= (xDivs * yDivs);
}
return ret;
}
/**
* Return a box in standard AWT coordinates. Useful for some graphical ops and rasterization.
*
* @return box
*/
public java.awt.geom.Rectangle2D.Double getDefaultBox() {
ReferencedEnvelope env = getEnvelope();
return new java.awt.geom.Rectangle2D.Double(env.getMinX(), env.getMinY(), env.getWidth(), env
.getHeight());
}
@Override
public ShapeValue getShape() {
return (ShapeValue) (shape == null ? getFullExtentValue() : shape);
}
public java.awt.geom.Rectangle2D.Double getNormalizedBox() {
ReferencedEnvelope env = getEnvelope();
return new java.awt.geom.Rectangle2D.Double(env.getMinX(), env.getMinY(), env.getWidth(), env
.getHeight());
}
/**
* return the envelope of the cell at x,y, irrespective of the activation layer
*/
public ReferencedEnvelope getCellEnvelope(int x, int y) {
double x1 = xOrigin + (cellLength * x);
double y1 = yOrigin + (cellHeight * y);
double x2 = x1 + cellLength;
double y2 = y1 + cellHeight;
return new ReferencedEnvelope(x1, x2, y1, y2, getCRS());
}
@Override
public double getCellArea() {
return getCellAreaMeters();
// Envelope env = getCellEnvelope(0, 0);
// return env.getHeight() * env.getWidth();
}
@Override
public String toString() {
return "grid(" + xDivs + "," + yDivs + ": " + envelope.getMinX() + " : " + envelope.getMaxX() + ", "
+ envelope.getMinY() + " : " + envelope.getMaxY() + " " + getCRS().getName() + ")";
}
public IGeometricShape getFullExtentValue() {
return new ShapeValue(getBoundary(), crs);
}
public ShapeValue getCellPolygon(int index) {
Cell cell = getCell(index);
return (ShapeValue) cell.getGeometry();
// /*
// * determine coordinates of granule. This should not get called if the activation layer
// * is null, so no check for now.
// */
// Pair xy = requireActivationLayer(true).getCell(index);
//
// /*
// * TODO reimplement to use nextCell on activation layer. Must enforce sequential access
// * to use effectively.
// */
//
// double x1 = envelope.getMinX() + (cellLength * xy.getFirst());
// double y1 = envelope.getMinY() + (cellHeight * xy.getSecond());
// double x2 = x1 + cellLength;
// double y2 = y1 + cellHeight;
//
// return new ShapeValue(x1, y1, x2, y2, crs);
}
public Cell getCell(int index) {
int[] xy = getXYOffsets(index);
return new CellImpl(xy[0], xy[1]);
}
public Cell getCell(int x, int y) {
return new CellImpl(x, y);
}
@Override
public int getCellCount() {
return xDivs * yDivs;
}
public Geometry getBoundary() {
if (boundary == null) {
ReferencedEnvelope env = getEnvelope();
boundary = ShapeValue.makeCell(env.getMinX(), env.getMinY(), env.getMaxX(), env.getMaxY());
}
return boundary;
}
@Override
public int getXCells() {
return xDivs;
}
@Override
public int getYCells() {
return yDivs;
}
public int getXMaxCell() {
return xDivs;
}
public int getXMinCell() {
return 0;
}
public int getYMinCell() {
return 0;
}
public int getYMaxCell() {
return yDivs;
}
public GridEnvelope getGridRange() {
return new GeneralGridEnvelope(new int[] { 0, 0 }, new int[] { xDivs, yDivs }, false);
}
public IGridMask getActivationLayer() throws KlabValidationException {
if (activationLayer == null) {
throw new KlabValidationException("no activation layer in grid extent");
}
return activationLayer;
}
public IGridMask requireActivationLayer(boolean active) {
if (activationLayer == null) {
activationLayer = new RasterActivationLayer(xDivs, yDivs, active, this);
((RasterActivationLayer)activationLayer).createMask();
}
return activationLayer;
}
public double getNSResolution() {
return cellHeight;
}
public double getEWResolution() {
return cellLength;
}
/*
* TODO check: this should be creating the largest extent that preserves
* the topology of the original extent and covers the other extent. Cell size should
* be the same as in the original extent, and the envelope should be the same or
* smaller.
*
* @see org.integratedmodelling.geospace.extents.ArealExtent#createMergedExtent(org.integratedmodelling.geospace.extents.ArealExtent, org.integratedmodelling.geospace.extents.ArealExtent, org.opengis.referencing.crs.CoordinateReferenceSystem, com.vividsolutions.jts.geom.Envelope, com.vividsolutions.jts.geom.Envelope, com.vividsolutions.jts.geom.Envelope)
*/
@Override
protected Area createMergedExtent(Area orextent, Area otextent, CoordinateReferenceSystem ccr, Envelope common, Envelope orenvnorm, Envelope otenvnorm)
throws KlabException {
/*
* for now, raster always wins
*/
if (otextent instanceof Grid && !(orextent instanceof Grid)) {
return makeRasterExtent((Grid) otextent, orextent, ccr, common);
}
if (orextent instanceof Grid && !(otextent instanceof Grid)) {
return makeRasterExtent((Grid) orextent, otextent, ccr, common);
}
/*
* if we get here, we must be merging two rasters
*/
if (!(orextent instanceof Grid && otextent instanceof Grid)) {
throw new KlabUnsupportedOperationException("RasterModel: cannot yet merge extents of different types");
}
Grid orext = (Grid) orextent;
Grid otext = (Grid) otextent;
/* we'll fix the resolution later */
Grid nwext = new Grid(ccr, common.getMinX(), common.getMinY(), common.getMaxX(), common
.getMaxY(), 1, 1);
// this is the constrained resolution; we recompute it below if we are free to choose
int xc = otext.getXCells();
int yc = otext.getYCells();
// phase errors in both directions due to choosing resolutions that do not
// match exactly. This is computed but not used right now. We could set what to
// do with it as a property: e.g., log a warning or even raise an exception if nonzero.
double errx = 0.0;
double erry = 0.0;
/* choose the smallest of the reprojected cells unless we're constrained to accept otextent */
Envelope cor = null;
Envelope cot = null;
try {
cor = orext.getCellEnvelope(0, 0).transform(ccr, true, 10);
cot = otext.getCellEnvelope(0, 0).transform(ccr, true, 10);
} catch (Exception e) {
throw new KlabValidationException(e);
}
double aor = cor.getHeight() * cor.getWidth();
double aot = cot.getHeight() * cot.getWidth();
/*
* We take the finest res
*/
Envelope cell = aor < aot ? cor : cot;
// System.out.println("cells are " + cor + " and " + cot + "; chosen " + cell + " because areas are "
// + aor + " and " + aot);
/* recompute the number of cells in the new extent */
xc = (int) Math.round(nwext.getEnvelope().getWidth() / cell.getWidth());
yc = (int) Math.round(nwext.getEnvelope().getHeight() / cell.getHeight());
errx = nwext.getEnvelope().getWidth() - (cell.getWidth() * xc);
erry = nwext.getEnvelope().getHeight() - (cell.getHeight() * yc);
// System.out.println("new cell size is " + xc + "," + yc);
// TODO use the error, make sure it's not larger than too much
// System.out.println("errors are " + errx + "," + erry);
nwext.setResolution(xc, yc);
System.out.println("extent is now " + nwext);
return nwext;
}
private Area makeRasterExtent(Grid grid, Area otextent, CoordinateReferenceSystem ccr, Envelope env)
throws KlabException {
/*
* TODO adjust the extent as necessary to make the grid extent reflect the coordinate system and area
*/
Grid ret = new Grid(ccr, env.getMinX(), env.getMinY(), env.getMaxX(), env.getMaxY(), grid
.getXCells(), grid.getYCells());
return ret;
}
@Override
protected Area createConstrainedExtent(Area orextent, Area otextent, CoordinateReferenceSystem ccr, Envelope common, Envelope orenvnorm, Envelope otenvnorm)
throws KlabException {
/*
* for now, raster always wins
*/
if (otextent instanceof Grid && !(orextent instanceof Grid)) {
return makeRasterExtent((Grid) otextent, orextent, ccr, common);
}
if (orextent instanceof Grid && !(otextent instanceof Grid)) {
return makeRasterExtent((Grid) orextent, otextent, ccr, common);
}
/*
* if we get here, we must be merging two rasters
*/
if (!(orextent instanceof Grid && otextent instanceof Grid)) {
throw new KlabUnsupportedOperationException("RasterModel: cannot yet merge extents of different types");
}
Grid otext = (Grid) otextent;
double xRatio = common.getWidth() / otenvnorm.getWidth();
double yRatio = common.getHeight() / otenvnorm.getHeight();
/* recompute the number of cells in the new extent */
int xc = (int) Math.round(otext.getXCells() * xRatio);
int yc = (int) Math.round(otext.getYCells() * yRatio);
/* we'll fix the resolution later */
Grid nwext = new Grid(ccr, common.getMinX(), common.getMinY(), common.getMaxX(), common
.getMaxY(), xc, yc);
System.out.println("constrained extent is now " + nwext);
return nwext;
}
public String getSignature() {
return "grid," + getWest() + "," + getSouth() + "," + getEast() + "," + getNorth() + ","
+ getXCells() + "," + getYCells() + Geospace.getCRSIdentifier(getCRS(), false);
}
// public IExtent getAggregatedExtent() {
// return new ShapeExtent((ShapeValue) getFullExtentValue());
// }
// @Override
// public Collection> getStateLocators(int index) {
//
// Pair xy = requireActivationLayer(true).getCell(index);
//
// int brow = 0, bcol = 0, trow = getYCells() -1, tcol = getXCells() - 1;
// int col = xy.getFirst(), row = xy.getSecond();
// int n = 0;
//
// ArrayList> ret = new ArrayList>();
//
// // n
// if (row < trow) {
// n = getIndex(col, row + 1);
// ret.add(new Pair("n", n));
// }
//
// // s
// if (row > brow) {
// n = getIndex(col, row - 1);
// ret.add(new Pair("s", n));
// }
//
// // e
// if (col > bcol) {
// n = getIndex(col - 1, row);
// ret.add(new Pair("e", n));
// }
//
// // w
// if (col < tcol) {
// n = getIndex(col + 1, row);
// ret.add(new Pair("w", n));
// }
//
// // nw
// if (row < trow && col < tcol) {
// n = getIndex(col + 1, row + 1);
// ret.add(new Pair("nw", n));
// }
//
// // ne
// if (row < trow && col > bcol) {
// n = getIndex(col - 1, row + 1);
// ret.add(new Pair("ne", n));
// }
//
// // se
// if (row > brow && col > bcol) {
// n = getIndex(col - 1, row - 1);
// ret.add(new Pair("se", n));
// }
//
// // sw
// if (row > brow && col < tcol) {
// n = getIndex(col + 1, row - 1);
// ret.add(new Pair("sw", n));
// }
//
// return ret;
// }
// public void setActivationLayer(IGridMask mask) {
// this.activationLayer = mask;
// }
// public boolean checkDomainDiscontinuity() throws ThinklabException {
// return false;
// }
public Grid intersection(IExtent extent) {
/*
* compute intersection of envelopes in our CSR. Ignore axis swap finally.
*/
ReferencedEnvelope ourenv = this.getEnvelope();
ReferencedEnvelope itsenv = null;
try {
itsenv = ((Area) extent).getEnvelope().transform(crs, true, 10);
} catch (Exception e) {
throw new KlabRuntimeException(e);
}
Envelope comenv = ourenv.intersection(itsenv);
if (comenv == null)
return null;
/*
* adjust intersection to be in phase with our cell size, ensuring it is fully contained.
*/
double startx = comenv.getMinX();
if (startx > ourenv.getMinX()) {
double delta = startx - ourenv.getMinX();
double nc = Math.floor(delta / cellLength);
startx = ourenv.getMinX() + (cellLength * nc);
}
double endx = comenv.getMaxX();
if (endx < ourenv.getMaxX()) {
double delta = ourenv.getMaxX() - endx;
double nc = Math.floor(delta / cellLength);
endx = ourenv.getMaxX() - (cellLength * nc);
}
double starty = comenv.getMinY();
if (starty > ourenv.getMinY()) {
double delta = starty - ourenv.getMinY();
double nc = Math.floor(delta / cellHeight);
starty = ourenv.getMinY() + (cellHeight * nc);
}
double endy = comenv.getMaxY();
if (endy < ourenv.getMaxY()) {
double delta = ourenv.getMaxY() - endy;
double nc = Math.floor(delta / cellHeight);
endx = ourenv.getMaxY() - (cellHeight * nc);
}
/*
* compute new grid extent resolution for intersection.
*/
int xc = (int) ((endx - startx) / cellLength);
int yc = (int) ((endy - starty) / cellHeight);
/*
* create and return new extent.
*/
return new Grid(crs, startx, starty, endx, endy, xc, yc);
}
public Collection getNeumannNeighbors(int xcell, int ycell) {
ArrayList ret = new ArrayList();
if (inRange(xcell - 1, ycell))
ret.add(getOffset(xcell - 1, ycell));
if (inRange(xcell + 1, ycell))
ret.add(getOffset(xcell + 1, ycell));
if (inRange(xcell, ycell - 1))
ret.add(getOffset(xcell, ycell - 1));
if (inRange(xcell, ycell + 1))
ret.add(getOffset(xcell, ycell + 1));
return ret;
}
public Collection getMooreNeighbors(int xcell, int ycell) {
ArrayList ret = (ArrayList) getNeumannNeighbors(xcell, ycell);
if (inRange(xcell - 1, ycell - 1))
ret.add(getOffset(xcell - 1, ycell - 1));
if (inRange(xcell + 1, ycell))
ret.add(getOffset(xcell + 1, ycell + 1));
if (inRange(xcell, ycell - 1))
ret.add(getOffset(xcell + 1, ycell - 1));
if (inRange(xcell, ycell + 1))
ret.add(getOffset(xcell - 1, ycell + 1));
return ret;
}
/**
* Get all grid cell indexes that are within n meters of the given point.
* @param xcell
* @param ycell
* @param meters
* @return index of cells within distance
*/
public Collection getCoordinatesWithinM(int xcell, int ycell, double meters) {
ArrayList ret = new ArrayList();
int xspan = (int) (meters / getCellWidthMeters());
int yspan = (int) (meters / getCellHeightMeters());
for (int x = 0; x < xspan; x++) {
for (int y = 0; y < yspan; y++) {
if (inRange(xcell - x, ycell - y))
ret.add(getOffset(xcell - x, ycell - y));
if (inRange(xcell + x, ycell - y))
ret.add(getOffset(xcell + x, ycell - y));
if (inRange(xcell - x, ycell + y))
ret.add(getOffset(xcell - x, ycell + y));
if (inRange(xcell + x, ycell + y))
ret.add(getOffset(xcell + x, ycell + y));
}
}
return ret;
}
public boolean inRange(int x, int y) {
return x >= 0 && x < getXCells() && y >= 0 && y < getYCells();
}
public boolean isCovered(int granule) {
if (activationLayer == null) {
return true;
}
return activationLayer.isActive(granule);
}
public double[] getCoordinatesAt(int x, int y) {
double x1 = getEnvelope().getMinX() + (cellLength * x);
double y1 = getEnvelope().getMinY() + (cellHeight * (getYCells() - y - 1));
return new double[] { x1 + (cellLength / 2), y1 + (cellHeight / 2) };
}
@Override
public double[] getCoordinates(int index) {
int[] xy = getXYOffsets(index);
return getCoordinatesAt(xy[0], xy[1]);
}
// @Override
// public Object adapt() {
// Map ret = new HashMap();
// ret.put("xdivs", xDivs);
// ret.put("ydivs", yDivs);
// if (shape != null) {
// ret.put("shape", shape.toString());
// }
// return ret;
// }
public int getIndex(int x, int y) {
if (xDivs != 0) {
return (y * xDivs) + x;
}
return 0;
}
public int[] getGridCoordinatesAt(double x, double y) {
int ofs = getOffsetFromWorldCoordinates(x, y);
return getXYOffsets(ofs);
}
@Override
public int getOffsetFromWorldCoordinates(double x, double y) {
if (x < getEnvelope().getMinX() || x > getEnvelope().getMaxX() || y < getEnvelope().getMinY()
|| y > getEnvelope().getMaxY())
return -1;
int xx = (int) (((x - getEnvelope().getMinX()) / (getEnvelope().getMaxX() - getEnvelope().getMinX()))
* xDivs);
int yy = (int) (((y - getEnvelope().getMinY()) / (getEnvelope().getMaxY() - getEnvelope().getMinY()))
* yDivs);
return getIndex(xx, yy);
}
@Override
public boolean isActive(int x, int y) {
return activationLayer == null ? true : activationLayer.isActive(x, y);
}
@Override
public double getMinX() {
return getWest();
}
@Override
public double getMaxX() {
return getEast();
}
@Override
public double getMinY() {
return getSouth();
}
@Override
public double getMaxY() {
return getNorth();
}
@Override
public Iterator iterator() {
return new Iterator| () {
int n = 0;
@Override
public boolean hasNext() {
return n < getCellCount();
}
@Override
public Cell next() {
return getCell(n++);
}
@Override
public void remove() {
}
};
}
/**
* Return a new grid with compatible geometry to this but limited to the
* bounding box of the passed shape, which must be included in our envelope.
* Also return an array with the x and y offsets of the grid within the
* original.
*
* @param shape
* @return
* @throws KlabException
*/
public IGrid cutToShape(ShapeValue shape) throws KlabException {
Envelope genv = new Envelope(getMinX(), getMaxX(), getMinY(), getMaxY());
Envelope senv = shape.getGeometry().getEnvelope().getEnvelopeInternal();
if (!genv.covers(senv)) {
return null;
}
/*
* adjusts envelope boundaries to cover original cells exactly
*/
double gxmin = senv.getMinX();
double gxmax = senv.getMaxX();
double dx = gxmax - gxmin;
double gymin = senv.getMinY();
double gymax = senv.getMaxY();
double dy = gymax - gymin;
int nx = (int) (dx / getCellWidth());
int ny = (int) (dy / getCellHeight());
if ((nx * getCellWidth()) < dx) {
nx++;
gxmin -= (getCellWidth() / 2);
gxmax += (getCellWidth() / 2);
}
if ((ny * getCellHeight()) < dy) {
ny++;
gymin -= (getCellHeight() / 2);
gymax += (getCellHeight() / 2);
}
int xofs = (int) ((gxmin - getMinX()) / getCellWidth());
int yofs = (int) ((gymin - getMinY()) / getCellHeight());
Grid ret = new Grid(shape, nx, ny);
ret.offsetInSupergrid = new int[] {xofs, yofs};
ret.superGridId = getSignature();
ret.createActivationLayer(shape);
return ret;
}
public String getParentGridId() {
return superGridId;
}
public int[] getOffsetInParentGrid() {
return offsetInSupergrid;
}
}
| | | |