gCubeSystem
/
ecological-engine-geospatial-extensions
15
0
Fork 0
Code Issues Pull Requests Releases Wiki Activity

MaxEnt implementation finished and almost tested 

git-svn-id: https://svn.d4science.research-infrastructures.eu/gcube/trunk/data-analysis/EcologicalEngineGeoSpatialExtension@100360 82a268e6-3cf1-43bd-a215-b396298e98cf
This commit is contained in:
Gianpaolo Coro 2014-10-02 12:24:14 +00:00
parent b6c2356e6d
commit 96fd4efcb0
8 changed files with 862 additions and 40 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

301
src/main/java/density/Maxent.java Normal file
View File

@ -0,0 +1,301 @@
package density;
import java.io.BufferedReader;
import java.io.File;
import java.io.FileReader;
import java.util.ArrayList;
import java.util.LinkedHashMap;
import java.util.List;
import java.util.Map;
import org.gcube.dataanalysis.ecoengine.utils.Transformations;
public class Maxent {
public String samplesFilePath;
public String environmentalLayersLocation;
public String outputDirectory;
public int maxIterations;
public double defaultPrevalence;
public int noDataValue;
public List<String> envVariables = new ArrayList<String>();
public static Params getDefaultParameters() {
Params p = new Params();
p.setResponsecurves(true);
p.setPictures(true);
p.setJackknife(false);
p.setOutputformat("Logistic");
p.setOutputfiletype("asc");
p.setRandomseed(false);
p.setLogscale(true);
p.setWarnings(true);
p.setTooltips(false);
p.setAskoverwrite(false);
p.setSkipifexists(false);
p.setRemoveduplicates(true);
p.setWriteclampgrid(false);
p.setWritemess(false);
p.setRandomtestpoints(0);
p.setBetamultiplier(1);
p.setMaximumbackground(10000);
p.setReplicates(1);
p.setReplicatetype("Crossvalidate");
p.setPerspeciesresults(false);
p.setWritebackgroundpredictions(false);
p.setBiasisbayesianprior(false);
p.setResponsecurvesexponent(false);
p.setLinear(true);
p.setQuadratic(true);
p.setProduct(true);
p.setThreshold(true);
p.setHinge(true);
p.setPolyhedral(true);
p.setAddsamplestobackground(true);
p.setAddallsamplestobackground(false);
p.setAutorun(false);
p.setAutofeature(true);
p.setDosqrtcat(false);
p.setWriteplotdata(false);
p.setFadebyclamping(false);
p.setExtrapolate(true);
p.setVisible(false);
p.setGivemaxaucestimate(true);
p.setDoclamp(true);
p.setOutputgrids(true);
p.setPlots(true);
p.setAppendtoresultsfile(false);
p.setParallelupdatefrequency(30);
p.setMaximumiterations(1000);
p.setConvergencethreshold(0.00001);
p.setAdjustsampleradius(0);
p.setThreads(2);
p.setLq2lqptthreshold(80);
p.setL2lqthreshold(10);
p.setHingethreshold(15);
p.setBeta_threshold(-1);
p.setBeta_categorical(-1);
p.setBeta_lqp(-1);
p.setBeta_hinge(-1);
p.setBiastype(0);
p.setLogfile("maxent.log");
p.setScientificpattern("#.#####E0");
p.setCache(false);
p.setCachefeatures(false);
p.setDefaultprevalence(0.5);
p.setVerbose(true);
p.setAllowpartialdata(false);
p.setPrefixes(true);
p.setPrintversion(false);
p.setNodata(-9999);
p.setNceas(false);
p.setMinclamping(false);
p.setManualreplicates(false);
p.setSamplesfile("CarcharodonPoints.csv");
p.setEnvironmentallayers("./");
p.setOutputdirectory("./maxentout/");
return p;
}
public void executeMaxEnt() {
File outDir = new File(outputDirectory);
if (!outDir.exists())
outDir.mkdir();
Params p = getDefaultParameters();
p.setSamplesfile(samplesFilePath);
p.setEnvironmentallayers(environmentalLayersLocation);
p.setOutputdirectory(outputDirectory);
p.setMaximumiterations(maxIterations);
p.setDefaultprevalence(defaultPrevalence);
p.setNodata(noDataValue);
Utils.applyStaticParams(p);
p.setSelections();
Runner runner = new Runner(p);
runner.start();
runner.end();
}
public Maxent(String samplesFilePath, String environmentalLayersLocation, String outputDirectory, int maxIterations, double defaultPrevalence, int noDataValue) {
this.samplesFilePath = samplesFilePath;
this.environmentalLayersLocation = environmentalLayersLocation;
this.outputDirectory = outputDirectory;
this.maxIterations = maxIterations;
this.defaultPrevalence = defaultPrevalence;
this.noDataValue = noDataValue;
File layersLocation = new File(environmentalLayersLocation);
File[] list = layersLocation.listFiles();
for (File f : list) {
if (f.getName().endsWith(".asc")) {
envVariables.add(f.getName());
}
}
}
public String getSpeciesName() throws Exception{
Map<String, String> value = getOutputValues("Species");
String species = value.get("Species");
return species;
}
public String getResult() throws Exception{
String species = getSpeciesName();
File f = new File(outputDirectory, species + ".asc");
if (f.exists()) {
return f.getAbsolutePath();
}
return null;
}
public String getWorldPlot() throws Exception{
File f = new File(outputDirectory, "plots/"+getSpeciesName()+".png");
if (f.exists()) {
return f.getAbsolutePath();
}
return null;
}
public String getOmissionPlot() throws Exception{
File f = new File(outputDirectory, "plots/"+getSpeciesName()+"_omission.png");
if (f.exists()) {
return f.getAbsolutePath();
}
return null;
}
public String getROCPlot() throws Exception{
File f = new File(outputDirectory, "plots/"+getSpeciesName()+"_roc.png");
if (f.exists()) {
return f.getAbsolutePath();
}
return null;
}
public Map<String, String> getVariablesContributions() throws Exception{
return getOutputValues(" contribution");
}
public Map<String, String> getVariablesPermutationsImportance() throws Exception{
return getOutputValues(" permutation importance");
}
public double getPrevalence() throws Exception{
return Double.parseDouble(getOutputValues(" (average of logistic output over background sites)").values().iterator().next());
}
public double getBestThr() throws Exception{
return Double.parseDouble(getOutputValues(" training omission, predicted area and threshold value logistic threshold").values().iterator().next());
}
public String getWarnings() {
File f = new File(outputDirectory, "maxent.log");
if (f.exists()) {
StringBuffer buffer = new StringBuffer();
try {
BufferedReader br = new BufferedReader(new FileReader(f));
String line = br.readLine();
while(line!=null) {
if (line.startsWith("Warning:"))
buffer.append(line+"\n");
line = br.readLine();
}
br.close();
return buffer.toString();
} catch (Exception e) {
e.printStackTrace();
}
}
return null;
}
private Map<String, String> getOutputValues(String adjective) throws Exception{
File f = new File(outputDirectory, "maxentResults.csv");
if (f.exists()) {
Map<String, String> contributions = new LinkedHashMap<String, String>();
try {
BufferedReader br = new BufferedReader(new FileReader(f));
String headers = br.readLine();
String values = br.readLine();
List<String> heads = Transformations.parseCVSString(headers, ",");
List<String> vals = Transformations.parseCVSString(values,",");
int i = 0;
for (String head : heads) {
if (head.contains(adjective)) {
int idx = head.indexOf(" ");
String var = head;
if (idx>-1)
var = head.substring(0,idx);
contributions.put(var, vals.get(i));
}
i++;
}
br.close();
return contributions;
} catch (Exception e) {
e.printStackTrace();
throw new Exception("No occurrence records in the selected bounding box");
}
}
return null;
}
private void delFiles(String dir){
System.gc();
File f = new File(dir);
File[] list = f.listFiles();
if (list != null) {
for (File l : list) {
if (l.isFile()){
for (int j=0;j<3;j++)
l.delete();
}
}
}
}
private void delDir(String dir){
System.gc();
File f = new File(dir);
f.delete();
}
public void clean() throws Exception {
delFiles(new File(outputDirectory, "plots").getAbsolutePath());
delDir(new File(outputDirectory, "plots").getAbsolutePath());
delFiles(outputDirectory);
delDir(outputDirectory);
}
public static void main(String[] args) throws Exception {
Maxent me = new Maxent("./maxenttestfolder/occurrence_species_id0045886b_2a7c_4ede_afc4_3157c694b893_occ.csv", "./maxenttestfolder/", "./maxenttestfolder/output/", 10000, 0.5, -9999);
me.executeMaxEnt();
System.out.println("Result: "+me.getResult());
System.out.println("ROC plot: "+me.getROCPlot());
System.out.println("World plot: "+me.getWorldPlot());
System.out.println("Best Threshold: "+me.getBestThr());
System.out.println("Prevalence: "+me.getPrevalence());
System.out.println("Variables Contribution: "+me.getVariablesContributions());
System.out.println("Variables Permutations: "+me.getVariablesPermutationsImportance());
System.out.println("Omission/Commission Plot: "+me.getOmissionPlot());
System.out.println("Warnings: "+me.getWarnings());
// me.clean();
}
}

85
src/main/java/org/gcube/dataanalysis/geo/algorithms/MaxEnt4NicheModellingTransducer.java
View File

@ -8,7 +8,6 @@ import java.io.FileReader;
import java.io.FileWriter; import java.io.FileWriter;
import java.io.IOException; import java.io.IOException;
import java.util.ArrayList; import java.util.ArrayList;
import java.util.HashMap;
import java.util.LinkedHashMap; import java.util.LinkedHashMap;
import java.util.List; import java.util.List;
import java.util.UUID; import java.util.UUID;
@ -34,8 +33,12 @@ import org.gcube.dataanalysis.ecoengine.interfaces.Transducerer;
import org.gcube.dataanalysis.ecoengine.utils.DatabaseUtils; import org.gcube.dataanalysis.ecoengine.utils.DatabaseUtils;
import org.gcube.dataanalysis.ecoengine.utils.IOHelper; import org.gcube.dataanalysis.ecoengine.utils.IOHelper;
import org.gcube.dataanalysis.ecoengine.utils.ResourceFactory; import org.gcube.dataanalysis.ecoengine.utils.ResourceFactory;
import org.gcube.dataanalysis.geo.connectors.asc.ASC; import org.gcube.dataanalysis.geo.connectors.asc.AscRaster;
import org.gcube.dataanalysis.geo.connectors.asc.AscRasterReader;
import org.gcube.dataanalysis.geo.matrixmodel.ASCConverter; import org.gcube.dataanalysis.geo.matrixmodel.ASCConverter;
import org.gcube.dataanalysis.geo.matrixmodel.RasterTable;
import org.gcube.dataanalysis.geo.matrixmodel.XYExtractor;
import org.gcube.dataanalysis.geo.utils.MapUtils;
import org.hibernate.SessionFactory; import org.hibernate.SessionFactory;
import density.Maxent; import density.Maxent;
@ -50,7 +53,6 @@ public class MaxEnt4NicheModellingTransducer implements Transducerer {
static String xUR = "BBox_UpperRightLong"; static String xUR = "BBox_UpperRightLong";
static String xRes = "XResolution"; static String xRes = "XResolution";
static String yRes = "YResolution"; static String yRes = "YResolution";
static String tableName = "OutputTableName";
static String tableLabel = "OutputTableLabel"; static String tableLabel = "OutputTableLabel";
static String speciesLabel = "SpeciesName"; static String speciesLabel = "SpeciesName";
static String OccurrencesTableNameParameter = "OccurrencesTable"; static String OccurrencesTableNameParameter = "OccurrencesTable";
@ -81,8 +83,10 @@ public class MaxEnt4NicheModellingTransducer implements Transducerer {
private String variablesContributions = ""; private String variablesContributions = "";
private String variablesPermutationsImportance = ""; private String variablesPermutationsImportance = "";
private String warnings = ""; private String warnings = "";
private String layerIdentities = "";
private File warningsFile=null; private File warningsFile=null;
private File projectionFile=null; private File projectionFile=null;
private File asciimapsFile=null;
LinkedHashMap<String, Image> producedImages = new LinkedHashMap<String, Image>(); LinkedHashMap<String, Image> producedImages = new LinkedHashMap<String, Image>();
@ -95,8 +99,7 @@ public class MaxEnt4NicheModellingTransducer implements Transducerer {
private String longitudeColumn; private String longitudeColumn;
private String latitudeColumn; private String latitudeColumn;
SessionFactory dbconnection = null;
@Override @Override
public void init() throws Exception { public void init() throws Exception {
AnalysisLogger.getLogger().debug("MaxEnt: Initialization"); AnalysisLogger.getLogger().debug("MaxEnt: Initialization");
@ -107,16 +110,15 @@ public class MaxEnt4NicheModellingTransducer implements Transducerer {
public String getDescription() { public String getDescription() {
return "A Maximum-Entropy model for species habitat modeling, based on the implementation by Shapire et al. v 3.3.3k, Princeton University, http://www.cs.princeton.edu/schapire/maxent/. " + return "A Maximum-Entropy model for species habitat modeling, based on the implementation by Shapire et al. v 3.3.3k, Princeton University, http://www.cs.princeton.edu/schapire/maxent/. " +
"In this adaptation for the D4Science infrastructure, the software can accept a table produced by the Species Product Discovery service and a set of environmental layers in various formats (NetCDF, WFS, WCS, ASC, GeoTiff) via direct links or GeoExplorer UUIDs. " + "In this adaptation for the D4Science infrastructure, the software can accept a table produced by the Species Product Discovery service and a set of environmental layers in various formats (NetCDF, WFS, WCS, ASC, GeoTiff) via direct links or GeoExplorer UUIDs. " +
"Also, the user can establish the bounding box and the spatial resolution (in deg) of the training and the projection. The application will adapt the layers to that resolution if this is higher than the native one." + "The user can also establish the bounding box and the spatial resolution (in deg) of the training and the projection. The application will adapt the layers to that resolution if this is higher than the native one." +
"The output contains: a thumbnail map of the projected model, the ROC curve, the Omission/Commission chart, the raw assigned values." + "The output contains: a thumbnail map of the projected model, the ROC curve, the Omission/Commission chart, a table containing the raw assigned values, a threshold to transform the table into a 0-1 probability distribution, a report of the importance of the used layers in the model, ASCII representations of the input layers to check their alignment." +
"Other processes can be later applied to the raw values to produce a GIS map (e.g. the Statistical Manager Points to Map process)."; "Other processes can be later applied to the raw values to produce a GIS map (e.g. the Statistical Manager Points to Map process) and results can be shared. Demo video: http://goo.gl/TYYnTO";
} }
@Override @Override
public List<StatisticalType> getInputParameters() { public List<StatisticalType> getInputParameters() {
// output table parameter // output table parameter
IOHelper.addRandomStringInput(inputs, tableName, "The db name of the table to produce", "maxent_");
IOHelper.addStringInput(inputs, tableLabel, "The name of the table to produce", "maxent_"); IOHelper.addStringInput(inputs, tableLabel, "The name of the table to produce", "maxent_");
IOHelper.addStringInput(inputs, speciesLabel, "The name of the species to model and the occurrence records refer to", "generic_species"); IOHelper.addStringInput(inputs, speciesLabel, "The name of the species to model and the occurrence records refer to", "generic_species");
IOHelper.addIntegerInput(inputs, maxIterations, "The number of learning iterations of the MaxEnt algorithm", "1000"); IOHelper.addIntegerInput(inputs, maxIterations, "The number of learning iterations of the MaxEnt algorithm", "1000");
@ -135,8 +137,8 @@ public class MaxEnt4NicheModellingTransducer implements Transducerer {
IOHelper.addDoubleInput(inputs, z, "Value of Z. Default is 0, that means processing will be at surface level or at the first avaliable Z value in the layer", "0"); IOHelper.addDoubleInput(inputs, z, "Value of Z. Default is 0, that means processing will be at surface level or at the first avaliable Z value in the layer", "0");
IOHelper.addIntegerInput(inputs, t, "Time Index. The default is the first time indexed dataset", "0"); IOHelper.addIntegerInput(inputs, t, "Time Index. The default is the first time indexed dataset", "0");
IOHelper.addDoubleInput(inputs, xRes, "Projection resolution on the X axis", "1"); IOHelper.addDoubleInput(inputs, xRes, "Projection resolution on the X axis in degrees", "1");
IOHelper.addDoubleInput(inputs, yRes, "Projection resolution on the Y axis", "1"); IOHelper.addDoubleInput(inputs, yRes, "Projection resolution on the Y axis in degrees", "1");
// layers to use in the model // layers to use in the model
PrimitiveTypesList listEnvLayers = new PrimitiveTypesList(String.class.getName(), PrimitiveTypes.STRING, Layers, "The list of environmental layers to use for enriching the points. Each entry is a layer Title or UUID or HTTP link. E.g. the title or the UUID (preferred) of a layer indexed in the e-Infrastructure on GeoNetwork - You can retrieve it from GeoExplorer. Otherwise you can supply the direct HTTP link of the layer. The format will be guessed from the link. The default is GeoTiff. Supports several standards (NETCDF, WFS, WCS, ASC, GeoTiff )", false); PrimitiveTypesList listEnvLayers = new PrimitiveTypesList(String.class.getName(), PrimitiveTypes.STRING, Layers, "The list of environmental layers to use for enriching the points. Each entry is a layer Title or UUID or HTTP link. E.g. the title or the UUID (preferred) of a layer indexed in the e-Infrastructure on GeoNetwork - You can retrieve it from GeoExplorer. Otherwise you can supply the direct HTTP link of the layer. The format will be guessed from the link. The default is GeoTiff. Supports several standards (NETCDF, WFS, WCS, ASC, GeoTiff )", false);
@ -202,7 +204,6 @@ public class MaxEnt4NicheModellingTransducer implements Transducerer {
AnalysisLogger.getLogger().debug("MaxEnt: : xRes " + xResValue); AnalysisLogger.getLogger().debug("MaxEnt: : xRes " + xResValue);
// get output table value // get output table value
tableNameValue = IOHelper.getInputParameter(config, tableName);
tableLabelValue = IOHelper.getInputParameter(config, tableLabel); tableLabelValue = IOHelper.getInputParameter(config, tableLabel);
AnalysisLogger.getLogger().debug("MaxEnt: : tableName " + tableNameValue); AnalysisLogger.getLogger().debug("MaxEnt: : tableName " + tableNameValue);
AnalysisLogger.getLogger().debug("MaxEnt: : tableLabel " + tableLabelValue); AnalysisLogger.getLogger().debug("MaxEnt: : tableLabel " + tableLabelValue);
@ -230,14 +231,20 @@ public class MaxEnt4NicheModellingTransducer implements Transducerer {
String localOccurrencesFile = new File(localtempFolder, occurrencesTableName).getAbsolutePath(); String localOccurrencesFile = new File(localtempFolder, occurrencesTableName).getAbsolutePath();
String localFinalOccurrencesFile = localOccurrencesFile + "_occ.csv"; String localFinalOccurrencesFile = localOccurrencesFile + "_occ.csv";
String localAsciiMapsFile = localOccurrencesFile + "_maps.txt";
AnalysisLogger.getLogger().debug("MaxEnt: local occurrence file to produce "+localFinalOccurrencesFile); AnalysisLogger.getLogger().debug("MaxEnt: local occurrence file to produce "+localFinalOccurrencesFile);
AnalysisLogger.getLogger().debug("MaxEnt: initializing connection"); AnalysisLogger.getLogger().debug("MaxEnt: initializing connection");
dbconnection = DatabaseUtils.initDBSession(config);
// prepare input data // prepare input data
AnalysisLogger.getLogger().debug("MaxEnt: creating local file from remote table "+occurrencesTableName); AnalysisLogger.getLogger().debug("MaxEnt: creating local file from remote table "+occurrencesTableName);
DatabaseUtils.createLocalFileFromRemoteTable(localOccurrencesFile, "(select " + longitudeColumn + " as longitude," + latitudeColumn + " as latitude from " + occurrencesTableName + ")", ",", config.getDatabaseUserName(), config.getDatabasePassword(), config.getDatabaseURL());
DatabaseUtils.createLocalFileFromRemoteTable(
localOccurrencesFile, "(select " + longitudeColumn + " as longitude," + latitudeColumn + " as latitude from " + occurrencesTableName + ")", ",",
config.getParam("DatabaseUserName"),
config.getParam("DatabasePassword"),
config.getParam("DatabaseURL"));
AnalysisLogger.getLogger().debug("MaxEnt: table "+occurrencesTableName+" was dumped in file: " + localOccurrencesFile); AnalysisLogger.getLogger().debug("MaxEnt: table "+occurrencesTableName+" was dumped in file: " + localOccurrencesFile);
// write an input file for maxent // write an input file for maxent
AnalysisLogger.getLogger().debug("MaxEnt: preparing input for maxent in file "+localFinalOccurrencesFile); AnalysisLogger.getLogger().debug("MaxEnt: preparing input for maxent in file "+localFinalOccurrencesFile);
@ -246,15 +253,26 @@ public class MaxEnt4NicheModellingTransducer implements Transducerer {
// convert all the layers into file // convert all the layers into file
int layersCount = 1; int layersCount = 1;
status = 30; status = 30;
BufferedWriter mapwriter = new BufferedWriter(new FileWriter(new File(localAsciiMapsFile)));
for (String layer : layers) { for (String layer : layers) {
ASCConverter converter = new ASCConverter(config); ASCConverter converter = new ASCConverter(config);
String layerfile = new File(localtempFolder, "layer"+layersCount+ ".asc").getAbsolutePath(); String layerfile = new File(localtempFolder, "layer"+layersCount+ ".asc").getAbsolutePath();
AnalysisLogger.getLogger().debug("MaxEnt: converting " + layer +" into "+layerfile);
String converted = converter.convertToASC(layer, layerfile, time, BBxLL, BBxUR, BByLL, BByUR, zValue, xResValue, yResValue); AnalysisLogger.getLogger().debug("MaxEnt: converting " + layer +" into "+layerfile);
XYExtractor extractor = new XYExtractor(config);
double[][] values = extractor.extractXYGrid(layer, time, BBxLL, BBxUR, BByLL, BByUR, zValue, xResValue, yResValue);
mapwriter.write(MapUtils.globalASCIIMap(values));
AnalysisLogger.getLogger().debug("MaxEnt: layer name " + extractor.layerName);
layerIdentities+=layer+"="+extractor.layerName+" ";
String converted = converter.convertToASC(layerfile, values, BBxLL,BByLL,xResValue,yResValue);
AnalysisLogger.getLogger().debug("MaxEnt: converted into ASC file " + converted + " check: " + new File(converted).exists()); AnalysisLogger.getLogger().debug("MaxEnt: converted into ASC file " + converted + " check: " + new File(converted).exists());
layersCount++; layersCount++;
} }
mapwriter.close();
status = 70; status = 70;
AnalysisLogger.getLogger().debug("MaxEnt: executing MaxEnt"); AnalysisLogger.getLogger().debug("MaxEnt: executing MaxEnt");
@ -272,8 +290,8 @@ public class MaxEnt4NicheModellingTransducer implements Transducerer {
bestThreshold = me.getBestThr(); bestThreshold = me.getBestThr();
prevalenceVal = me.getPrevalence(); prevalenceVal = me.getPrevalence();
variablesContributions = me.getVariablesContributions().toString(); variablesContributions = me.getVariablesContributions().toString().replace("{", "").replace("}","");
variablesPermutationsImportance = me.getVariablesPermutationsImportance().toString(); variablesPermutationsImportance = me.getVariablesPermutationsImportance().toString().replace("{", "").replace("}","");
warnings = me.getWarnings(); warnings = me.getWarnings();
String worldFile = me.getWorldPlot(); String worldFile = me.getWorldPlot();
@ -294,6 +312,7 @@ public class MaxEnt4NicheModellingTransducer implements Transducerer {
} }
projectionFile = new File(me.getResult()); projectionFile = new File(me.getResult());
asciimapsFile = new File(localAsciiMapsFile);
AnalysisLogger.getLogger().debug("MaxEnt: Best Threshold: "+bestThreshold); AnalysisLogger.getLogger().debug("MaxEnt: Best Threshold: "+bestThreshold);
AnalysisLogger.getLogger().debug("MaxEnt: Prevalence: "+prevalenceVal); AnalysisLogger.getLogger().debug("MaxEnt: Prevalence: "+prevalenceVal);
AnalysisLogger.getLogger().debug("MaxEnt: Variables Contribution: "+variablesContributions); AnalysisLogger.getLogger().debug("MaxEnt: Variables Contribution: "+variablesContributions);
@ -302,10 +321,15 @@ public class MaxEnt4NicheModellingTransducer implements Transducerer {
AnalysisLogger.getLogger().debug("MaxEnt: Warnings file: "+warningsFile.getAbsolutePath() +" exists " + warningsFile.exists() ); AnalysisLogger.getLogger().debug("MaxEnt: Warnings file: "+warningsFile.getAbsolutePath() +" exists " + warningsFile.exists() );
AnalysisLogger.getLogger().debug("MaxEnt: Projection file: "+projectionFile.getAbsolutePath()+" exists " + projectionFile.exists() ); AnalysisLogger.getLogger().debug("MaxEnt: Projection file: "+projectionFile.getAbsolutePath()+" exists " + projectionFile.exists() );
status = 80; status = 80;
AnalysisLogger.getLogger().debug("MaxEnt: Writing the table "+tableNameValue); AnalysisLogger.getLogger().debug("MaxEnt: Generating table");
//TO DO : write a table
//write a table
AscRasterReader reader = new AscRasterReader();
AscRaster ascFile = reader.readRaster(projectionFile.getAbsolutePath());
RasterTable table = new RasterTable(BBxLL, BBxUR, BByLL, BByUR, zValue, xResValue, yResValue, ascFile.getInvertedAxisData(), config);
table.dumpGeoTable();
tableNameValue = table.getTablename();
AnalysisLogger.getLogger().debug("MaxEnt: Generated table "+tableNameValue);
// me.clean(); // me.clean();
status = 90; status = 90;
AnalysisLogger.getLogger().debug("MaxEnt: Elapsed: Whole operation completed in " + ((double) (System.currentTimeMillis() - t0) / 1000d) + "s"); AnalysisLogger.getLogger().debug("MaxEnt: Elapsed: Whole operation completed in " + ((double) (System.currentTimeMillis() - t0) / 1000d) + "s");
@ -349,8 +373,6 @@ public class MaxEnt4NicheModellingTransducer implements Transducerer {
@Override @Override
public void shutdown() { public void shutdown() {
AnalysisLogger.getLogger().debug("MaxEnt: Shutdown"); AnalysisLogger.getLogger().debug("MaxEnt: Shutdown");
if (dbconnection != null)
dbconnection.close();
} }
@Override @Override
@ -362,19 +384,14 @@ public class MaxEnt4NicheModellingTransducer implements Transducerer {
map.put("Best Threshold", new PrimitiveType(String.class.getName(), "" + bestThreshold, PrimitiveTypes.STRING, "Best Threshold", "Best threshold for transforming MaxEnt values into 0/1 probability assignments")); map.put("Best Threshold", new PrimitiveType(String.class.getName(), "" + bestThreshold, PrimitiveTypes.STRING, "Best Threshold", "Best threshold for transforming MaxEnt values into 0/1 probability assignments"));
map.put("Estimated Prevalence", new PrimitiveType(String.class.getName(), "" + prevalenceVal, PrimitiveTypes.STRING, "Estimated Prevalence", "The a posteriori estimated prevalence of the species")); map.put("Estimated Prevalence", new PrimitiveType(String.class.getName(), "" + prevalenceVal, PrimitiveTypes.STRING, "Estimated Prevalence", "The a posteriori estimated prevalence of the species"));
// map.put("Variables names in the layers", new PrimitiveType(String.class.getName(), layerIdentities, PrimitiveTypes.STRING, "Variables names in the layers", "Variables names in the layers"));
map.put("Variables Contributions", new PrimitiveType(String.class.getName(), variablesContributions, PrimitiveTypes.STRING, "Variables contributions", "The contribution of each variable to the MaxEnt values estimates")); map.put("Variables Contributions", new PrimitiveType(String.class.getName(), variablesContributions, PrimitiveTypes.STRING, "Variables contributions", "The contribution of each variable to the MaxEnt values estimates"));
map.put("Variables Permutations Importance", new PrimitiveType(String.class.getName(), variablesPermutationsImportance, PrimitiveTypes.STRING, "Variables Permutations Importance", "The importance of the permutations of the variables during the training")); map.put("Variables Permutations Importance", new PrimitiveType(String.class.getName(), variablesPermutationsImportance, PrimitiveTypes.STRING, "Variables Permutations Importance", "The importance of the permutations of the variables during the training"));
if (warningsFile!=null) if (warningsFile!=null)
map.put("Warnings", new PrimitiveType(File.class.getName(), warningsFile, PrimitiveTypes.FILE, "Warnings", "The warnings from the underlying MaxEnt model")); map.put("Warnings generated by the MaxEnt process", new PrimitiveType(File.class.getName(), warningsFile, PrimitiveTypes.FILE, "Warnings generated by the MaxEnt process", "The warnings from the underlying MaxEnt model"));
map.put("Projection File", new PrimitiveType(File.class.getName(), projectionFile, PrimitiveTypes.FILE, "Projection file", "The file containing the projection of the model")); map.put("ASCII Maps of the environmental layers for checking features aligments", new PrimitiveType(File.class.getName(), asciimapsFile, PrimitiveTypes.FILE, "ASCII Maps of the environmental layers for checking features aligments", "ASCII Maps of the environmental layers for checking features aligments"));
try {
AnalysisLogger.getLogger().debug("Test: "+new File(config.getConfigPath(),"testsspecies.png").getAbsolutePath());
producedImages.put("test",ImageTools.toImage(ImageIO.read(new File(config.getConfigPath(),"testsspecies.png"))));
} catch (IOException e) {
e.printStackTrace();
}
PrimitiveType images = new PrimitiveType(LinkedHashMap.class.getName(), producedImages, PrimitiveTypes.IMAGES, "Model performance", "Model performance and projection"); PrimitiveType images = new PrimitiveType(LinkedHashMap.class.getName(), producedImages, PrimitiveTypes.IMAGES, "Model performance", "Model performance and projection");

9
src/main/java/org/gcube/dataanalysis/geo/connectors/asc/AscRaster.java
View File

@ -156,6 +156,15 @@ public class AscRaster
return data; return data;
} }
public double[][] getInvertedAxisData()
{
double[][] dataInv = new double[data.length][data[0].length];
for (int i=data.length-1;i>=0;i--){
dataInv[data.length-1-i]=data[i];
}
return dataInv;
}
public void setValue( int row, int column, double value ) public void setValue( int row, int column, double value )
{ {
if( row < rows && column < cols ) if( row < rows && column < cols )

4
src/main/java/org/gcube/dataanalysis/geo/connectors/asc/AscRasterWriter.java
View File

@ -45,8 +45,10 @@ public class AscRasterWriter {
} }
o.println( "NODATA_value " + r.getNDATA() ); o.println( "NODATA_value " + r.getNDATA() );
for( double[] row : r.getData() ) double[][] values = r.getData();
for(int k=values.length-1;k>=0;k-- )
{ {
double[] row =values[k];
StringBuffer b = new StringBuffer(); StringBuffer b = new StringBuffer();
for( int i = 0; i < row.length; i++ ) for( int i = 0; i < row.length; i++ )
{ {

9
src/main/java/org/gcube/dataanalysis/geo/matrixmodel/ASCConverter.java
View File

@ -26,10 +26,9 @@ public class ASCConverter {
} }
public String convertToASC(String layerTitle, String outFilePath, int timeInstant, double x1, double x2, double y1, double y2, double z, double xResolution, double yResolution) throws Exception { public String convertToASC(String outFilePath, double[][] values, double x1, double y1, double xResolution, double yResolution ) throws Exception {
try { try {
double[][] values = extractor.extractXYGrid(layerTitle, timeInstant, x1, x2, y1, y2, z, xResolution, yResolution);
AscRaster raster = null; AscRaster raster = null;
if (xResolution == yResolution) if (xResolution == yResolution)
raster = new AscRaster(values, xResolution, -1, -1, x1, y1); raster = new AscRaster(values, xResolution, -1, -1, x1, y1);
@ -51,6 +50,12 @@ public class ASCConverter {
} }
public String convertToASC(String layerTitle, String outFilePath, int timeInstant, double x1, double x2, double y1, double y2, double z, double xResolution, double yResolution) throws Exception {
double[][] values = extractor.extractXYGrid(layerTitle, timeInstant, x1, x2, y1, y2, z, xResolution, yResolution);
return convertToASC(outFilePath, values, x1, y1, xResolution,yResolution);
}
public static void main(String[] args) throws Exception{ public static void main(String[] args) throws Exception{
AlgorithmConfiguration config = new AlgorithmConfiguration(); AlgorithmConfiguration config = new AlgorithmConfiguration();

13
src/main/java/org/gcube/dataanalysis/geo/test/projections/TestExtractionXYMatrixFromTable.java
View File

@ -65,6 +65,17 @@ public class TestExtractionXYMatrixFromTable {
config.setParam(TableMatrixRepresentation.filterParameter, ""); config.setParam(TableMatrixRepresentation.filterParameter, "");
} }
public static void sliceMaxEnt(AlgorithmConfiguration config) throws Exception {
config.setParam(TableMatrixRepresentation.tableNameParameter, "rstrf31af9ff13de42e583327e4ca51c38ef");
config.setParam(TableMatrixRepresentation.xDimensionColumnParameter, "x");
config.setParam(TableMatrixRepresentation.yDimensionColumnParameter, "y");
config.setParam(TableMatrixRepresentation.timeDimensionColumnParameter, "time");
config.setParam(TableMatrixRepresentation.valueDimensionColumnParameter, "fvalue");
config.setParam(TableMatrixRepresentation.filterParameter, "");
}
public static void main(String[] args) throws Exception { public static void main(String[] args) throws Exception {
AnalysisLogger.setLogger("./cfg/" + AlgorithmConfiguration.defaultLoggerFile); AnalysisLogger.setLogger("./cfg/" + AlgorithmConfiguration.defaultLoggerFile);
AlgorithmConfiguration config = new AlgorithmConfiguration(); AlgorithmConfiguration config = new AlgorithmConfiguration();
@ -76,7 +87,7 @@ public class TestExtractionXYMatrixFromTable {
config.setParam("DatabaseURL", "jdbc:postgresql://localhost/testdb"); config.setParam("DatabaseURL", "jdbc:postgresql://localhost/testdb");
config.setParam("DatabaseDriver", "org.postgresql.Driver"); config.setParam("DatabaseDriver", "org.postgresql.Driver");
config.setGcubeScope("/gcube/devsec/devVRE"); config.setGcubeScope("/gcube/devsec/devVRE");
sliceMapCreated2(config); sliceMaxEnt(config);
double resolution = 1; double resolution = 1;
FileWriter fw = new FileWriter(new File("maps.txt")); FileWriter fw = new FileWriter(new File("maps.txt"));

12
src/main/java/org/gcube/dataanalysis/geo/test/projections/TestXYExtractionConnectors.java
View File

@ -1,6 +1,7 @@
package org.gcube.dataanalysis.geo.test.projections; package org.gcube.dataanalysis.geo.test.projections;
import java.util.ArrayList; import java.io.File;
import java.io.FileWriter;
import java.util.List; import java.util.List;
import org.gcube.contentmanagement.lexicalmatcher.utils.AnalysisLogger; import org.gcube.contentmanagement.lexicalmatcher.utils.AnalysisLogger;
@ -51,6 +52,7 @@ public class TestXYExtractionConnectors {
}; };
static String[] urlToTest = { static String[] urlToTest = {
"https://dl.dropboxusercontent.com/u/12809149/layer1.asc",
"http://thredds.research-infrastructures.eu/thredds/fileServer/public/netcdf/ph.asc", "http://thredds.research-infrastructures.eu/thredds/fileServer/public/netcdf/ph.asc",
"http://thredds.research-infrastructures.eu/thredds/fileServer/public/netcdf/calcite.asc", "http://thredds.research-infrastructures.eu/thredds/fileServer/public/netcdf/calcite.asc",
"https://dl.dropboxusercontent.com/u/12809149/wind1.tif", "https://dl.dropboxusercontent.com/u/12809149/wind1.tif",
@ -67,9 +69,11 @@ public class TestXYExtractionConnectors {
"http://thredds.research-infrastructures.eu/thredds/fileServer/public/netcdf/cloudmean.asc", "http://thredds.research-infrastructures.eu/thredds/fileServer/public/netcdf/cloudmean.asc",
"http://geoserver-dev.d4science-ii.research-infrastructures.eu/geoserver/wcs/wcs?service=wcs&version=1.0.0&request=GetCoverage&coverage=aquamaps:WorldClimBio2&CRS=EPSG:4326&bbox=-180,0,180,90&width=1&height=1&format=geotiff&RESPONSE_CRS=EPSG:4326", "http://geoserver-dev.d4science-ii.research-infrastructures.eu/geoserver/wcs/wcs?service=wcs&version=1.0.0&request=GetCoverage&coverage=aquamaps:WorldClimBio2&CRS=EPSG:4326&bbox=-180,0,180,90&width=1&height=1&format=geotiff&RESPONSE_CRS=EPSG:4326",
"http://geoserver2.d4science.research-infrastructures.eu/geoserver" "http://geoserver2.d4science.research-infrastructures.eu/geoserver"
}; };
static String[] layernamesTest = { static String[] layernamesTest = {
"layer1",
"ph", "ph",
"calcite", "calcite",
"wind", "wind",
@ -99,7 +103,7 @@ public class TestXYExtractionConnectors {
config.setParam("DatabaseURL", "jdbc:postgresql://localhost/testdb"); config.setParam("DatabaseURL", "jdbc:postgresql://localhost/testdb");
config.setParam("DatabaseDriver", "org.postgresql.Driver"); config.setParam("DatabaseDriver", "org.postgresql.Driver");
config.setGcubeScope("/gcube/devsec/devVRE"); config.setGcubeScope("/gcube/devsec/devVRE");
FileWriter fw = new FileWriter(new File("mapsconnectors.txt"));
for (int t=0;t<urlToTest.length;t++){ for (int t=0;t<urlToTest.length;t++){
String layerURL = urlToTest[t]; String layerURL = urlToTest[t];
@ -135,8 +139,12 @@ public class TestXYExtractionConnectors {
// System.out.println(MapUtils.globalASCIIMap(values,step,step)); // System.out.println(MapUtils.globalASCIIMap(values,step,step));
System.out.println(MapUtils.globalASCIIMap(matrix)); System.out.println(MapUtils.globalASCIIMap(matrix));
fw.write(MapUtils.globalASCIIMap(matrix));
} }
fw.close();
} }
} }

469
src/main/java/ucar/nc2/dt/grid/GridCoordinate2D.java Normal file
View File

@ -0,0 +1,469 @@
/*
* Copyright (c) 1998 - 2009. University Corporation for Atmospheric Research/Unidata
* Portions of this software were developed by the Unidata Program at the
* University Corporation for Atmospheric Research.
*
* Access and use of this software shall impose the following obligations
* and understandings on the user. The user is granted the right, without
* any fee or cost, to use, copy, modify, alter, enhance and distribute
* this software, and any derivative works thereof, and its supporting
* documentation for any purpose whatsoever, provided that this entire
* notice appears in all copies of the software, derivative works and
* supporting documentation. Further, UCAR requests that the user credit
* UCAR/Unidata in any publications that result from the use of this
* software or in any product that includes this software. The names UCAR
* and/or Unidata, however, may not be used in any advertising or publicity
* to endorse or promote any products or commercial entity unless specific
* written permission is obtained from UCAR/Unidata. The user also
* understands that UCAR/Unidata is not obligated to provide the user with
* any support, consulting, training or assistance of any kind with regard
* to the use, operation and performance of this software nor to provide
* the user with any updates, revisions, new versions or "bug fixes."
*
* THIS SOFTWARE IS PROVIDED BY UCAR/UNIDATA "AS IS" AND ANY EXPRESS OR
* IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
* WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
* DISCLAIMED. IN NO EVENT SHALL UCAR/UNIDATA BE LIABLE FOR ANY SPECIAL,
* INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES WHATSOEVER RESULTING
* FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN ACTION OF CONTRACT,
* NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF OR IN CONNECTION
* WITH THE ACCESS, USE OR PERFORMANCE OF THIS SOFTWARE.
*/
package ucar.nc2.dt.grid;
import ucar.nc2.dataset.CoordinateAxis2D;
import ucar.nc2.dataset.CoordinateAxis;
import ucar.nc2.dt.GridCoordSystem;
import ucar.ma2.ArrayDouble;
import ucar.ma2.MAMath;
import java.io.IOException;
/**
* 2D Coordinate System has lat(x,y) and lon(x,y).
* This class implements finding the index (i,j) from (lat, lon) coord.
* This is for "one-off" computation, not a systematic lookup table for all points in a pixel array.
* Hueristically searches the 2D space for the cell tha contains the point.
*
* @author caron
* @since Jul 10, 2009
*/
public class GridCoordinate2D {
static private boolean debug = false;
static private org.slf4j.Logger log = org.slf4j.LoggerFactory.getLogger(GridCoordinate2D.class);
private CoordinateAxis2D latCoord, lonCoord;
private ArrayDouble.D2 latEdge, lonEdge;
private MAMath.MinMax latMinMax, lonMinMax;
int nrows, ncols;
public GridCoordinate2D(CoordinateAxis2D latCoord, CoordinateAxis2D lonCoord) {
this.latCoord = latCoord;
this.lonCoord = lonCoord;
assert latCoord.getRank() == 2;
assert lonCoord.getRank() == 2;
int[] shape = latCoord.getShape();
nrows = shape[0];
ncols = shape[1];
}
private void findBounds() {
if (lonMinMax != null) return;
if ((lonEdge==null) && (latEdge==null)){
lonEdge = CoordinateAxis2D.makeXEdges(lonCoord.getMidpoints());
latEdge = CoordinateAxis2D.makeYEdges(latCoord.getMidpoints());
// assume missing values have been converted to NaNs
latMinMax = MAMath.getMinMax(latEdge);
lonMinMax = MAMath.getMinMax(lonEdge);
}
if (debug)
System.out.printf("Bounds (%d %d): lat= (%f,%f) lon = (%f,%f) %n", nrows, ncols, latMinMax.min, latMinMax.max, lonMinMax.min, lonMinMax.max);
}
// brute force
public boolean findCoordElementForce(double wantLat, double wantLon, int[] rectIndex) {
findBounds();
if (wantLat < latMinMax.min) return false;
if (wantLat > latMinMax.max) return false;
if (wantLon < lonMinMax.min) return false;
if (wantLon > lonMinMax.max) return false;
boolean saveDebug = debug;
debug = false;
for (int row=0; row<nrows; row++) {
for (int col=0; col<ncols; col++) {
rectIndex[0] = row;
rectIndex[1] = col;
if (contains(wantLat, wantLon, rectIndex)) {
debug = saveDebug;
return true;
}
}
}
//debug = saveDebug;
return false;
}
public boolean findCoordElement(double wantLat, double wantLon, int[] rectIndex) {
return findCoordElementNoForce(wantLat, wantLon,rectIndex);
}
/**
* Find the best index for the given lat,lon point.
* @param wantLat lat of point
* @param wantLon lon of point
* @param rectIndex return (row,col) index, or best guess here. may not be null
*
* @return false if not in the grid.
*/
public boolean findCoordElementNoForce(double wantLat, double wantLon, int[] rectIndex) {
findBounds();
if (wantLat < latMinMax.min) return false;
if (wantLat > latMinMax.max) return false;
if (wantLon < lonMinMax.min) return false;
if (wantLon > lonMinMax.max) return false;
double gradientLat = (latMinMax.max - latMinMax.min) / nrows;
double gradientLon = (lonMinMax.max - lonMinMax.min) / ncols;
double diffLat = wantLat - latMinMax.min;
double diffLon = wantLon - lonMinMax.min;
// initial guess
rectIndex[0] = (int) Math.round(diffLat / gradientLat); // row
rectIndex[1] =(int) Math.round(diffLon / gradientLon); // col
int count = 0;
while (true) {
count++;
if (debug) System.out.printf("%nIteration %d %n", count);
if (contains(wantLat, wantLon, rectIndex))
return true;
if (!jump2(wantLat, wantLon, rectIndex)) return false;
// bouncing around
if (count > 10) {
// last ditch attempt
return incr(wantLat, wantLon, rectIndex);
//if (!ok)
// log.error("findCoordElement didnt converge lat,lon = "+wantLat+" "+ wantLon);
//return ok;
}
}
}
/**
* Is the point (lat,lon) contained in the (row, col) rectangle ?
*
* @param wantLat lat of point
* @param wantLon lon of point
* @param rectIndex rectangle row, col, will be clipped to [0, nrows), [0, ncols)
* @return true if contained
*/
private boolean containsOld(double wantLat, double wantLon, int[] rectIndex) {
rectIndex[0] = Math.max( Math.min(rectIndex[0], nrows-1), 0);
rectIndex[1] = Math.max( Math.min(rectIndex[1], ncols-1), 0);
int row = rectIndex[0];
int col = rectIndex[1];
if (debug) System.out.printf(" (%d,%d) contains (%f,%f) in (lat=%f %f) (lon=%f %f) ?%n",
rectIndex[0], rectIndex[1], wantLat, wantLon,
latEdge.get(row, col), latEdge.get(row + 1, col), lonEdge.get(row, col), lonEdge.get(row, col + 1));
if (wantLat < latEdge.get(row, col)) return false;
if (wantLat > latEdge.get(row + 1, col)) return false;
if (wantLon < lonEdge.get(row, col)) return false;
if (wantLon > lonEdge.get(row, col + 1)) return false;
return true;
}
/**
* Is the point (lat,lon) contained in the (row, col) rectangle ?
*
* @param wantLat lat of point
* @param wantLon lon of point
* @param rectIndex rectangle row, col, will be clipped to [0, nrows), [0, ncols)
* @return true if contained
*/
/*
http://mathforum.org/library/drmath/view/54386.html
Given any three points on the plane (x0,y0), (x1,y1), and
(x2,y2), the area of the triangle determined by them is
given by the following formula:
1 | x0 y0 1 |
A = - | x1 y1 1 |,
2 | x2 y2 1 |
where the vertical bars represent the determinant.
the value of the expression above is:
(.5)(x1*y2 - y1*x2 -x0*y2 + y0*x2 + x0*y1 - y0*x1)
The amazing thing is that A is positive if the three points are
taken in a counter-clockwise orientation, and negative otherwise.
To be inside a rectangle (or any convex body), as you trace
around in a clockwise direction from p1 to p2 to p3 to p4 and
back to p1, the "areas" of triangles p1 p2 p, p2 p3 p, p3 p4 p,
and p4 p1 p must all be positive. If you don't know the
orientation of your rectangle, then they must either all be
positive or all be negative.
this method works for arbitrary convex regions oo the plane.
*/
private boolean contains(double wantLat, double wantLon, int[] rectIndex) {
rectIndex[0] = Math.max( Math.min(rectIndex[0], nrows-1), 0);
rectIndex[1] = Math.max( Math.min(rectIndex[1], ncols-1), 0);
int row = rectIndex[0];
int col = rectIndex[1];
double x1 = lonEdge.get(row, col);
double y1 = latEdge.get(row, col);
double x2 = lonEdge.get(row, col+1);
double y2 = latEdge.get(row, col+1);
double x3 = lonEdge.get(row+1, col+1);
double y3 = latEdge.get(row+1, col+1);
double x4 = lonEdge.get(row+1, col);
double y4 = latEdge.get(row+1, col);
// must all have same determinate sign
/*
boolean sign = detIsPositive(x1, y1, x2, y2, wantLon, wantLat);
if (sign != detIsPositive(x2, y2, x3, y3, wantLon, wantLat)) return false;
if (sign != detIsPositive(x3, y3, x4, y4, wantLon, wantLat)) return false;
if (sign != detIsPositive(x4, y4, x1, y1, wantLon, wantLat)) return false;
*/
return true;
}
private boolean detIsPositive(double x0, double y0, double x1, double y1, double x2, double y2) {
double det = (x1*y2 - y1*x2 -x0*y2 + y0*x2 + x0*y1 - y0*x1);
if (det == 0)
System.out.printf("determinate = 0%n");
return det > 0;
}
private boolean jumpOld(double wantLat, double wantLon, int[] rectIndex) {
int row = Math.max( Math.min(rectIndex[0], nrows-1), 0);
int col = Math.max( Math.min(rectIndex[1], ncols-1), 0);
double gradientLat = latEdge.get(row + 1, col) - latEdge.get(row, col);
double gradientLon = lonEdge.get(row, col + 1) - lonEdge.get(row, col);
double diffLat = wantLat - latEdge.get(row, col);
double diffLon = wantLon - lonEdge.get(row, col);
int drow = (int) Math.round(diffLat / gradientLat);
int dcol = (int) Math.round(diffLon / gradientLon);
if (debug) System.out.printf(" jump from %d %d (grad=%f %f) (diff=%f %f) (delta=%d %d)",
row, col, gradientLat, gradientLon,
diffLat, diffLon, drow, dcol);
if ((drow == 0) && (dcol == 0)) {
if (debug) System.out.printf("%n incr:");
return incr(wantLat, wantLon, rectIndex);
} else {
rectIndex[0] = Math.max( Math.min(row + drow, nrows-1), 0);
rectIndex[1] = Math.max( Math.min(col + dcol, ncols-1), 0);
if (debug) System.out.printf(" to (%d %d)%n", rectIndex[0], rectIndex[1]);
if ((row == rectIndex[0]) && (col == rectIndex[1])) return false; // nothing has changed
}
return true;
}
/**
* jump to a new row,col
* @param wantLat want this lat
* @param wantLon want this lon
* @param rectIndex starting row, col and replaced by new guess
* @return true if new guess, false if failure
*/
/*
choose x, y such that (matrix multiply) :
(wantx) = (fxx fxy) (x)
(wanty) (fyx fyy) (y)
where fxx = d(fx)/dx ~= delta lon in lon direction
where fxy = d(fx)/dy ~= delta lon in lat direction
where fyx = d(fy)/dx ~= delta lat in lon direction
where fyy = d(fy)/dy ~= delta lat in lat direction
2 linear equations in 2 unknowns, solve in usual way
*/
private boolean jump2(double wantLat, double wantLon, int[] rectIndex) {
int row = Math.max( Math.min(rectIndex[0], nrows-1), 0);
int col = Math.max( Math.min(rectIndex[1], ncols-1), 0);
double dlatdy = latEdge.get(row + 1, col) - latEdge.get(row, col);
double dlatdx = latEdge.get(row, col+1) - latEdge.get(row, col);
double dlondx = lonEdge.get(row, col + 1) - lonEdge.get(row, col);
double dlondy = lonEdge.get(row + 1, col) - lonEdge.get(row, col);
double diffLat = wantLat - latEdge.get(row, col);
double diffLon = wantLon - lonEdge.get(row, col);
// solve for dlon
double dx = (diffLon - dlondy * diffLat / dlatdy) / (dlondx - dlatdx * dlondy / dlatdy);
// double dy = (diffLat - dlatdx * diffLon / dlondx) / (dlatdy - dlatdx * dlondy / dlondx);
double dy = (diffLat - dlatdx * dx) / dlatdy;
if (debug) System.out.printf(" jump from %d %d (dlondx=%f dlondy=%f dlatdx=%f dlatdy=%f) (diffLat,Lon=%f %f) (deltalat,Lon=%f %f)",
row, col, dlondx, dlondy, dlatdx, dlatdy,
diffLat, diffLon, dy, dx);
int drow = (int) Math.round(dy);
int dcol = (int) Math.round(dx);
if ((drow == 0) && (dcol == 0)) {
if (debug) System.out.printf("%n incr:");
return incr(wantLat, wantLon, rectIndex);
} else {
rectIndex[0] = Math.max( Math.min(row + drow, nrows-1), 0);
rectIndex[1] = Math.max( Math.min(col + dcol, ncols-1), 0);
if (debug) System.out.printf(" to (%d %d)%n", rectIndex[0], rectIndex[1]);
if ((row == rectIndex[0]) && (col == rectIndex[1])) return false; // nothing has changed
}
return true;
}
private boolean incr(double wantLat, double wantLon, int[] rectIndex) {
int row = rectIndex[0];
int col = rectIndex[1];
double diffLat = wantLat - latEdge.get(row, col);
double diffLon = wantLon - lonEdge.get(row, col);
if (Math.abs(diffLat) > Math.abs(diffLon)) { // try lat first
rectIndex[0] = row + ((diffLat > 0) ? 1 : -1);
if (contains(wantLat, wantLon, rectIndex)) return true;
rectIndex[1] = col + ((diffLon > 0) ? 1 : -1);
if (contains(wantLat, wantLon, rectIndex)) return true;
} else {
rectIndex[1] = col + ((diffLon > 0) ? 1 : -1);
if (contains(wantLat, wantLon, rectIndex)) return true;
rectIndex[0] = row + ((diffLat > 0) ? 1 : -1);
if (contains(wantLat, wantLon, rectIndex)) return true;
}
// back to original, do box search
rectIndex[0] = row;
rectIndex[1] = col;
return box9(wantLat, wantLon, rectIndex);
}
// we think its got to be in one of the 9 boxes around rectIndex
private boolean box9(double wantLat, double wantLon, int[] rectIndex) {
int row = rectIndex[0];
int minrow = Math.max(row-1, 0);
int maxrow = Math.min(row+1, nrows);
int col = rectIndex[1];
int mincol= Math.max(col-1, 0);
int maxcol = Math.min(col+1, ncols);
if (debug) System.out.printf("%n box9:");
for (int i=minrow; i<=maxrow; i++)
for (int j=mincol; j<=maxcol; j++) {
rectIndex[0] = i;
rectIndex[1] = j;
if (contains(wantLat, wantLon, rectIndex)) return true;
}
return false;
}
public static void doOne(GridCoordinate2D g2d, double wantLat, double wantLon) {
int[] result = new int[2];
if (g2d.findCoordElementForce(wantLat, wantLon, result))
System.out.printf("Brute (%f %f) == (%d %d) %n", wantLat, wantLon, result[0], result[1]);
else {
System.out.printf("Brute (%f %f) FAIL", wantLat, wantLon);
return;
}
if (g2d.findCoordElement(wantLat, wantLon, result))
System.out.printf("(%f %f) == (%d %d) %n", wantLat, wantLon, result[0], result[1]);
else
System.out.printf("(%f %f) FAIL %n", wantLat, wantLon);
System.out.printf("----------------------------------------%n");
}
public static void test1() throws IOException {
String filename = "D:/work/asaScience/EGM200_3.ncml";
GridDataset gds = GridDataset.open(filename);
GeoGrid grid = gds.findGridByName("u_wind");
GridCoordSystem gcs = grid.getCoordinateSystem();
CoordinateAxis lonAxis = gcs.getXHorizAxis();
assert lonAxis instanceof CoordinateAxis2D;
CoordinateAxis latAxis = gcs.getYHorizAxis();
assert latAxis instanceof CoordinateAxis2D;
GridCoordinate2D g2d = new GridCoordinate2D((CoordinateAxis2D) latAxis, (CoordinateAxis2D) lonAxis);
doOne(g2d, 35.0, -6.0);
doOne(g2d, 34.667302, -5.008376); // FAIL
doOne(g2d, 34.667303, -6.394240);
doOne(g2d, 36.6346, -5.0084);
doOne(g2d, 36.6346, -6.394240);
gds.close();
}
public static void test2() throws IOException {
String filename = "C:/data/fmrc/apex_fmrc/Run_20091025_0000.nc";
GridDataset gds = GridDataset.open(filename);
GeoGrid grid = gds.findGridByName("temp");
GridCoordSystem gcs = grid.getCoordinateSystem();
CoordinateAxis lonAxis = gcs.getXHorizAxis();
assert lonAxis instanceof CoordinateAxis2D;
CoordinateAxis latAxis = gcs.getYHorizAxis();
assert latAxis instanceof CoordinateAxis2D;
GridCoordinate2D g2d = new GridCoordinate2D((CoordinateAxis2D) latAxis, (CoordinateAxis2D) lonAxis);
doOne(g2d, 40.166959,-73.954234);
gds.close();
}
public static void test3() throws IOException {
String filename = "/data/testdata/cdmUnitTest/fmrc/rtofs/ofs.20091122/ofs_atl.t00z.F024.grb.grib2";
GridDataset gds = GridDataset.open(filename);
GeoGrid grid = gds.findGridByName("Sea_Surface_Height_Relative_to_Geoid");
GridCoordSystem gcs = grid.getCoordinateSystem();
CoordinateAxis lonAxis = gcs.getXHorizAxis();
assert lonAxis instanceof CoordinateAxis2D;
CoordinateAxis latAxis = gcs.getYHorizAxis();
assert latAxis instanceof CoordinateAxis2D;
GridCoordinate2D g2d = new GridCoordinate2D((CoordinateAxis2D) latAxis, (CoordinateAxis2D) lonAxis);
doOne(g2d, -15.554099426977835, -0.7742870290336263);
gds.close();
}
public static void main(String[] args) throws IOException {
test3();
}
}