package org.gcube.dataanalysis.geo.connectors.netcdf; import java.util.ArrayList; import java.util.Formatter; import java.util.LinkedHashMap; import java.util.List; import org.gcube.contentmanagement.graphtools.utils.MathFunctions; import org.gcube.contentmanagement.lexicalmatcher.utils.AnalysisLogger; import org.gcube.dataanalysis.ecoengine.utils.Tuple; import org.gcube.dataanalysis.geo.utils.VectorOperations; import ucar.ma2.Array; import ucar.ma2.ArrayByte; import ucar.ma2.ArrayDouble; import ucar.ma2.ArrayFloat; import ucar.ma2.ArrayInt; import ucar.ma2.ArrayLong; import ucar.ma2.Range; import ucar.ma2.StructureData; import ucar.ma2.StructureMembers.Member; import ucar.nc2.constants.FeatureType; import ucar.nc2.dataset.CoordinateAxis; import ucar.nc2.dataset.CoordinateAxis1DTime; import ucar.nc2.dt.GridCoordSystem; import ucar.nc2.dt.GridDatatype; import ucar.nc2.dt.grid.GridDataset; import ucar.nc2.ft.FeatureCollection; import ucar.nc2.ft.FeatureDataset; import ucar.nc2.ft.FeatureDatasetFactoryManager; import ucar.nc2.ft.PointFeatureCollection; import ucar.nc2.ft.PointFeatureIterator; import ucar.nc2.ft.point.PointDatasetImpl; import ucar.nc2.ft.point.standard.StandardPointCollectionImpl; import ucar.unidata.geoloc.LatLonPoint; import ucar.unidata.geoloc.LatLonPointImpl; import ucar.unidata.geoloc.LatLonRect; public class NetCDFDataExplorer { // http://thredds.research-infrastructures.eu:8080/thredds/catalog/public/netcdf/catalog.xml public static String timePrefix = "time:"; public NetCDFDataExplorer(String openDapLink, String layer) { calcZRange(openDapLink, layer); } public List retrieveDataFromNetCDF(String openDapLink, String layer, int time, List> triplets) { try { List values = new ArrayList(); if (isGridDataset(openDapLink)) { AnalysisLogger.getLogger().debug("Managing Grid File"); return manageGridDataset(layer, openDapLink, time, triplets); } /* * else if (isPointDataset(openDapLink)) { AnalysisLogger.getLogger().debug("Managing Points File"); } */ else AnalysisLogger.getLogger().debug("Warning: the NETCDF file is of an unknown type"); return values; } catch (Exception e) { AnalysisLogger.getLogger().debug("ERROR: " + e.getMessage()); AnalysisLogger.getLogger().debug(e); // e.printStackTrace(); return null; } } public double minZ = 0; public double maxZ = 0; public void calcZRange(String openDapLink, String layer) { try { if (isGridDataset(openDapLink)) { gds = ucar.nc2.dt.grid.GridDataset.open(openDapLink); List gridTypes = gds.getGrids(); for (GridDatatype gdt : gridTypes) { // AnalysisLogger.getLogger().debug("Inside File - layer name: " + gdt.getFullName()); AnalysisLogger.getLogger().debug("Inside File - layer name: " + gdt.getName()); if (layer.equalsIgnoreCase(gdt.getName())) { CoordinateAxis zAxis = gdt.getCoordinateSystem().getVerticalAxis(); minZ = zAxis.getMinValue(); maxZ = zAxis.getMaxValue(); break; } } } } catch (Exception e) { e.printStackTrace(); AnalysisLogger.getLogger().debug("NetCDF Explorer Error:" + e.getLocalizedMessage()); AnalysisLogger.getLogger().debug(e); } } // A GridDatatype is like a specialized Variable that explicitly handles X,Y,Z,T dimensions GridDataset gds; public List manageGridDataset(String layer, String filename, int time, List> triplets) throws Exception { List values = new ArrayList(); if (gds == null) gds = ucar.nc2.dt.grid.GridDataset.open(filename); List gridTypes = gds.getGrids(); for (GridDatatype gdt : gridTypes) { AnalysisLogger.getLogger().debug("Inside File - layer name: " + gdt.getName() + " layer to find " + layer); // if the layer is an HTTP link then take the first innser layer if (layer.equalsIgnoreCase(gdt.getName()) || layer.toLowerCase().startsWith("http:")) { AnalysisLogger.getLogger().debug("Found layer " + layer + " inside file"); GridDatatype grid = gds.findGridDatatype(gdt.getName()); CoordinateAxis zAxis = gdt.getCoordinateSystem().getVerticalAxis(); CoordinateAxis xAxis = gdt.getCoordinateSystem().getXHorizAxis(); CoordinateAxis yAxis = gdt.getCoordinateSystem().getYHorizAxis(); double resolutionZ = 0; try { resolutionZ = Math.abs((double) (zAxis.getMaxValue() - zAxis.getMinValue()) / (double) zAxis.getShape()[0]); AnalysisLogger.getLogger().debug("Zmin:" + zAxis.getMinValue() + " Zmax:" + zAxis.getMaxValue()); } catch (Exception e) { } GridCoordSystem gcs = grid.getCoordinateSystem(); int tsize = triplets.size(); long t01 = System.currentTimeMillis(); LatLonRect llr = null; AnalysisLogger.getLogger().debug("Extracting subset..."); GridDatatype gdtsub = grid.makeSubset(new Range(time, time), null, llr, 1, 1, 1); Array data = gdtsub.readVolumeData(time); // note order is t, z, y, x int[] shapeD = data.getShape(); int zD = 0; int xD = 0; int yD = 0; if (shapeD.length > 2) { zD = shapeD[0]; yD = shapeD[1]; xD = shapeD[2]; } else if (shapeD.length > 1) { yD = shapeD[0]; xD = shapeD[1]; } // double resolutionX = Math.abs((double) (xAxis.getMaxValue() - xAxis.getMinValue()) / (double) xAxis.getShape()[0]); // double resolutionY = Math.abs((double) (yAxis.getMaxValue() - yAxis.getMinValue()) / (double) yAxis.getShape()[0]); double resolutionX = Math.abs((double) (xAxis.getMaxValue() - xAxis.getMinValue()) / (double) xD); double resolutionY = Math.abs((double) (yAxis.getMaxValue() - yAxis.getMinValue()) / (double) yD); AnalysisLogger.getLogger().debug("Shape: Z:" + zD + " X:" + xD + " Y:" + yD); AnalysisLogger.getLogger().debug("Layer Information Retrieval ELAPSED Time: " + (System.currentTimeMillis() - t01)); int rank = data.getRank(); AnalysisLogger.getLogger().debug("Rank of the layer: " + rank); ArrayFloat.D3 data3Float = null; ArrayDouble.D3 data3Double = null; ArrayInt.D3 data3Int = null; ArrayLong.D3 data3Long = null; ArrayFloat.D2 data2Float = null; ArrayDouble.D2 data2Double = null; ArrayInt.D2 data2Int = null; ArrayLong.D2 data2Long = null; if (data.getRank() == 3) { if (data instanceof ArrayFloat.D3) data3Float = (ArrayFloat.D3) data; else if (data instanceof ArrayInt.D3) data3Int = (ArrayInt.D3) data; else if (data instanceof ArrayDouble.D3) data3Double = (ArrayDouble.D3) data; else if (data instanceof ArrayDouble.D3) data3Double = (ArrayDouble.D3) data; else if (data instanceof ArrayLong.D3) data3Long = (ArrayLong.D3) data; else if (data instanceof ArrayByte.D3) data3Double = (ArrayDouble.D3) VectorOperations.arrayByte3DArrayDouble((ArrayByte) data); else throw new Exception("Layer data format not supported"); } else { if (data instanceof ArrayFloat.D2) data2Float = (ArrayFloat.D2) data; else if (data instanceof ArrayInt.D2) data2Int = (ArrayInt.D2) data; else if (data instanceof ArrayDouble.D2) data2Double = (ArrayDouble.D2) data; else if (data instanceof ArrayLong.D2) data2Long = (ArrayLong.D2) data; else if (data instanceof ArrayByte.D2) data2Double = (ArrayDouble.D2) VectorOperations.arrayByte2DArrayDouble((ArrayByte) data); else throw new Exception("Layer data format not supported"); } double xmin = xAxis.getMinValue(); double xmax = xAxis.getMaxValue(); double ymax = yAxis.getMaxValue(); double ymin = yAxis.getMinValue(); int xmaxidx = (int) Math.round((xmax - xmin) / resolutionX); int ymaxidx = (int) Math.round((ymax - ymin) / resolutionY); boolean is0_360 = false; // if (((xmax == 360) && (xmin == 0)) || ((xmax == 359.5) && (xmin == 0.5))) { // if ((xmin>=0) || (ymin == -77.0104751586914 && ymax==89.94786834716797)) { AnalysisLogger.getLogger().debug("X dimension: " + xD + " Xmin:" + xmin + " Xmax:" + xmax + " Xmaxidx:" + xmaxidx+" XRes: "+resolutionX); AnalysisLogger.getLogger().debug("Y dimension: " + yD + " Ymin:" + ymin + " Ymax:" + ymax + " Ymaxidx:" + ymaxidx+" YRes: "+resolutionY); if ((xmin >= 0)) { xmax = 180; xmin = -180; is0_360 = true; } AnalysisLogger.getLogger().debug("Assigning "+tsize+" grid elements to the NetCDF values"); for (int i = 0; i < tsize; i++) { int zint = 0; int xint = 0; int yint = 0; Tuple triplet = triplets.get(i); double x = triplet.getElements().get(0); double y = triplet.getElements().get(1); if (x == 180) x = -180; if (y == 90) y = -90; double z = 0; if (triplet.getElements().size() > 1) z = triplet.getElements().get(2); if (resolutionZ > 0) { if ((zAxis.getMinValue() <= z) && (zAxis.getMaxValue() >= z)) zint = Math.abs((int) Math.round((z - zAxis.getMinValue()) / resolutionZ)); } if (y < ymin) y = ymin; if (x < xmin) x = xmin; if (y > ymax) y = ymax; if (x > xmax) x = xmax; try{ int[] idxbb = gcs.findXYindexFromLatLon(75,10, null); int[] idxo = gcs.findXYindexFromLatLon(0,0, null); LatLonPoint inverseOrigin = gcs.getLatLon(idxo[0],idxo[1]); // LatLonPoint inverseBB = gcs.getLatLon(idxbb[0],idxbb[1]); //correction to origin offset x = x - inverseOrigin.getLongitude(); y = y - inverseOrigin.getLatitude(); if (i==0) AnalysisLogger.getLogger().debug("bb: " + idxbb[0] +","+idxbb[1]+" origin: "+idxo[0]+","+idxo[1]+" middle "+xD/2+","+yD/2+" shift "+(idxo[0]-(xD/2))+" inverse shift on origin "+inverseOrigin); }catch(Exception e){ AnalysisLogger.getLogger().debug("Error getting x,y corrections "+e.getLocalizedMessage()); e.printStackTrace(); } int[] idx = gcs.findXYindexFromLatLon(y,x, null); xint = idx[0]; yint = idx[1]; if (yint < 0) { yint = 0; } if (xint < 0) { xint = 0; } if (xint > xD - 1) xint = xD - 1; if (yint > yD - 1) yint = yD - 1; /* * if ((xmin <= x) && (xmax >= x)) // xint = (int) Math.round((x - xmin) / resolutionX); { if (is0_360) { if (x < 0) xint = (int) Math.round((x - xmin + xmax) / resolutionX); else xint = (int) Math.round((x) / resolutionX); } else { xint = (int) Math.round((x-xmin) / resolutionX); } } * * if ((yAxis.getMinValue() <= y) && (yAxis.getMaxValue() >= y)) { yint = (int) Math.round((ymax - y) / resolutionY); } * * * if (xint > xD - 1) xint = xD - 1; if (yint > yD - 1) yint = yD - 1; */ Double val = Double.NaN; if (zint > zD - 1) zint = zD - 1; if (data3Float != null) val = Double.valueOf(data3Float.get(zint, yint, xint)); else if (data3Int != null) val = Double.valueOf(data3Int.get(zint, yint, xint)); else if (data3Double != null) val = Double.valueOf(data3Double.get(zint, yint, xint)); else if (data3Long != null) val = Double.valueOf(data3Long.get(zint, yint, xint)); else if (data2Float != null) val = Double.valueOf(data2Float.get(yint, xint)); else if (data2Int != null) val = Double.valueOf(data2Int.get(yint, xint)); else if (data2Double != null) val = Double.valueOf(data2Double.get(yint, xint)); else if (data2Long != null) val = Double.valueOf(data2Long.get(yint, xint)); /*AnalysisLogger.getLogger().debug("Choice "+ (data3Float!=null)+","+ (data3Int!=null)+","+ (data3Double!=null)+","+ (data3Long!=null)+","+ (data2Float!=null)+","+ (data2Int!=null)+","+ (data2Double!=null)+","+ (data2Long!=null)); */ // AnalysisLogger.getLogger().debug("Assigning "+val+" to "+x+","+y+" ["+xint+","+yint+"]"); // AnalysisLogger.getLogger().debug("checking "+data2Float.get(yint, xint)+" vs "); // try{AnalysisLogger.getLogger().debug("checking2 "+data2Float.get(xint,yint));}catch(Exception e){} values.add(val); } break; } } return values; } 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; } public static GridDatatype getGrid(String layer, String netcdffile) throws Exception{ AnalysisLogger.getLogger().debug("Opening File : " + netcdffile); AnalysisLogger.getLogger().debug("Searching for layer: " + layer); GridDataset gds = ucar.nc2.dt.grid.GridDataset.open(netcdffile); List gridTypes = gds.getGrids(); StringBuffer sb = new StringBuffer(); for (GridDatatype gdt : gridTypes) { AnalysisLogger.getLogger().debug("Inside File - layer name: " + gdt.getName()); sb.append(gdt.getName()+" "); if (layer.equals(gdt.getName())) { AnalysisLogger.getLogger().debug("Found layer " + layer + " inside file"); GridDatatype grid = gds.findGridDatatype(gdt.getName()); return grid; } } throw new java.lang.Exception("No layer with name "+layer+" is available in the NetCDF file. Possible values are "+sb.toString()); } // A GridDatatype is like a specialized Variable that explicitly handles X,Y,Z,T dimensions public static LinkedHashMap manageGridDataset(String layer, String filename, double x, double y, double z) throws Exception { LinkedHashMap valuesMap = new LinkedHashMap(); GridDataset gds = ucar.nc2.dt.grid.GridDataset.open(filename); List gridTypes = gds.getGrids(); for (GridDatatype gdt : gridTypes) { AnalysisLogger.getLogger().debug("Inside File - layer name: " + gdt.getName()); if (layer.equalsIgnoreCase(gdt.getName())) { AnalysisLogger.getLogger().debug("Found layer " + layer + " inside file"); GridDatatype grid = gds.findGridDatatype(gdt.getName()); GridCoordSystem gcs = grid.getCoordinateSystem(); long timeSteps = 0; java.util.Date[] dates = null; if (gcs.hasTimeAxis1D()) { CoordinateAxis1DTime tAxis1D = gcs.getTimeAxis1D(); dates = tAxis1D.getTimeDates(); timeSteps = dates.length; } else if (gcs.hasTimeAxis()) { CoordinateAxis tAxis = gcs.getTimeAxis(); timeSteps = tAxis.getSize(); } CoordinateAxis zAxis = gdt.getCoordinateSystem().getVerticalAxis(); double resolutionZ = Math.abs((double) (zAxis.getMaxValue() - zAxis.getMinValue()) / (double) zAxis.getShape()[0]); int zint = 0; if (resolutionZ > 0) { if ((zAxis.getMinValue() <= z) && (zAxis.getMaxValue() >= z)) zint = Math.abs((int) Math.round((z - zAxis.getMinValue()) / resolutionZ)); } AnalysisLogger.getLogger().debug("Z index to take: " + zint); int[] xy = gcs.findXYindexFromLatLon(x, y, null); for (int j = 0; j < timeSteps; j++) { try { Array data = grid.readDataSlice(j, zint, xy[1], xy[0]); // note order is t, z, y, x Double val = takeFirstDouble(data); if (!val.isNaN()) { String date = "" + j; if (dates != null) date = dates[j].toString(); valuesMap.put(timePrefix + date, Double.parseDouble("" + val)); } } catch (Exception e) { AnalysisLogger.getLogger().debug("Error in getting grid values in (" + x + "," + y + "," + z + "= with zint: " + zint + " resolution: " + resolutionZ + " and shape: " + zAxis.getShape()[0]); } } break; } } return valuesMap; } public static Double takeFirstDouble(Array data) { long datal = data.getSize(); Double val = Double.NaN; try { for (int k = 0; k < datal; k++) { Double testVal = data.getDouble(k); if (!testVal.isNaN()) { val = testVal; break; } } } catch (Exception ee) { AnalysisLogger.getLogger().debug("NetCDFDataExplorer-> WARNING: Error in getting value: " + ee.getLocalizedMessage()); } return val; } // A GridDatatype is like a specialized Variable that explicitly handles X,Y,Z,T dimensions public LinkedHashMap managePointsDataset(String layer, String filename, double x, double y) throws Exception { LinkedHashMap valuesMap = new LinkedHashMap(); float tolerance = 0.25f; Formatter errlog = new Formatter(); FeatureDataset fdataset = FeatureDatasetFactoryManager.open(FeatureType.POINT, filename, null, errlog); PointDatasetImpl ds = (PointDatasetImpl) fdataset; List lfc = ds.getPointFeatureCollectionList(); for (FeatureCollection fc : lfc) { StandardPointCollectionImpl spf = (StandardPointCollectionImpl) fc; PointFeatureIterator iter = null; while ((y - tolerance > -90) && (x - tolerance > -180) && (y + tolerance < 90) && (x + tolerance < 180)) { LatLonRect rect = new LatLonRect(new LatLonPointImpl(y - tolerance, x - tolerance), new LatLonPointImpl(y + tolerance, x + tolerance)); PointFeatureCollection coll = spf.subset(rect, null); iter = coll.getPointFeatureIterator(100 * 1000); // 100Kb buffer if (iter.getCount() == 0) iter.finish(); else break; tolerance = tolerance + 0.25f; AnalysisLogger.getLogger().debug("NetCDFDataExplorer-> -> tolerance = " + tolerance); } if (iter != null) { try { while (iter.hasNext()) { ucar.nc2.ft.PointFeature pf = iter.next(); AnalysisLogger.getLogger().debug("NetCDFDataExplorer-> -> EarthLoc: " + pf.getLocation()); AnalysisLogger.getLogger().debug("NetCDFDataExplorer-> -> EarthTime: " + pf.getObservationTime()); StructureData sd = pf.getData(); List mems = sd.getMembers(); for (Member m : mems) { String unit = m.getUnitsString(); if ((unit != null) && (unit.length() > 0)) { AnalysisLogger.getLogger().debug("NetCDFDataExplorer-> -> description: " + m.getDescription()); AnalysisLogger.getLogger().debug("NetCDFDataExplorer-> -> data param: " + m.getDataParam()); AnalysisLogger.getLogger().debug("NetCDFDataExplorer-> -> name: " + m.getName()); AnalysisLogger.getLogger().debug("NetCDFDataExplorer-> -> unit: " + m.getUnitsString()); AnalysisLogger.getLogger().debug("NetCDFDataExplorer-> -> type: " + m.getDataType()); Array arr = sd.getArray(m.getName()); AnalysisLogger.getLogger().debug("NetCDFDataExplorer-> -> is Time: " + m.getDataType()); Double val = takeFirstDouble(arr); AnalysisLogger.getLogger().debug("NetCDFDataExplorer-> -> extracted value: " + val); } } AnalysisLogger.getLogger().debug("NetCDFDataExplorer-> -> EarthTime: "); } } finally { iter.finish(); } } break; } return valuesMap; } // A GridDatatype is like a specialized Variable that explicitly handles X,Y,Z,T dimensions public static boolean isGridDataset(String filename) { try { AnalysisLogger.getLogger().debug("Analyzing file " + filename); Formatter errlog = new Formatter(); FeatureDataset fdataset = FeatureDatasetFactoryManager.open(FeatureType.GRID, filename, null, errlog); if (fdataset == null) { // System.out.printf("GRID Parse failed --> %s\n", errlog); AnalysisLogger.getLogger().debug("NetCDFDataExplorer-> NOT GRID"); return false; } else return true; } catch (Throwable e) { e.printStackTrace(); return false; } } // A GridDatatype is like a specialized Variable that explicitly handles X,Y,Z,T dimensions public static boolean isPointDataset(String filename) { try { Formatter errlog = new Formatter(); FeatureDataset fdataset = FeatureDatasetFactoryManager.open(FeatureType.POINT, filename, null, errlog); if (fdataset == null) { AnalysisLogger.getLogger().debug("NetCDFDataExplorer-> NOT POINT"); return false; } else return true; } catch (Exception e) { return false; } } public static boolean isDataset(String filename) throws Exception { boolean isdataset = false; try { Formatter errlog = new Formatter(); FeatureType[] fts = FeatureType.values(); for (int i = 0; i < fts.length; i++) { FeatureDataset fdataset = FeatureDatasetFactoryManager.open(fts[i], filename, null, errlog); if (fdataset == null) { // System.out.printf(fts[i]+": Parse failed --> %s\n",errlog); } else { AnalysisLogger.getLogger().debug("NetCDFDataExplorer-> " + fts[i] + " OK!"); isdataset = true; } } } catch (Exception e) { } return isdataset; } public static double adjX(double x) { /* * if (x < -180) x = -180; if (x > 180) x = 180; */ return x; } public static double adjY(double y) { /* * if (y < -90) y = -90; if (y > 90) y = 90; */ return y; } public static double getMinX(GridCoordSystem gcs) { CoordinateAxis xAxis = gcs.getXHorizAxis(); return adjX(xAxis.getMinValue()); } public static double getMaxX(GridCoordSystem gcs) { CoordinateAxis xAxis = gcs.getXHorizAxis(); return adjX(xAxis.getMaxValue()); } public static double getMinY(GridCoordSystem gcs) { CoordinateAxis yAxis = gcs.getYHorizAxis(); return adjY(yAxis.getMinValue()); } public static double getMaxY(GridCoordSystem gcs) { CoordinateAxis yAxis = gcs.getYHorizAxis(); return adjY(yAxis.getMaxValue()); } public static double getResolution(String layer, String file) throws Exception{ GridDatatype gdt = getGrid(layer, file); double minY = NetCDFDataExplorer.getMinY(gdt.getCoordinateSystem()); double maxY = NetCDFDataExplorer.getMaxY(gdt.getCoordinateSystem()); GridDatatype gdtsub = gdt.makeSubset(new Range(0, 0), null, null, 1, 1, 1); Array data = gdtsub.readVolumeData(0); // note order is t, z, y, x int[] shapeD = data.getShape(); int yD = 0; if (shapeD.length > 2) yD = shapeD[1]; else yD = shapeD[0]; double resolutionY = Math.abs((double) (maxY - minY) / (double) yD); resolutionY = MathFunctions.roundDecimal(resolutionY, 4); return resolutionY; } }