From 96fd4efcb07a13319e97b7e560cd2cd0afe0fe82 Mon Sep 17 00:00:00 2001 From: Gianpaolo Coro Date: Thu, 2 Oct 2014 12:24:14 +0000 Subject: [PATCH] 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 --- src/main/java/density/Maxent.java | 301 +++++++++++ .../MaxEnt4NicheModellingTransducer.java | 85 ++-- .../geo/connectors/asc/AscRaster.java | 9 + .../geo/connectors/asc/AscRasterWriter.java | 4 +- .../geo/matrixmodel/ASCConverter.java | 9 +- .../TestExtractionXYMatrixFromTable.java | 13 +- .../TestXYExtractionConnectors.java | 12 +- .../ucar/nc2/dt/grid/GridCoordinate2D.java | 469 ++++++++++++++++++ 8 files changed, 862 insertions(+), 40 deletions(-) create mode 100644 src/main/java/density/Maxent.java create mode 100644 src/main/java/ucar/nc2/dt/grid/GridCoordinate2D.java diff --git a/src/main/java/density/Maxent.java b/src/main/java/density/Maxent.java new file mode 100644 index 0000000..f75bc0c --- /dev/null +++ b/src/main/java/density/Maxent.java @@ -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 envVariables = new ArrayList(); + + 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 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 getVariablesContributions() throws Exception{ + return getOutputValues(" contribution"); + } + + public Map 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 getOutputValues(String adjective) throws Exception{ + + File f = new File(outputDirectory, "maxentResults.csv"); + if (f.exists()) { + Map contributions = new LinkedHashMap(); + try { + BufferedReader br = new BufferedReader(new FileReader(f)); + String headers = br.readLine(); + String values = br.readLine(); + List heads = Transformations.parseCVSString(headers, ","); + List 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(); + } + +} diff --git a/src/main/java/org/gcube/dataanalysis/geo/algorithms/MaxEnt4NicheModellingTransducer.java b/src/main/java/org/gcube/dataanalysis/geo/algorithms/MaxEnt4NicheModellingTransducer.java index d3c4453..0c8cd72 100644 --- a/src/main/java/org/gcube/dataanalysis/geo/algorithms/MaxEnt4NicheModellingTransducer.java +++ b/src/main/java/org/gcube/dataanalysis/geo/algorithms/MaxEnt4NicheModellingTransducer.java @@ -8,7 +8,6 @@ import java.io.FileReader; import java.io.FileWriter; import java.io.IOException; import java.util.ArrayList; -import java.util.HashMap; import java.util.LinkedHashMap; import java.util.List; 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.IOHelper; 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.RasterTable; +import org.gcube.dataanalysis.geo.matrixmodel.XYExtractor; +import org.gcube.dataanalysis.geo.utils.MapUtils; import org.hibernate.SessionFactory; import density.Maxent; @@ -50,7 +53,6 @@ public class MaxEnt4NicheModellingTransducer implements Transducerer { static String xUR = "BBox_UpperRightLong"; static String xRes = "XResolution"; static String yRes = "YResolution"; - static String tableName = "OutputTableName"; static String tableLabel = "OutputTableLabel"; static String speciesLabel = "SpeciesName"; static String OccurrencesTableNameParameter = "OccurrencesTable"; @@ -81,8 +83,10 @@ public class MaxEnt4NicheModellingTransducer implements Transducerer { private String variablesContributions = ""; private String variablesPermutationsImportance = ""; private String warnings = ""; + private String layerIdentities = ""; private File warningsFile=null; private File projectionFile=null; + private File asciimapsFile=null; LinkedHashMap producedImages = new LinkedHashMap(); @@ -95,8 +99,7 @@ public class MaxEnt4NicheModellingTransducer implements Transducerer { private String longitudeColumn; private String latitudeColumn; - SessionFactory dbconnection = null; - + @Override public void init() throws Exception { AnalysisLogger.getLogger().debug("MaxEnt: Initialization"); @@ -107,16 +110,15 @@ public class MaxEnt4NicheModellingTransducer implements Transducerer { 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/. " + "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 output contains: a thumbnail map of the projected model, the ROC curve, the Omission/Commission chart, the raw assigned values." + - "Other processes can be later applied to the raw values to produce a GIS map (e.g. the Statistical Manager Points to Map process)."; + "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, 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) and results can be shared. Demo video: http://goo.gl/TYYnTO"; } @Override public List getInputParameters() { // 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, 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"); @@ -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.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, yRes, "Projection resolution on the Y 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 in degrees", "1"); // 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); @@ -202,7 +204,6 @@ public class MaxEnt4NicheModellingTransducer implements Transducerer { AnalysisLogger.getLogger().debug("MaxEnt: : xRes " + xResValue); // get output table value - tableNameValue = IOHelper.getInputParameter(config, tableName); tableLabelValue = IOHelper.getInputParameter(config, tableLabel); AnalysisLogger.getLogger().debug("MaxEnt: : tableName " + tableNameValue); AnalysisLogger.getLogger().debug("MaxEnt: : tableLabel " + tableLabelValue); @@ -230,14 +231,20 @@ public class MaxEnt4NicheModellingTransducer implements Transducerer { String localOccurrencesFile = new File(localtempFolder, occurrencesTableName).getAbsolutePath(); String localFinalOccurrencesFile = localOccurrencesFile + "_occ.csv"; - + String localAsciiMapsFile = localOccurrencesFile + "_maps.txt"; + AnalysisLogger.getLogger().debug("MaxEnt: local occurrence file to produce "+localFinalOccurrencesFile); AnalysisLogger.getLogger().debug("MaxEnt: initializing connection"); - dbconnection = DatabaseUtils.initDBSession(config); - + // prepare input data 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); // write an input file for maxent 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 int layersCount = 1; status = 30; + BufferedWriter mapwriter = new BufferedWriter(new FileWriter(new File(localAsciiMapsFile))); + for (String layer : layers) { ASCConverter converter = new ASCConverter(config); 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()); layersCount++; } + + mapwriter.close(); status = 70; AnalysisLogger.getLogger().debug("MaxEnt: executing MaxEnt"); @@ -272,8 +290,8 @@ public class MaxEnt4NicheModellingTransducer implements Transducerer { bestThreshold = me.getBestThr(); prevalenceVal = me.getPrevalence(); - variablesContributions = me.getVariablesContributions().toString(); - variablesPermutationsImportance = me.getVariablesPermutationsImportance().toString(); + variablesContributions = me.getVariablesContributions().toString().replace("{", "").replace("}",""); + variablesPermutationsImportance = me.getVariablesPermutationsImportance().toString().replace("{", "").replace("}",""); warnings = me.getWarnings(); String worldFile = me.getWorldPlot(); @@ -294,6 +312,7 @@ public class MaxEnt4NicheModellingTransducer implements Transducerer { } projectionFile = new File(me.getResult()); + asciimapsFile = new File(localAsciiMapsFile); AnalysisLogger.getLogger().debug("MaxEnt: Best Threshold: "+bestThreshold); AnalysisLogger.getLogger().debug("MaxEnt: Prevalence: "+prevalenceVal); 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: Projection file: "+projectionFile.getAbsolutePath()+" exists " + projectionFile.exists() ); status = 80; - AnalysisLogger.getLogger().debug("MaxEnt: Writing the table "+tableNameValue); - - //TO DO : write a table + AnalysisLogger.getLogger().debug("MaxEnt: Generating 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(); status = 90; 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 public void shutdown() { AnalysisLogger.getLogger().debug("MaxEnt: Shutdown"); - if (dbconnection != null) - dbconnection.close(); } @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("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 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) - 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")); - - 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(); - } + 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")); PrimitiveType images = new PrimitiveType(LinkedHashMap.class.getName(), producedImages, PrimitiveTypes.IMAGES, "Model performance", "Model performance and projection"); diff --git a/src/main/java/org/gcube/dataanalysis/geo/connectors/asc/AscRaster.java b/src/main/java/org/gcube/dataanalysis/geo/connectors/asc/AscRaster.java index 0072a1c..23e9cc2 100644 --- a/src/main/java/org/gcube/dataanalysis/geo/connectors/asc/AscRaster.java +++ b/src/main/java/org/gcube/dataanalysis/geo/connectors/asc/AscRaster.java @@ -156,6 +156,15 @@ public class AscRaster 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 ) { if( row < rows && column < cols ) diff --git a/src/main/java/org/gcube/dataanalysis/geo/connectors/asc/AscRasterWriter.java b/src/main/java/org/gcube/dataanalysis/geo/connectors/asc/AscRasterWriter.java index 8f5a492..2f40e84 100644 --- a/src/main/java/org/gcube/dataanalysis/geo/connectors/asc/AscRasterWriter.java +++ b/src/main/java/org/gcube/dataanalysis/geo/connectors/asc/AscRasterWriter.java @@ -45,8 +45,10 @@ public class AscRasterWriter { } 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(); for( int i = 0; i < row.length; i++ ) { diff --git a/src/main/java/org/gcube/dataanalysis/geo/matrixmodel/ASCConverter.java b/src/main/java/org/gcube/dataanalysis/geo/matrixmodel/ASCConverter.java index b6ff27d..32a26d3 100644 --- a/src/main/java/org/gcube/dataanalysis/geo/matrixmodel/ASCConverter.java +++ b/src/main/java/org/gcube/dataanalysis/geo/matrixmodel/ASCConverter.java @@ -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 { - double[][] values = extractor.extractXYGrid(layerTitle, timeInstant, x1, x2, y1, y2, z, xResolution, yResolution); AscRaster raster = null; if (xResolution == yResolution) 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{ AlgorithmConfiguration config = new AlgorithmConfiguration(); diff --git a/src/main/java/org/gcube/dataanalysis/geo/test/projections/TestExtractionXYMatrixFromTable.java b/src/main/java/org/gcube/dataanalysis/geo/test/projections/TestExtractionXYMatrixFromTable.java index d1ae399..7c29009 100644 --- a/src/main/java/org/gcube/dataanalysis/geo/test/projections/TestExtractionXYMatrixFromTable.java +++ b/src/main/java/org/gcube/dataanalysis/geo/test/projections/TestExtractionXYMatrixFromTable.java @@ -65,6 +65,17 @@ public class TestExtractionXYMatrixFromTable { 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 { AnalysisLogger.setLogger("./cfg/" + AlgorithmConfiguration.defaultLoggerFile); AlgorithmConfiguration config = new AlgorithmConfiguration(); @@ -76,7 +87,7 @@ public class TestExtractionXYMatrixFromTable { config.setParam("DatabaseURL", "jdbc:postgresql://localhost/testdb"); config.setParam("DatabaseDriver", "org.postgresql.Driver"); config.setGcubeScope("/gcube/devsec/devVRE"); - sliceMapCreated2(config); + sliceMaxEnt(config); double resolution = 1; FileWriter fw = new FileWriter(new File("maps.txt")); diff --git a/src/main/java/org/gcube/dataanalysis/geo/test/projections/TestXYExtractionConnectors.java b/src/main/java/org/gcube/dataanalysis/geo/test/projections/TestXYExtractionConnectors.java index b76e19e..65ed708 100644 --- a/src/main/java/org/gcube/dataanalysis/geo/test/projections/TestXYExtractionConnectors.java +++ b/src/main/java/org/gcube/dataanalysis/geo/test/projections/TestXYExtractionConnectors.java @@ -1,6 +1,7 @@ package org.gcube.dataanalysis.geo.test.projections; -import java.util.ArrayList; +import java.io.File; +import java.io.FileWriter; import java.util.List; import org.gcube.contentmanagement.lexicalmatcher.utils.AnalysisLogger; @@ -51,6 +52,7 @@ public class TestXYExtractionConnectors { }; 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/calcite.asc", "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://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" + }; static String[] layernamesTest = { + "layer1", "ph", "calcite", "wind", @@ -99,7 +103,7 @@ public class TestXYExtractionConnectors { config.setParam("DatabaseURL", "jdbc:postgresql://localhost/testdb"); config.setParam("DatabaseDriver", "org.postgresql.Driver"); config.setGcubeScope("/gcube/devsec/devVRE"); - + FileWriter fw = new FileWriter(new File("mapsconnectors.txt")); for (int t=0;t 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 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(); + } + +}