2014-02-24 12:52:57 +01:00
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 ;
2014-02-26 16:55:31 +01:00
import org.gcube.dataanalysis.geo.utils.VectorOperations ;
2014-02-24 12:52:57 +01:00
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: " ;
2014-02-26 16:55:31 +01:00
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 ) {
2014-02-24 12:52:57 +01:00
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 ) ;
2014-02-26 01:02:21 +01:00
// e.printStackTrace();
2014-02-24 12:52:57 +01:00
return null ;
}
}
2014-02-26 16:55:31 +01:00
public double minZ = 0 ;
public double maxZ = 0 ;
private 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 ( ) ) ;
if ( layer . equalsIgnoreCase ( gdt . getFullName ( ) ) ) {
CoordinateAxis zAxis = gdt . getCoordinateSystem ( ) . getVerticalAxis ( ) ;
minZ = zAxis . getMinValue ( ) ;
maxZ = zAxis . getMaxValue ( ) ;
break ;
}
}
}
} catch ( Exception e ) {
AnalysisLogger . getLogger ( ) . debug ( " NetCDF Explorer Error: " + e . getLocalizedMessage ( ) ) ;
}
}
2014-02-24 12:52:57 +01:00
// A GridDatatype is like a specialized Variable that explicitly handles X,Y,Z,T dimensions
2014-02-26 16:55:31 +01:00
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 {
2014-02-24 12:52:57 +01:00
List < Double > values = new ArrayList < Double > ( ) ;
2014-02-26 16:55:31 +01:00
if ( gds = = null )
gds = ucar . nc2 . dt . grid . GridDataset . open ( filename ) ;
2014-02-24 12:52:57 +01:00
List < GridDatatype > gridTypes = gds . getGrids ( ) ;
for ( GridDatatype gdt : gridTypes ) {
AnalysisLogger . getLogger ( ) . debug ( " Inside File - layer name: " + gdt . getFullName ( ) ) ;
if ( layer . equalsIgnoreCase ( gdt . getFullName ( ) ) ) {
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 )
2014-02-26 16:55:31 +01:00
data3Double = ( ArrayDouble . D3 ) VectorOperations . arrayByte3DArrayDouble ( ( ArrayByte ) data ) ;
2014-02-24 12:52:57 +01:00
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 )
2014-02-26 16:55:31 +01:00
data2Double = ( ArrayDouble . D2 ) VectorOperations . arrayByte2DArrayDouble ( ( ArrayByte ) data ) ;
2014-02-24 12:52:57 +01:00
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 . getFullName ( ) ) ;
if ( layer . equalsIgnoreCase ( gdt . getFullName ( ) ) ) {
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 ) {
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 ( ) ) ;
}
}