471 lines
18 KiB
Java
471 lines
18 KiB
Java
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.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.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<Double> retrieveDataFromNetCDF(String openDapLink, String layer, int time, List<Tuple<Double>> triplets, double xL, double xR, double yL, double yR) {
|
|
try {
|
|
List<Double> values = new ArrayList<Double>();
|
|
if (isGridDataset(openDapLink)) {
|
|
AnalysisLogger.getLogger().debug("Managing Grid File");
|
|
return manageGridDataset(layer, openDapLink, time, triplets, xL, xR, yL, yR);
|
|
}
|
|
/*
|
|
* 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<GridDatatype> 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){
|
|
AnalysisLogger.getLogger().debug("NetCDF Explorer Error:"+e.getLocalizedMessage());
|
|
}
|
|
}
|
|
|
|
// A GridDatatype is like a specialized Variable that explicitly handles X,Y,Z,T dimensions
|
|
GridDataset gds;
|
|
public List<Double> manageGridDataset(String layer, String filename, int time, List<Tuple<Double>> triplets, double xL, double xR, double yL, double yR) throws Exception {
|
|
List<Double> values = new ArrayList<Double>();
|
|
if (gds==null)
|
|
gds = ucar.nc2.dt.grid.GridDataset.open(filename);
|
|
|
|
List<GridDatatype> 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());
|
|
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){};
|
|
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]);
|
|
|
|
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];
|
|
}
|
|
|
|
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();
|
|
if (((xmax==360) && (xmin==0)) || ((xmax==359.5) && (xmin==0.5))){
|
|
xmax = 180;
|
|
xmin=-180;
|
|
}
|
|
AnalysisLogger.getLogger().debug("X dimension: "+xD+" Xmin:"+ xmax+" Xmax:"+xmin);
|
|
|
|
for (int i = 0; i < tsize; i++) {
|
|
int zint = 0;
|
|
int xint = 0;
|
|
int yint = 0;
|
|
Tuple<Double> 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));
|
|
}
|
|
|
|
// AnalysisLogger.getLogger().debug("Z Index: "+zint);
|
|
/*
|
|
GridCoordSystem gcs = grid.getCoordinateSystem();
|
|
int[] xy = gcs.findXYindexFromLatLon(x, y, null);
|
|
Array datas=grid.readDataSlice(time, zint, xy[1], xy[0]);
|
|
*/
|
|
if ((xmin <= x) && (xmax >= x))
|
|
xint = (int) Math.round((x - xmin) / resolutionX);
|
|
if ((yAxis.getMinValue() <= y) && (yAxis.getMaxValue() >= y))
|
|
yint = (int) Math.round((y - yAxis.getMinValue()) / resolutionY);
|
|
Double val = Double.NaN;
|
|
if (xint > xD-1)
|
|
xint = xD-1;
|
|
if (yint > yD-1)
|
|
yint = yD-1;
|
|
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));
|
|
|
|
values.add(val);
|
|
}
|
|
break;
|
|
}
|
|
}
|
|
return values;
|
|
}
|
|
|
|
// A GridDatatype is like a specialized Variable that explicitly handles X,Y,Z,T dimensions
|
|
public static LinkedHashMap<String, Double> manageGridDataset(String layer, String filename, double x, double y, double z) throws Exception {
|
|
LinkedHashMap<String, Double> valuesMap = new LinkedHashMap<String, Double>();
|
|
GridDataset gds = ucar.nc2.dt.grid.GridDataset.open(filename);
|
|
List<GridDatatype> 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<String, String> managePointsDataset(String layer, String filename, double x, double y) throws Exception {
|
|
LinkedHashMap<String, String> valuesMap = new LinkedHashMap<String, String>();
|
|
float tolerance = 0.25f;
|
|
Formatter errlog = new Formatter();
|
|
FeatureDataset fdataset = FeatureDatasetFactoryManager.open(FeatureType.POINT, filename, null, errlog);
|
|
PointDatasetImpl ds = (PointDatasetImpl) fdataset;
|
|
List<FeatureCollection> 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<Member> 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());
|
|
}
|
|
}
|