From 48319abfee2c59c2e8865531b142e9f90bfbbedf Mon Sep 17 00:00:00 2001
From: "benjamin.jakimow" <benjamin.jakimow@geo.hu-berlin.de>
Date: Fri, 14 Jun 2019 09:00:10 +0200
Subject: [PATCH] fixes issues related to netCDF data

Signed-off-by: benjamin.jakimow <benjamin.jakimow@geo.hu-berlin.de>
---
 eotimeseriesviewer/pixelloader.py | 12 ++++--------
 eotimeseriesviewer/timeseries.py  | 19 +++++++++++++------
 tests/test_main.py                |  1 -
 tests/test_pixelloader.py         | 11 +++++++++++
 tests/test_timeseries.py          |  4 ++++
 5 files changed, 32 insertions(+), 15 deletions(-)

diff --git a/eotimeseriesviewer/pixelloader.py b/eotimeseriesviewer/pixelloader.py
index b378c29a..4eaaea80 100644
--- a/eotimeseriesviewer/pixelloader.py
+++ b/eotimeseriesviewer/pixelloader.py
@@ -186,9 +186,6 @@ def transformPoint2Px(trans, pt, gt):
     return geo2px(QgsPointXY(x, y), gt)
 
 
-
-
-
 def doLoaderTask(taskWrapper:QgsTask, dump):
 
     #assert isinstance(taskWrapper, QgsTask)
@@ -210,14 +207,13 @@ def doLoaderTask(taskWrapper:QgsTask, dump):
     gt = ds.GetGeoTransform()
     result.resGeoTransformation = gt
     result.resCrsWkt = ds.GetProjection()
+    if result.resCrsWkt == '':
+        result.resCrsWkt = QgsRasterLayer(ds.GetDescription()).crs().toWkt()
     crsSrc = osr.SpatialReference(result.resCrsWkt)
 
-
-    #convert Geometries into pixel indices to be extracted
+    # convert Geometries into pixel indices to be extracted
     PX_SUBSETS = []
 
-
-
     for geom in task.geometries:
         crsRequest = osr.SpatialReference()
 
@@ -332,7 +328,7 @@ def doLoaderTask(taskWrapper:QgsTask, dump):
 
 
 
-    #finally, ensure that there is on 2D array only
+    # finally, ensure that there is on 2D array only
     for i in range(len(PROFILE_DATA)):
         d = PROFILE_DATA[i]
         if isinstance(d, np.ndarray):
diff --git a/eotimeseriesviewer/timeseries.py b/eotimeseriesviewer/timeseries.py
index b9698f15..30921982 100644
--- a/eotimeseriesviewer/timeseries.py
+++ b/eotimeseriesviewer/timeseries.py
@@ -400,24 +400,31 @@ class TimeSeriesSource(object):
             assert self.mDate is not None, 'Unable to find acquisition date of {}'.format(self.mUri)
 
             self.mDrv = dataset.GetDriver().ShortName
-            self.mGT = dataset.GetGeoTransform()
-            self.mWKT = dataset.GetProjection()
-            self.mCRS = QgsCoordinateReferenceSystem(self.mWKT)
 
             self.mWL, self.mWLU = extractWavelengths(dataset)
 
 
             self.nb, self.nl, self.ns = dataset.RasterCount, dataset.RasterYSize, dataset.RasterXSize
             self.mGeoTransform = dataset.GetGeoTransform()
+            self.mMetaData = collections.OrderedDict()
+            for domain in dataset.GetMetadataDomainList():
+                self.mMetaData[domain] = dataset.GetMetadata_Dict(domain)
+
+            self.mWKT = dataset.GetProjection()
+            if self.mWKT == '':
+                # no CRS? try with QGIS API
+                lyr = QgsRasterLayer(self.mUri)
+                if lyr.crs().isValid():
+                    self.mWKT = lyr.crs().toWkt()
+
+            self.mCRS = QgsCoordinateReferenceSystem(self.mWKT)
+
             px_x = float(abs(self.mGeoTransform[1]))
             px_y = float(abs(self.mGeoTransform[5]))
             self.mGSD = (px_x, px_y)
             self.mDataType = dataset.GetRasterBand(1).DataType
             self.mSid = sensorID(self.nb, px_x, px_y, self.mDataType, self.mWL, self.mWLU)
 
-            self.mMetaData = collections.OrderedDict()
-            for domain in dataset.GetMetadataDomainList():
-                self.mMetaData[domain] = dataset.GetMetadata_Dict(domain)
 
             self.mUL = QgsPointXY(*px2geo(QPoint(0, 0), self.mGeoTransform, pxCenter=False))
             self.mLR = QgsPointXY(*px2geo(QPoint(self.ns + 1, self.nl + 1), self.mGeoTransform, pxCenter=False))
diff --git a/tests/test_main.py b/tests/test_main.py
index e7640e42..57bce4c2 100644
--- a/tests/test_main.py
+++ b/tests/test_main.py
@@ -149,7 +149,6 @@ class TestInit(TestCase):
         if SHOW_GUI:
             QGIS_APP.exec_()
 
-
     def test_TimeSeriesViewerNoSource(self):
 
         from eotimeseriesviewer.main import TimeSeriesViewer
diff --git a/tests/test_pixelloader.py b/tests/test_pixelloader.py
index 76ec79a0..8ec77242 100644
--- a/tests/test_pixelloader.py
+++ b/tests/test_pixelloader.py
@@ -142,6 +142,17 @@ class PixelLoaderTest(unittest.TestCase):
 
         self.assertTrue(len(RESULTS) == len(tasks))
 
+    def test_pixelLoader_netCDF(self):
+
+        p = r'Q:\Processing_BJ\99_OSARIS_Testdata\Loibl-2019-OSARIS-Ala-Archa\Amplitudes\20151207--20151231-amplitude.grd'
+        if os.path.isfile(p):
+            lyr = QgsRasterLayer(p)
+            task = PixelLoaderTask(p, [SpatialPoint.fromMapLayerCenter(lyr)])
+
+            result = doLoaderTask(None, task.toDump())
+
+            s = ""
+
     def test_pixelLoader(self):
         from eotimeseriesviewer.pixelloader import doLoaderTask, PixelLoaderTask, INFO_OUT_OF_IMAGE, INFO_NO_DATA
         from eotimeseriesviewer import px2geo
diff --git a/tests/test_timeseries.py b/tests/test_timeseries.py
index 4e1ec1a5..ff2e9b0e 100644
--- a/tests/test_timeseries.py
+++ b/tests/test_timeseries.py
@@ -148,6 +148,10 @@ class TestInit(unittest.TestCase):
             tss = TimeSeriesSource.create(p)
             self.assertIsInstance(tss, TimeSeriesSource)
 
+        p = r'Q:\Processing_BJ\99_OSARIS_Testdata\Loibl-2019-OSARIS-Ala-Archa\Amplitudes\20151207--20151231-amplitude.grd'
+        if os.path.isfile(p):
+            tss = TimeSeriesSource.create(p)
+            self.assertIsInstance(tss, TimeSeriesSource)
 
         if False:
             webSources = [QgsRasterLayer(wcs, 'test', 'wcs')]
-- 
GitLab