CSIP_Georeferencing.java [src/java/methods] Revision: default Date:
/*
* To change this template, choose Tools | Templates
* and open the template in the editor.
*/
package methods;
import com.vividsolutions.jts.geom.Geometry;
import com.vividsolutions.jts.geom.MultiPolygon;
import com.vividsolutions.jts.geom.Polygon;
import com.vividsolutions.jts.io.ParseException;
import csip.SessionLogger;
import java.io.File;
import java.io.IOException;
import java.net.MalformedURLException;
import java.net.URL;
import java.util.ArrayList;
import java.util.List;
import java.util.concurrent.Callable;
import java.util.concurrent.ExecutionException;
import java.util.concurrent.ExecutorService;
import java.util.concurrent.Executors;
import java.util.concurrent.FutureTask;
import lamps.CSIP_Const;
import lamps.utils.Geospatial;
import lamps.utils.Partition;
import methods.objects.ERU;
import org.geotools.data.simple.SimpleFeatureCollection;
import org.geotools.data.simple.SimpleFeatureIterator;
import org.geotools.feature.FeatureCollection;
import org.geotools.geometry.jts.JTS;
import org.geotools.geometry.jts.ReferencedEnvelope;
import org.geotools.resources.coverage.IntersectUtils;
import org.opengis.feature.simple.SimpleFeature;
/**
*
* @author hokipka
*/
public class CSIP_Georeferencing {
/**
*
* @param dataDir
* @param env
* @param AOI
* @param multigeo
* @param LOG
* @return
* @throws MalformedURLException
* @throws IOException
* @throws InterruptedException
* @throws ExecutionException
*/
public static List<ERU> intersectStatesCMZs(File dataDir, ReferencedEnvelope env, List<Geometry> AOI, Geometry multigeo, boolean cmzs_switch, SessionLogger LOG) throws MalformedURLException, IOException, InterruptedException, ExecutionException, ParseException {
List<ERU> local_hrus = new ArrayList<>();
ArrayList<String> states = new ArrayList<>();
ArrayList<String> CMZs = new ArrayList<>();
ArrayList<String> County = new ArrayList<>();
ArrayList<SimpleFeature> geo_list_states = new ArrayList<>();
ArrayList<SimpleFeature> geo_list_cmzs = new ArrayList<>();
ArrayList<SimpleFeature> geo_list_county = new ArrayList<>();
//SimpleFeatureType simp_type = featureIterator.next().getType();
File Statesfile = new File(dataDir, CSIP_Const.USStatesshp);
URL dir_url = Statesfile.toURI().toURL();
FeatureCollection featureCollection = (SimpleFeatureCollection) Geospatial.featureCollectionFromShapefile(dir_url);
SimpleFeatureIterator featureIterator = (SimpleFeatureIterator) featureCollection.features();
Geometry extent = JTS.toGeometry(env);
if (multigeo != null) {
extent = multigeo;
}
while (featureIterator.hasNext()) {
SimpleFeature feature = featureIterator.next();
//simp_type = feature.getType();
Geometry geometry = (Geometry) feature.getDefaultGeometry();
geometry = Geospatial.roundCoordinates(geometry.convexHull());
if (extent.within(geometry)) {
states.add(feature.getProperty("STATE").getValue().toString());
if (!geo_list_states.contains(feature)) {
geo_list_states.add(feature);
}
break;
} else if (extent.intersects(geometry)) {
Geometry intersection = geometry.intersection(extent);
if (intersection instanceof MultiPolygon) {
MultiPolygon mp = (MultiPolygon) intersection;
for (int i = 0; i < mp.getNumGeometries(); i++) {
com.vividsolutions.jts.geom.Polygon g = (com.vividsolutions.jts.geom.Polygon) mp.getGeometryN(i);
//Geometry gIntersection = IntersectUtils.intersection(g, JTS.toGeometry(env));
if (!states.contains(feature.getProperty("STATE").getValue().toString())) {
states.add(feature.getProperty("STATE").getValue().toString());
}
if (!geo_list_states.contains(feature)) {
geo_list_states.add(feature);
}
//all_States_parts.union(gIntersection);
}
} else if (intersection instanceof Polygon) {
com.vividsolutions.jts.geom.Polygon g = (com.vividsolutions.jts.geom.Polygon) intersection;
//Geometry gIntersection = IntersectUtils.intersection(g, JTS.toGeometry(env));
if (!states.contains(feature.getProperty("STATE").getValue().toString())) {
states.add(feature.getProperty("STATE").getValue().toString());
}
if (!geo_list_states.contains(feature)) {
geo_list_states.add(feature);
}
//all_States_parts.union(gIntersection);
}
}
}
//w_log.println("===> Done rough State check. " + states.size());
LOG.info("===> Done rough State check. " + states.size());
for (int i = 0; i < states.size(); i++) {
//w_log.println(i + " " + states.get(i));
LOG.info(i + " " + states.get(i));
}
//w_log.flush();
featureIterator.close();
LOG.info("===> CMZs Switch : " + cmzs_switch);
if (!cmzs_switch) {
File CMZfile = new File(dataDir, CSIP_Const.CMZshp);
dir_url = CMZfile.toURI().toURL();
featureCollection = (SimpleFeatureCollection) Geospatial.featureCollectionFromShapefile(dir_url);
featureIterator = (SimpleFeatureIterator) featureCollection.features();
//List<Geometry> all_CMZs = new ArrayList<>();
while (featureIterator.hasNext()) {
SimpleFeature feature = featureIterator.next();
//simp_type = feature.getType();
Geometry geometry = (Geometry) feature.getDefaultGeometry();
geometry = Geospatial.roundCoordinates(geometry.convexHull());
if (extent.within(geometry)) {
String no_start_zero = feature.getProperty("CMZ").getValue().toString().replace("CMZ ", "");
if (feature.getProperty("CMZ").getValue().toString().replace("CMZ ", "").startsWith("0")) {
no_start_zero = feature.getProperty("CMZ").getValue().toString().replace("CMZ ", "").replace("0", "");
}
CMZs.add(no_start_zero);
if (!geo_list_cmzs.contains(feature)) {
geo_list_cmzs.add(feature);
}
break;
} else if (extent.intersects(geometry)) {
Geometry intersection = geometry.intersection(extent);
if (intersection instanceof MultiPolygon) {
// MultiPolygon mp = (MultiPolygon) intersection;
// double intersection_area_portion = 0;
// for (int i = 0; i < mp.getNumGeometries(); i++) {
// com.vividsolutions.jts.geom.Polygon g = (com.vividsolutions.jts.geom.Polygon) mp.getGeometryN(i);
// Geometry gIntersection = IntersectUtils.intersection(g, JTS.toGeometry(env));
// intersection_area_portion = (gIntersection.getArea() / AOI.get(i).getArea()) * 100;
// }
String no_start_zero = feature.getProperty("CMZ").getValue().toString().replace("CMZ ", "");
if (feature.getProperty("CMZ").getValue().toString().replace("CMZ ", "").startsWith("0")) {
no_start_zero = feature.getProperty("CMZ").getValue().toString().replace("CMZ ", "").replace("0", "");
}
if (!CMZs.contains(no_start_zero)) {
CMZs.add(no_start_zero);
}
//all_CMZs.add(gIntersection);
if (!geo_list_cmzs.contains(feature)) {
geo_list_cmzs.add(feature);
}
} else if (intersection instanceof Polygon) {
//com.vividsolutions.jts.geom.Polygon g = (com.vividsolutions.jts.geom.Polygon) intersection;
//Geometry gIntersection = IntersectUtils.intersection(g, JTS.toGeometry(env));
String no_start_zero = feature.getProperty("CMZ").getValue().toString().replace("CMZ ", "");
if (feature.getProperty("CMZ").getValue().toString().replace("CMZ ", "").startsWith("0")) {
no_start_zero = feature.getProperty("CMZ").getValue().toString().replace("CMZ ", "").replace("0", "");
}
if (!CMZs.contains(no_start_zero)) {
CMZs.add(no_start_zero);
}
//all_CMZs.add(gIntersection);
if (!geo_list_cmzs.contains(feature)) {
geo_list_cmzs.add(feature);
}
}
}
}
featureIterator.close();
}
//w_log.println("===> Done rough CMZ check. " + CMZs.size());
LOG.info("===> Done rough CMZ check. " + CMZs.size());
for (int i = 0; i < CMZs.size(); i++) {
//w_log.println(i + " " + CMZs.get(i));
LOG.info(i + " " + CMZs.get(i));
}
//w_log.flush();
//--------------Counties
featureIterator.close();
File Countyfile = new File(dataDir, CSIP_Const.USCountiesshp);
dir_url = Countyfile.toURI().toURL();
featureCollection = (SimpleFeatureCollection) Geospatial.featureCollectionFromShapefile(dir_url);
featureIterator = (SimpleFeatureIterator) featureCollection.features();
while (featureIterator.hasNext()) {
SimpleFeature feature = featureIterator.next();
Geometry geometry = (Geometry) feature.getDefaultGeometry();
geometry = Geospatial.roundCoordinates(geometry.convexHull());
if (extent.within(geometry)) {
County.add(feature.getProperty("NAME").getValue().toString());
if (!geo_list_county.contains(feature)) {
geo_list_county.add(feature);
}
break;
} else if (extent.intersects(geometry)) {
Geometry intersection = geometry.intersection(extent);
if (intersection instanceof MultiPolygon) {
if (!County.contains(feature.getProperty("NAME").getValue().toString())) {
County.add(feature.getProperty("NAME").getValue().toString());
}
if (!geo_list_county.contains(feature)) {
geo_list_county.add(feature);
}
} else if (intersection instanceof Polygon) {
if (!County.contains(feature.getProperty("NAME").getValue().toString())) {
County.add(feature.getProperty("NAME").getValue().toString());
}
if (!geo_list_county.contains(feature)) {
geo_list_county.add(feature);
}
}
}
}
featureIterator.close();
LOG.info("===> Done rough County check. " + County.size());
for (int i = 0; i < County.size(); i++) {
LOG.info(i + " " + County.get(i));
}
if (CMZs.size() < 2 && states.size() < 2 && County.size() < 2) {
if (CMZs.size() < 1) {
LOG.info("===> No CMZ ");
CMZs.add("XX");
}
for (int i = 0; i < AOI.size(); i++) {
ERU new_hru = new ERU();
new_hru.ID = i;
new_hru.state = states;
new_hru.cmz = CMZs;
new_hru.county = County;
new_hru.multiCMZs = false;
new_hru.multiStates = false;
new_hru.geom = AOI.get(i);
new_hru.irrigated = false;
new_hru.irrigated_area = 0.0;
local_hrus.add(new_hru);
}
} else {
local_hrus = concurrentOverlay(AOI, geo_list_cmzs, geo_list_states, geo_list_county, LOG);
}
LOG.info("===> Done State & CMZ & County check for all " + local_hrus.size());
return local_hrus;
}
/**
*
* @param AOI
* @param list_cmzs
* @param list_states
* @param list_county
* @param LOG
* @return
* @throws InterruptedException
* @throws ExecutionException
*/
public static List<ERU> concurrentOverlay(List<Geometry> AOI, ArrayList<SimpleFeature> list_cmzs, ArrayList<SimpleFeature> list_states, ArrayList<SimpleFeature> list_county, SessionLogger LOG) throws InterruptedException, ExecutionException {
final int threadNum = Partition.getThreadCount(AOI.size());
final int first = 0;
final int last = AOI.size() - 1;
final List<Geometry> AOI_fin = AOI;
final ArrayList<SimpleFeature> fin_coll_states = list_states;
final ArrayList<SimpleFeature> fin_coll_cmz = list_cmzs;
final ArrayList<SimpleFeature> fin_coll_county = list_county;
final SessionLogger fSL = LOG;
final int p = (last + 1) / threadNum;
int s = (last + 1) - p * threadNum;
List<ERU> allHRUs = new ArrayList();
// Prepare to execute and store the Futures
ExecutorService executor = Executors.newFixedThreadPool(threadNum);
List<FutureTask<List<ERU>>> taskList = new ArrayList<>();
for (int j = 0; j < threadNum; j++) {
final int ii = j;
FutureTask<List<ERU>> futureTask_1 = new FutureTask<>(new Callable<List<ERU>>() {
int inlast = (ii == threadNum - 1 ? last : p + (ii * p));
int infirst = (ii == 0 ? first : (p * ii) + 1);
@Override
public List<ERU> call() throws Exception {
return CSIP_Georeferencing.overlay(infirst, inlast, AOI_fin, fin_coll_cmz, fin_coll_states, fin_coll_county, ii, fSL);
}
});
taskList.add(futureTask_1);
executor.execute(futureTask_1);
}
// Wait until all results are available and combine them at the same time
for (int j = 0; j < threadNum; j++) {
FutureTask<List<ERU>> futureTask = taskList.get(j);
for (int i = 0; i < futureTask.get().size(); i++) {
allHRUs.add(futureTask.get().get(i));
}
}
executor.shutdown();
return allHRUs;
}
/**
*
* @param first
* @param last
* @param AOI
* @param list_cmzs
* @param list_states
* @param list_county
* @param part
* @param LOG
* @return
*/
public static List<ERU> overlay(int first, int last, List<Geometry> AOI, ArrayList<SimpleFeature> list_cmzs, ArrayList<SimpleFeature> list_states, ArrayList<SimpleFeature> list_county, int part, SessionLogger LOG) throws ParseException {
int count = 0;
List<ERU> allHRUs = new ArrayList();
int part_all = last - first;
for (int i = first; i <= last; i++) {
//LOG.info("===> #AOI " + i);
ArrayList<String> CMZs = new ArrayList<>();
ArrayList<String> states = new ArrayList<>();
ArrayList<String> counties = new ArrayList<>();
double geometry_AOI_area_size = AOI.get(i).getArea();
double majority_size = 12;
if (geometry_AOI_area_size > 5000000) {
majority_size = 1;
}
count += 1;
if (count == 100) {
count = 0;
double progress = ((double) i - first) / ((double) part_all);
progress = (double) progress * 100.00;
LOG.info(part + ".-" + ((int) progress) + "%...");
}
//LOG.info("===> Going to real check CMZs");
if (list_cmzs.size() > 0) {
for (SimpleFeature feature : list_cmzs) {
Geometry geometry = (Geometry) feature.getDefaultGeometry();
geometry = Geospatial.roundCoordinates(geometry.convexHull());
if (AOI.get(i).within(geometry)) {
String no_start_zero = feature.getProperty("CMZ").getValue().toString().replace("CMZ ", "");
if (feature.getProperty("CMZ").getValue().toString().replace("CMZ ", "").startsWith("0")) {
no_start_zero = feature.getProperty("CMZ").getValue().toString().replace("CMZ ", "").replace("0", "");
}
CMZs.add(no_start_zero);
break;
} else if (AOI.get(i).intersects(geometry) || AOI.get(i).crosses(geometry) || AOI.get(i).overlaps(geometry) || AOI.get(i).touches(geometry)) {
Geometry intersection = AOI.get(i).intersection(geometry);
if (intersection instanceof MultiPolygon) {
MultiPolygon mp = (MultiPolygon) intersection;
double sum_area_part = 0;
for (int uu = 0; uu < mp.getNumGeometries(); uu++) {
com.vividsolutions.jts.geom.Polygon g = (com.vividsolutions.jts.geom.Polygon) mp.getGeometryN(uu);
Geometry gIntersection = IntersectUtils.intersection(g, AOI.get(i));
sum_area_part = sum_area_part + ((gIntersection.getArea() / AOI.get(i).getArea()) * 100);
}
String no_start_zero = feature.getProperty("CMZ").getValue().toString().replace("CMZ ", "");
if (feature.getProperty("CMZ").getValue().toString().replace("CMZ ", "").startsWith("0")) {
no_start_zero = feature.getProperty("CMZ").getValue().toString().replace("CMZ ", "").replace("0", "");
}
if (!CMZs.contains(no_start_zero) && sum_area_part > majority_size) {
CMZs.add(no_start_zero);
}
}
if (intersection instanceof Polygon) {
com.vividsolutions.jts.geom.Polygon g = (com.vividsolutions.jts.geom.Polygon) intersection;
Geometry gIntersection = IntersectUtils.intersection(g, AOI.get(i));
double intersection_area_portion = (gIntersection.getArea() / AOI.get(i).getArea()) * 100;
String no_start_zero = feature.getProperty("CMZ").getValue().toString().replace("CMZ ", "");
if (feature.getProperty("CMZ").getValue().toString().replace("CMZ ", "").startsWith("0")) {
no_start_zero = feature.getProperty("CMZ").getValue().toString().replace("CMZ ", "").replace("0", "");
}
if (!CMZs.contains(no_start_zero) && intersection_area_portion > majority_size) {
CMZs.add(no_start_zero);
}
}
}
}
}
//LOG.info("===> Done real check CMZs");
//LOG.info("===> Going to real check States");
states = Geospatial.overlayResultList(AOI.get(i), list_states, "STATE", majority_size, LOG);
//LOG.info("===> Done real check States");
//LOG.info("===> Going to real check Counties");
counties = Geospatial.overlayResultList(AOI.get(i), list_county, "NAME", majority_size, LOG);
//LOG.info("===> Done real check Counties");
ERU new_hru = new ERU();
new_hru.ID = i;
new_hru.state = states;
new_hru.cmz = CMZs;
new_hru.county = counties;
if (new_hru.county.isEmpty()) {
for (SimpleFeature feature : list_county) {
Geometry geometry = (Geometry) feature.getDefaultGeometry();
if (AOI.get(i).getCentroid().within(geometry.getEnvelope())) {
counties.add(feature.getProperty("NAME").getValue().toString().trim());
break;
}
}
new_hru.county = counties;
if (new_hru.county.isEmpty()) {
LOG.info("===> No County " + new_hru.ID);
}
}
if (new_hru.state.isEmpty()) {
for (SimpleFeature feature : list_states) {
Geometry geometry = (Geometry) feature.getDefaultGeometry();
if (AOI.get(i).getCentroid().within(geometry.getEnvelope())) {
states.add(feature.getProperty("STATE").getValue().toString().trim());
break;
} else if (AOI.get(i).intersects(geometry.getEnvelope()) || AOI.get(i).crosses(geometry.getEnvelope()) || AOI.get(i).overlaps(geometry.getEnvelope()) || AOI.get(i).touches(geometry.getEnvelope())) {
Geometry intersection = AOI.get(i).intersection(geometry.getEnvelope());
if (intersection instanceof MultiPolygon) {
MultiPolygon mp = (MultiPolygon) intersection;
double sum_area_part = 0;
for (int uu = 0; uu < mp.getNumGeometries(); uu++) {
com.vividsolutions.jts.geom.Polygon g = (com.vividsolutions.jts.geom.Polygon) mp.getGeometryN(uu);
Geometry gIntersection = IntersectUtils.intersection(g, AOI.get(i));
sum_area_part = sum_area_part + ((gIntersection.getArea() / AOI.get(i).getArea()) * 100);
}
if (!states.contains(feature.getProperty("STATE").getValue().toString())) {
states.add(feature.getProperty("STATE").getValue().toString());
}
}
if (intersection instanceof Polygon) {
com.vividsolutions.jts.geom.Polygon g = (com.vividsolutions.jts.geom.Polygon) intersection;
Geometry gIntersection = IntersectUtils.intersection(g, AOI.get(i));
//double sum_area_part = 0;
//sum_area_part = sum_area_part + ((gIntersection.getArea() / AOI.get(i).getArea()) * 100);
if (!states.contains(feature.getProperty("STATE").getValue().toString())) {
states.add(feature.getProperty("STATE").getValue().toString());
}
}
}
}
new_hru.state = states;
if (new_hru.state.isEmpty()) {
LOG.info("===> No State " + new_hru.ID);
}
}
if (new_hru.cmz.isEmpty()) {
//LOG.info("===> in more overlay check CMZ is empty 1. !");
for (SimpleFeature feature : list_cmzs) {
Geometry geometry = (Geometry) feature.getDefaultGeometry();
if (AOI.get(i).getCentroid().within(geometry.getEnvelope())) {
String no_start_zero = feature.getProperty("CMZ").getValue().toString().replace("CMZ ", "");
if (feature.getProperty("CMZ").getValue().toString().replace("CMZ ", "").startsWith("0")) {
no_start_zero = feature.getProperty("CMZ").getValue().toString().replace("CMZ ", "").replace("0", "");
}
CMZs.add(no_start_zero);
break;
}
}
new_hru.cmz = CMZs;
if (new_hru.cmz.isEmpty()) {
//LOG.info("===> in more overlay check CMZ is empty 2. !");
LOG.info("===> No CMZ " + new_hru.ID);
CMZs.add("XX");
}
}
new_hru.state = states;
new_hru.cmz = CMZs;
new_hru.county = counties;
new_hru.multiCMZs = false;
new_hru.multiStates = false;
new_hru.geom = AOI.get(i);
new_hru.irrigated = false;
new_hru.irrigated_area = 0.0;
if (states.size() > 1) {
new_hru.multiStates = true;
}
if (CMZs.size() > 1) {
new_hru.multiCMZs = true;
}
allHRUs.add(new_hru);
}
return allHRUs;
}
}