GeometryInput.java [src/java/lamps/io] Revision: default Date:
/*
* To change this license header, choose License Headers in Project Properties.
* To change this template file, choose Tools | Templates
* and open the template in the editor.
*/
package lamps.io;
import com.vividsolutions.jts.geom.Geometry;
import com.vividsolutions.jts.geom.GeometryFactory;
import com.vividsolutions.jts.geom.LinearRing;
import com.vividsolutions.jts.geom.MultiPolygon;
import com.vividsolutions.jts.geom.Polygon;
import com.vividsolutions.jts.io.WKTReader;
import csip.SessionLogger;
import java.io.BufferedInputStream;
import java.io.File;
import java.io.FileInputStream;
import java.io.FileOutputStream;
import java.io.InputStream;
import java.io.OutputStream;
import java.net.URL;
import java.util.ArrayList;
import java.util.Collection;
import java.util.Iterator;
import java.util.List;
import java.util.zip.ZipEntry;
import java.util.zip.ZipInputStream;
import javax.xml.parsers.DocumentBuilder;
import javax.xml.parsers.DocumentBuilderFactory;
import lamps.CSIP_Const;
import lamps.utils.Geospatial;
import lamps.utils.Partition;
import methods.CSIP_Georeferencing;
import methods.objects.ERU;
import org.geotools.data.FeatureSource;
import org.geotools.data.FileDataStore;
import org.geotools.data.FileDataStoreFinder;
import org.geotools.feature.FeatureCollection;
import org.geotools.feature.FeatureIterator;
import org.geotools.geojson.feature.FeatureJSON;
import org.geotools.geometry.jts.JTS;
import org.geotools.geometry.jts.JTSFactoryFinder;
import org.geotools.geometry.jts.ReferencedEnvelope;
import org.geotools.kml.v22.KMLConfiguration;
import org.geotools.referencing.CRS;
import org.geotools.xml.Parser;
import org.opengis.feature.simple.SimpleFeature;
import org.opengis.feature.type.FeatureType;
import org.opengis.referencing.crs.CRSAuthorityFactory;
import org.opengis.referencing.crs.CoordinateReferenceSystem;
import org.opengis.referencing.operation.MathTransform;
import org.w3c.dom.Document;
import org.w3c.dom.Element;
import org.w3c.dom.Node;
import org.w3c.dom.NodeList;
/**
*
* @author od
*/
public class GeometryInput {
/**
*
* @param dataDir
* @param hru_url
* @param Input_Geometries
* @param outputDir
* @param unique_states
* @param unique_cmzs
* @param envelope
* @param LOG
* @return
* @throws Exception
*/
public static List<ERU> read(File dataDir, URL hru_url, List<Geometry> Input_Geometries, File outputDir, List<String> unique_states, List<String> unique_cmzs, List<Geometry> envelope, boolean cmzs_switch, SessionLogger LOG) throws Exception {
//NASS CropScape max. extend coordinates
double min_x = 28000000;
double max_x = -28000000;
double min_y = 3000000;
double max_y = 280000;
List<Integer> LAMPS_IDs = new ArrayList<>();
List<ERU> local_hrus = new ArrayList<>();
Geometry multigeo = null;
CoordinateReferenceSystem srcCRS = null;
ReferencedEnvelope env = null;
GeometryFactory geometryFactory = JTSFactoryFinder.getGeometryFactory(null);
//shape file detection
if (hru_url.toURI().toString().endsWith(".shp")) {
//w_log.println(" .shp File");
LOG.info(" .shp File");
//w_log.flush();
//target coordinate system in Albers Equal Area projection
CoordinateReferenceSystem targetCRS = CRS.parseWKT(CSIP_Const.WKT_CooSyst);
File dir = new File(hru_url.toURI());
//shapefile reading and original coodinates system detection
FileDataStore dataStore = FileDataStoreFinder.getDataStore(dir);
FeatureSource shapefileSource = dataStore.getFeatureSource();
FeatureType schema = shapefileSource.getSchema();
srcCRS = schema.getCoordinateReferenceSystem();
dataStore.dispose();
MathTransform transform = CRS.findMathTransform(srcCRS, targetCRS, true);
//shapefile into a featureCollection
FeatureCollection featureCollection = (FeatureCollection) Geospatial.featureCollectionFromShapefile(hru_url);
int number_polys = featureCollection.size();
FeatureIterator iterator = featureCollection.features();
//loop over each feature/geometry/polygon to get the extend coordinates
int o = 0;
while (iterator.hasNext() && o < number_polys) {
SimpleFeature feature = (SimpleFeature) iterator.next();
Geometry geometry = (Geometry) feature.getDefaultGeometry();
Geometry targetGeometry = JTS.transform(geometry, transform);
if (!targetGeometry.isValid()) {
targetGeometry = targetGeometry.convexHull();
}
int last = targetGeometry.getCoordinates().length - 1;
if (!targetGeometry.getCoordinates()[0].equals(targetGeometry.getCoordinates()[last])) {
targetGeometry = targetGeometry.convexHull();
}
if (targetGeometry.getBoundary().getEnvelopeInternal().getMaxX() > max_x) {
max_x = targetGeometry.getEnvelopeInternal().getMaxX();
}
if (targetGeometry.getBoundary().getEnvelopeInternal().getMinX() < min_x) {
min_x = targetGeometry.getEnvelopeInternal().getMinX();
}
if (targetGeometry.getBoundary().getEnvelopeInternal().getMaxY() > max_y) {
max_y = targetGeometry.getEnvelopeInternal().getMaxY();
}
if (targetGeometry.getBoundary().getEnvelopeInternal().getMinY() < min_y) {
min_y = targetGeometry.getEnvelopeInternal().getMinY();
}
if (!Input_Geometries.contains(targetGeometry)) {
Input_Geometries.add(targetGeometry);
LAMPS_IDs.add(o);
o++;
}
}
iterator.close();
Polygon[] all_polys = new Polygon[Input_Geometries.size()];
for (int i = 0; i < Input_Geometries.size(); i++) {
Geometry targetGeometry = Input_Geometries.get(i);
LinearRing lr = geometryFactory.createLinearRing(targetGeometry.convexHull().getCoordinates());
all_polys[i] = geometryFactory.createPolygon(lr, null);
}
MultiPolygon multpoly = geometryFactory.createMultiPolygon(all_polys);
multigeo = multpoly.convexHull();
env = new ReferencedEnvelope(min_x, max_x, min_y, max_y, targetCRS);
}
if (hru_url.toURI().toString().endsWith(".kmz")) {
//w_log.println(" .kmz File");
LOG.info(" .kmz File");
//w_log.flush();
File ff = new File(outputDir, "kmztokml.kml");
OutputStream out = new FileOutputStream(ff);
//System.out.println(" " + hru_url.getPath());
FileInputStream fin = new FileInputStream(hru_url.getPath());
BufferedInputStream bin = new BufferedInputStream(fin);
ZipInputStream zin = new ZipInputStream(bin);
ZipEntry ze = null;
while ((ze = zin.getNextEntry()) != null) {
//if (ze.getName().equals("doc.kml")) {
byte[] buffer = new byte[8192];
int len;
while ((len = zin.read(buffer)) != -1) {
out.write(buffer, 0, len);
}
out.close();
break;
//}
}
hru_url = ff.toURI().toURL();
fin.close();
bin.close();
zin.close();
}
if (hru_url.toURI().toString().endsWith(".kml")) {
LOG.info(" .kml File");
WKTReader reader = new WKTReader();
srcCRS = CRS.parseWKT(CSIP_Const.WKT_KML);
CoordinateReferenceSystem targetCRS = CRS.parseWKT(CSIP_Const.WKT_CooSyst);
MathTransform transform = CRS.findMathTransform(srcCRS, targetCRS, true);
InputStream isn = new FileInputStream(new File(hru_url.toURI()));
boolean success = false;
try {
Parser parser = new Parser(new KMLConfiguration());
SimpleFeature f = (SimpleFeature) parser.parse(isn);
Collection placemarks = (Collection) f.getAttribute("Feature");
Iterator place_iter = placemarks.iterator();
while (place_iter.hasNext()) {
SimpleFeature feature = (SimpleFeature) place_iter.next();
Geometry geometry = (Geometry) feature.getDefaultGeometry();
Geometry targetGeometry = JTS.transform(geometry, transform);
}
success = true;
throw new Exception();
} catch (Exception e) {
} finally {
if (!success) {
LOG.info(" old kml read method 1.0 ");
} else {
LOG.info(" new kml read method 2.2 ");
}
}
isn = new FileInputStream(new File(hru_url.toURI()));
if (!success) {
isn = new FileInputStream(new File(hru_url.toURI()));
DocumentBuilderFactory dbFactory = DocumentBuilderFactory.newInstance();
DocumentBuilder dBuilder = dbFactory.newDocumentBuilder();
Document doc = dBuilder.parse(isn);
doc.getDocumentElement().normalize();
NodeList PlacemarkList = doc.getElementsByTagName("Placemark");
for (int i = 0; i < PlacemarkList.getLength(); i++) {
Node n = PlacemarkList.item(i);
if (n.getNodeType() == Node.ELEMENT_NODE) {
Element eElement = (Element) n.getChildNodes();
NodeList Place_coords = eElement.getElementsByTagName("coordinates");
if (Place_coords.getLength() > 1) {
Geometry all = null;
for (int j = 0; j < Place_coords.getLength(); j++) {
String[] coords = Place_coords.item(j).getTextContent().trim().split(" ");
//String[] coords = nNode.getFirstChild().getTextContent().trim().split(" ");
String comma_changed = "";
for (int ai = 0; ai < coords.length; ai++) {
String[] each_coords = coords[ai].trim().split(",");
comma_changed = comma_changed + each_coords[0] + " " + each_coords[1] + ",";
if (coords.length - 1 == ai) {
comma_changed = comma_changed + each_coords[0] + " " + each_coords[1];
}
}
Geometry tShape = reader.read("POLYGON((" + comma_changed + "))");
Geometry targetGeometry = JTS.transform(tShape, transform);
if (all == null) {
all = targetGeometry;
} else {
all = all.union(targetGeometry);
}
if (!all.isValid()) {
all = all.convexHull();
}
}
int last = all.getCoordinates().length - 1;
if (!all.getCoordinates()[0].equals(all.getCoordinates()[last])) {
all = all.convexHull();
}
if (all.getBoundary().getEnvelopeInternal().getMaxX() > max_x) {
max_x = all.getEnvelopeInternal().getMaxX();
}
if (all.getBoundary().getEnvelopeInternal().getMinX() < min_x) {
min_x = all.getEnvelopeInternal().getMinX();
}
if (all.getBoundary().getEnvelopeInternal().getMaxY() > max_y) {
max_y = all.getEnvelopeInternal().getMaxY();
}
if (all.getBoundary().getEnvelopeInternal().getMinY() < min_y) {
min_y = all.getEnvelopeInternal().getMinY();
}
Input_Geometries.add(all);
LAMPS_IDs.add(i);
} else {
//ok placemark
String[] coords = Place_coords.item(0).getTextContent().trim().split(" ");
String comma_changed = "";
for (int ci = 0; ci < coords.length; ci++) {
String[] each_coords = coords[ci].trim().split(",");
comma_changed = comma_changed + each_coords[0] + " " + each_coords[1] + ",";
if (coords.length - 1 == ci) {
comma_changed = comma_changed + each_coords[0] + " " + each_coords[1];
}
}
Geometry tShape = reader.read("POLYGON((" + comma_changed + "))");
Geometry targetGeometry = JTS.transform(tShape, transform);
if (!targetGeometry.isValid()) {
targetGeometry = targetGeometry.convexHull();
}
int last = targetGeometry.getCoordinates().length - 1;
if (!targetGeometry.getCoordinates()[0].equals(targetGeometry.getCoordinates()[last])) {
targetGeometry = targetGeometry.convexHull();
}
if (targetGeometry.getBoundary().getEnvelopeInternal().getMaxX() > max_x) {
max_x = targetGeometry.getEnvelopeInternal().getMaxX();
}
if (targetGeometry.getBoundary().getEnvelopeInternal().getMinX() < min_x) {
min_x = targetGeometry.getEnvelopeInternal().getMinX();
}
if (targetGeometry.getBoundary().getEnvelopeInternal().getMaxY() > max_y) {
max_y = targetGeometry.getEnvelopeInternal().getMaxY();
}
if (targetGeometry.getBoundary().getEnvelopeInternal().getMinY() < min_y) {
min_y = targetGeometry.getEnvelopeInternal().getMinY();
}
Input_Geometries.add(targetGeometry);
LAMPS_IDs.add(i);
}
}
}
} else {
isn = new FileInputStream(new File(hru_url.toURI()));
Parser parser = new Parser(new KMLConfiguration());
SimpleFeature f = (SimpleFeature) parser.parse(isn);
Collection placemarks = (Collection) f.getAttribute("Feature");
Iterator place_iter = placemarks.iterator();
while (place_iter.hasNext()) {
SimpleFeature feature = (SimpleFeature) place_iter.next();
Geometry geometry = (Geometry) feature.getDefaultGeometry();
Geometry targetGeometry = JTS.transform(geometry, transform);
targetGeometry = Geospatial.roundCoordinates(targetGeometry);
if (!targetGeometry.isValid()) {
targetGeometry = targetGeometry.convexHull();
}
int last = targetGeometry.getCoordinates().length - 1;
if (!targetGeometry.getCoordinates()[0].equals(targetGeometry.getCoordinates()[last])) {
targetGeometry = targetGeometry.convexHull();
}
if (targetGeometry.getBoundary().getEnvelopeInternal().getMaxX() > max_x) {
max_x = targetGeometry.getEnvelopeInternal().getMaxX();
}
if (targetGeometry.getBoundary().getEnvelopeInternal().getMinX() < min_x) {
min_x = targetGeometry.getEnvelopeInternal().getMinX();
}
if (targetGeometry.getBoundary().getEnvelopeInternal().getMaxY() > max_y) {
max_y = targetGeometry.getEnvelopeInternal().getMaxY();
}
if (targetGeometry.getBoundary().getEnvelopeInternal().getMinY() < min_y) {
min_y = targetGeometry.getEnvelopeInternal().getMinY();
}
Input_Geometries.add(targetGeometry);
//LOG.info(" " + feature.getID()); // ------------------------------- GEtting ID for HRU-ID
LAMPS_IDs.add(Integer.parseInt(feature.getID()));
}
}
env = new ReferencedEnvelope(min_x, max_x, min_y, max_y, targetCRS);
isn.close();
}
if (hru_url.toURI().toString().endsWith(".geojson")) {
LOG.info(" .geojson File");
InputStream isn = new FileInputStream(new File(hru_url.toURI()));
CRSAuthorityFactory factory = CRS.getAuthorityFactory(true);
srcCRS = factory.createCoordinateReferenceSystem("EPSG:4326");
CoordinateReferenceSystem targetCRS = CRS.parseWKT(CSIP_Const.WKT_CooSyst);
MathTransform transform = CRS.findMathTransform(srcCRS, targetCRS, true);
FeatureJSON fj = new FeatureJSON();
FeatureIterator<SimpleFeature> features = fj.streamFeatureCollection(isn);
int x = 0;
while (features.hasNext()) {
SimpleFeature feature = (SimpleFeature) features.next();
Geometry geometry = (Geometry) feature.getDefaultGeometry();
Geometry targetGeometry = JTS.transform(geometry, transform);
if (!targetGeometry.isValid()) {
targetGeometry = targetGeometry.convexHull();
}
int last = targetGeometry.getCoordinates().length - 1;
if (!targetGeometry.getCoordinates()[0].equals(targetGeometry.getCoordinates()[last])) {
targetGeometry = targetGeometry.convexHull();
}
if (targetGeometry.getBoundary().getEnvelopeInternal().getMaxX() > max_x) {
max_x = targetGeometry.getEnvelopeInternal().getMaxX();
}
if (targetGeometry.getBoundary().getEnvelopeInternal().getMinX() < min_x) {
min_x = targetGeometry.getEnvelopeInternal().getMinX();
}
if (targetGeometry.getBoundary().getEnvelopeInternal().getMaxY() > max_y) {
max_y = targetGeometry.getEnvelopeInternal().getMaxY();
}
if (targetGeometry.getBoundary().getEnvelopeInternal().getMinY() < min_y) {
min_y = targetGeometry.getEnvelopeInternal().getMinY();
}
if (!Input_Geometries.contains(targetGeometry)) {
Input_Geometries.add(targetGeometry);
//LOG.info(" " + feature.getID()); // ------------------------------- GEtting ID for HRU-ID
LAMPS_IDs.add(x);
x++;
}
}
features.close();
env = new ReferencedEnvelope(min_x, max_x, min_y, max_y, targetCRS);
isn.close();
}
//start here for the sequential run based on 10000 HRUs
List<List<Geometry>> tenk_geo_lists = new ArrayList<>();
int hrus_pre_count = Input_Geometries.size();
if (hrus_pre_count > 2000) {
int end = Input_Geometries.size() - 1;
int ten_times = (end + 1) / 1000;
int rest = (end + 1) - ten_times * 1000;
if (rest > 0) {
ten_times++;
}
LOG.info("===> Georeferencing input times : " + ten_times);
tenk_geo_lists = Partition.chopIntoParts(Input_Geometries, ten_times);
} else {
tenk_geo_lists.add(new ArrayList<>(Input_Geometries));
}
LOG.info("===> number of list parts : " + tenk_geo_lists.size());
int number_of_chopParts = 0;
List<ERU> local_hrus_pre = new ArrayList<>();
for (int numberLists = 0; numberLists < tenk_geo_lists.size(); numberLists++) {
number_of_chopParts += tenk_geo_lists.get(numberLists).size();
LOG.info(" ** ");
LOG.info("===> total number of geometries in loop : " + number_of_chopParts);
LOG.info(" ** ");
//calling the georeferencing method to get a states and CMZ linkage for each polygon
local_hrus_pre = CSIP_Georeferencing.intersectStatesCMZs(dataDir, env, tenk_geo_lists.get(numberLists), multigeo, cmzs_switch, LOG);
for (ERU hruinlist : local_hrus_pre) {
local_hrus.add(hruinlist);
}
}
//calling the georeferencing method to get a states and CMZ linkage for each polygon
//local_hrus = CSIP_Georeferencing.intersectStatesCMZs(dataDir, env, Input_Geometries, multigeo, LOG);
//calling the georeferencing method to get a states and CMZ linkage for each polygon
Geometry extent = JTS.toGeometry(env);
if (multigeo != null) {
extent = multigeo;
}
envelope.add(extent);
Input_Geometries.clear();
for (int i = 0; i < local_hrus.size(); i++) {
local_hrus.get(i).ID = LAMPS_IDs.get(i);
Input_Geometries.add(local_hrus.get(i).geom);
ArrayList<String> cmzs = local_hrus.get(i).cmz;
for (int j = 0; j < cmzs.size(); j++) {
if (!unique_cmzs.contains(cmzs.get(j))) {
unique_cmzs.add(cmzs.get(j));
}
}
ArrayList<String> states = local_hrus.get(i).state;
for (int j = 0; j < states.size(); j++) {
if (!unique_states.contains(states.get(j))) {
unique_states.add(states.get(j));
}
}
}
LOG.info("number of entities: " + local_hrus.size() + " in: " + unique_cmzs.size() + " CMZ(s) " + unique_cmzs + " (" + unique_states.size() + " State(s) " + unique_states + ")");
return local_hrus;
}
}