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;
    }

}