# 366 #369Merge branches 'develop' and 'issue366_arome_mapper'

into issue366_arome_mapper

nansat/domain.py

@@ -17,16 +17,16 @@
 from __future__ import division, absolute_import
 
 import re
-import numpy as np
+import warnings
 from xml.etree.ElementTree import ElementTree
 
+import numpy as np
+
 from nansat.tools import add_logger, initial_bearing, haversine, gdal, osr, ogr
 from nansat.nsr import NSR
 from nansat.vrt import VRT
-
-import warnings
 from nansat.exceptions import NansatProjectionError
-
+from nansat.warnings import NansatFutureWarning
 
 class Domain(object):
     """Container for geographical reference of a raster
@@ -55,23 +55,13 @@ class Domain(object):
         -te xmin ymin xmax ymax
         -lle lonmin latmin lonmax latmax
     ds : GDAL dataset
-    lat : Numpy array
-        Grid with latitudes
-    lon : Numpy array
-        Grid with longitudes
-    name : string, optional
-        Name to be added to the Domain object
-    logLevel : int, optional
-        level of logging
 
     Examples
     --------
         >>> d = Domain(srs, ext) #size, extent and spatial reference is given by strings
         >>> d = Domain(ds=GDALDataset) #size, extent copied from input GDAL dataset
         >>> d = Domain(srs, ds=GDALDataset) # spatial reference is given by srs,
             but size and extent is determined from input GDAL dataset
-        >>> d = Domain(lon=lonGrid, lat=latGrid) # Size, extent and spatial reference is given
-            by two grids of longitude and latitude
 
     Notes
     -----
@@ -134,17 +124,8 @@ class Domain(object):
     logger = None
     name = None
 
-    def __init__(self, srs=None, ext=None, ds=None, lon=None,
-                 lat=None, name='', logLevel=None):
+    def __init__(self, srs=None, ext=None, ds=None, **kwargs):
         """Create Domain from GDALDataset or string options or lat/lon grids"""
-        # set default attributes
-        self.logger = add_logger('Nansat', logLevel)
-        self.name = name
-
-        self.logger.debug('ds: %s' % str(ds))
-        self.logger.debug('srs: %s' % srs)
-        self.logger.debug('ext: %s' % ext)
-
         # If too much information is given raise error
         if ds is not None and srs is not None and ext is not None:
             raise ValueError('Ambiguous specification of both dataset, srs- and ext-strings.')
@@ -153,7 +134,6 @@ def __init__(self, srs=None, ext=None, ds=None, lon=None,
         # ds
         # ds and srs
         # srs and ext
-        # lon and lat
 
         # if only a dataset is given:
         #     copy geo-reference from the dataset
@@ -187,14 +167,42 @@ def __init__(self, srs=None, ext=None, ds=None, lon=None,
                                                geo_transform=geo_transform,
                                                projection=srs.wkt,
                                                gcps=[], gcp_projection='')
-        elif lat is not None and lon is not None:
+        elif 'lat' in kwargs and 'lon' in kwargs:
+            warnings.warn('Domain(lon=lon, lat=lat) will be deprectaed!'
+                          'Use Domain.from_lonlat()', NansatFutureWarning)
             # create self.vrt from given lat/lon
-            self.vrt = VRT.from_lonlat(lon, lat)
+            self.vrt = VRT.from_lonlat(kwargs['lon'], kwargs['lat'])
         else:
             raise ValueError('"dataset" or "srsString and extentString" '
                               'or "dataset and srsString" are required')
 
-        self.logger.debug('vrt.dataset: %s' % str(self.vrt.dataset))
+    @classmethod
+    def from_lonlat(cls, lon, lat, add_gcps=True):
+        """Create Domain object from input longitudes, latitudes arrays
+
+        Parameters
+        ----------
+        lon : numpy.ndarray
+            longitudes
+        lat : numpy.ndarray
+            latitudes
+        add_gcps : bool
+            Add GCPs from lon/lat arrays.
+
+        Returns
+        -------
+            d : Domain
+
+        Examples
+        --------
+            >>> lon, lat = np.meshgrid(range(10), range(10))
+            >>> d1 = Domain.from_lonlat(lon, lat)
+            >>> d2 = Domain.from_lonlat(lon, lat, add_gcps=False) # add only geolocation arrays
+
+        """
+        d = cls.__new__(cls)
+        d.vrt = VRT.from_lonlat(lon, lat, add_gcps)
+        return d
 
     def __repr__(self):
         """Creates string with basic info about the Domain object
@@ -206,12 +214,12 @@ def __repr__(self):
 
         """
         corners_temp = '\t (%6.2f, %6.2f)  (%6.2f, %6.2f)\n'
-
+        wkt, src = self.vrt.get_projection()
         out_str = 'Domain:[%d x %d]\n' % self.shape()[::-1]
         out_str += self.OUTPUT_SEPARATOR
         corners = self.get_corners()
-        out_str += 'Projection:\n'
-        out_str += (NSR(self.vrt.get_projection()).ExportToPrettyWkt(1) + '\n')
+        out_str += 'Projection(%s):\n' % src
+        out_str += (NSR(wkt).ExportToPrettyWkt(1) + '\n')
         out_str += self.OUTPUT_SEPARATOR
         out_str += 'Corners (lon, lat):\n'
         out_str += corners_temp % (corners[0][0], corners[1][0], corners[0][2], corners[1][2])
@@ -378,7 +386,7 @@ def get_geolocation_grids(self, stepSize=1, dstSRS=NSR()):
         y_vec = list(range(0, self.vrt.dataset.RasterYSize, step_size))
         x_grid, y_grid = np.meshgrid(x_vec, y_vec)
 
-        if hasattr(self.vrt, 'geolocation') and len(self.vrt.geolocation.data) > 0:
+        if self.vrt.geolocation is not None and len(self.vrt.geolocation.data) > 0:
             # if the vrt dataset has geolocationArray
             # read lon,lat grids from geolocationArray
             lon_grid, lat_grid = self.vrt.geolocation.get_geolocation_grids()

nansat/nansat.py

@@ -250,7 +250,7 @@ def _init_from_domain(self, domain, array=None, parameters=None, log_level=30):
 
         """
         self._init_empty('', log_level)
-        self.vrt = VRT.from_gdal_dataset(domain.vrt.dataset)
+        self.vrt = VRT.from_gdal_dataset(domain.vrt.dataset, geolocation=domain.vrt.geolocation)
         self.mapper = ''
         if array is not None:
             self.add_band(array=array, parameters=parameters)

nansat/nsr.py

@@ -13,7 +13,7 @@
 # but WITHOUT ANY WARRANTY without even the implied warranty of
 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
 from __future__ import absolute_import, unicode_literals
-
+import sys
 import osr
 
 from nansat.exceptions import NansatProjectionError
@@ -56,6 +56,10 @@ class NSR(osr.SpatialReference, object):
     """
     def __init__(self, srs=0):
         """Create Spatial Reference System from input parameter"""
+        if sys.version_info.major == 2:
+            str_types = (str, unicode)
+        else:
+            str_types = str
         # create SRS
         osr.SpatialReference.__init__(self)
 
@@ -64,7 +68,7 @@ def __init__(self, srs=0):
         if srs is 0:
             # generate default WGS84 SRS
             status = self.ImportFromWkt(osr.SRS_WKT_WGS84)
-        elif isinstance(srs, str):
+        elif isinstance(srs, str_types):
             # parse as proj4 string
             status = self.ImportFromProj4(str(srs))
             if status > 0:

nansat/tests/test_domain.py

@@ -10,19 +10,24 @@
 #               under the terms of GNU General Public License, v.3
 #               http://www.gnu.org/licenses/gpl-3.0.html
 #------------------------------------------------------------------------------
-import unittest
 import os
+import sys
+import warnings
+import unittest
+
 import numpy as np
+from mock import patch, PropertyMock, Mock, MagicMock, DEFAULT
 
 from nansat.nsr import NSR
 from nansat.vrt import VRT
 from nansat.domain import Domain
 from nansat.tools import gdal, ogr
-import sys
 from nansat.tests import nansat_test_data as ntd
-from mock import patch, PropertyMock, Mock, MagicMock, DEFAULT
-
 from nansat.exceptions import NansatProjectionError
+from nansat.warnings import NansatFutureWarning
+
+warnings.simplefilter("always", NansatFutureWarning)
+warnings.simplefilter("always", UserWarning)
 
 EXTENT_TE_TS = "-te 25 70 35 72 -ts 500 500"
 EXTENT_DICT_TE_TS = {'te': [25.0, 70.0, 35.0, 72.0], 'ts': [500.0, 500.0]}
@@ -94,20 +99,42 @@ def test_init_from_srs_and_ext_lle(self, mock__get_geotransform, mock__convert_e
                    ext="-lle 25 70 35 72 -ts 500 500")
         self.assertEqual(type(d), Domain)
 
-    def test_init_from_lonlat(self):
+    def test_init_lonlat(self):
         lat, lon = np.mgrid[-90:90:0.5, -180:180:0.5]
-        d = Domain(lon=lon, lat=lat)
+        with warnings.catch_warnings(record=True) as recorded_warnings:
+            d = Domain(lon=lon, lat=lat)
+
+        nansat_warning_raised = False
+        for rw in recorded_warnings:
+            if rw.category == NansatFutureWarning:
+                nansat_warning_raised = True
+        self.assertTrue(nansat_warning_raised)
+
         self.assertEqual(type(d), Domain)
         self.assertEqual(d.shape(), lat.shape)
 
+    def test_init_from_lonlat(self):
+        lat, lon = np.mgrid[-10:10:0.5, -20:20:2]
+        d = Domain.from_lonlat(lon=lon, lat=lat)
+        self.assertEqual(type(d), Domain)
+        self.assertEqual(d.shape(), lat.shape)
+        self.assertEqual(len(d.vrt.dataset.GetGCPs()), 100)
+
+    def test_init_from_lonlat_no_gcps(self):
+        lat, lon = np.mgrid[-10:10:0.5, -20:20:2]
+        d = Domain.from_lonlat(lon=lon, lat=lat, add_gcps=False)
+        self.assertEqual(type(d), Domain)
+        self.assertEqual(d.shape(), lat.shape)
+        self.assertEqual(len(d.vrt.dataset.GetGCPs()), 0)
+
     @patch.object(Domain, 'get_corners',
         return_value=(np.array([ 25.,  25.,  35.,  35.]), np.array([ 72.,  70.,  72.,  70.])))
     def test_repr(self, mock_get_corners):
         d = Domain(4326, "-te 25 70 35 72 -ts 500 500")
         result = d.__repr__()
         test = ('Domain:[500 x 500]\n'
                 '----------------------------------------\n'
-                'Projection:\nGEOGCS["WGS 84",\n'
+                'Projection(dataset):\nGEOGCS["WGS 84",\n'
                 '    DATUM["WGS_1984",\n'
                 '        SPHEROID["WGS 84",6378137,298.257223563]],\n'
                 '    PRIMEM["Greenwich",0],\n'
@@ -367,12 +394,6 @@ def test_contains(self):
         self.assertFalse(Norway.contains(WestCoast))
         self.assertFalse(Paris.contains(Norway))
 
-    def test_transform_ts(self):
-        result = Domain._transform_ts(1.5, 1.0, [750.0, 500.0])
-        self.assertIsInstance(result, tuple)
-        self.assertEqual(len(result), 4)
-        for el in result:
-            self.assertIsInstance(el, float)
     def test_get_border_postgis(self):
         d = Domain(4326, "-te 25 70 35 72 -ts 500 500")
         result = d.get_border_postgis()
@@ -412,6 +433,7 @@ def test_get_pixelsize_meters(self):
         x, y = d.get_pixelsize_meters()
         self.assertEqual(int(x), 500)
         self.assertEqual(int(y), 500)
+
     def test_get_geotransform(self):
         input_1 = {'te': [25.0, 70.0, 35.0, 72.0], 'ts': [500.0, 500.0]}
         test_1 = ([25.0, 0.02, 0.0, 72.0, 0.0, -0.004], 500, 500)
@@ -436,7 +458,7 @@ def test_transform_tr(self):
             self.assertEqual(param_err.message,
                              '"-tr" is too large. width is 4.0, height is 1.3 ')
 
-    def test_transform_ts(self):
+    def test_transform_ts2(self):
         result = Domain._transform_ts(1.5, 1.0, [750.0, 500.0])
         self.assertIsInstance(result, tuple)
         self.assertEquals(len(result), 4)
@@ -490,18 +512,6 @@ def test_reproject_gcps_auto(self):
         self.assertEqual(gcpproj.GetAttrValue('PROJECTION'),
                         'Stereographic')
 
-    def test_repr(self):
-        dom = Domain(4326, "-te 4.5 60 6 61 -ts 750 500")
-        result = dom.__repr__()
-        test = 'Domain:[750 x 500]\n----------------------------------------\nProjection:\nGEOGC' \
-               'S["WGS 84",\n    DATUM["WGS_1984",\n        SPHEROID["WGS 84",6378137,298.257223' \
-               '563]],\n    PRIMEM["Greenwich",0],\n    UNIT["degree",0.0174532925199433]]\n-----' \
-               '-----------------------------------\nCorners (lon, lat):\n\t (  4.50,  61.00)  ' \
-               '(  6.00,  61.00)\n\t (  4.50,  60.00)  (  6.00,  60.00)\n'
-
-        self.assertIsInstance(result, str)
-        self.assertEqual(result, test)
-
     @patch.multiple(Domain, get_border_geometry=DEFAULT, __init__ = Mock(return_value=None))
     def test_intersects(self, get_border_geometry):
         other_domain = MagicMock()

nansat/tests/test_nansat.py

@@ -75,6 +75,14 @@ def test_from_domain_array(self):
         self.assertEqual(n.name, '')
         self.assertEqual(n.path, '')
 
+    def test_from_domain_nansat(self):
+        n1 = Nansat(self.test_file_gcps, log_level=40, mapper=self.default_mapper)
+        n2 = Nansat.from_domain(n1, n1[1])
+
+        self.assertEqual(type(n2), Nansat)
+        self.assertEqual(len(n2.bands()), 1)
+        self.assertEqual(type(n2[1]), np.ndarray)
+
     def test_add_band(self):
         d = Domain(4326, "-te 25 70 35 72 -ts 500 500")
         arr = np.random.randn(500, 500)
@@ -742,7 +750,7 @@ def test_repr_basic(self):
         self.assertIn('SourceFilename', n_repr)
         self.assertIn('/vsimem/', n_repr)
         self.assertIn('500 x 500', n_repr)
-        self.assertIn('Projection:', n_repr)
+        self.assertIn('Projection(dataset):', n_repr)
         self.assertIn('25', n_repr)
         self.assertIn('72', n_repr)
         self.assertIn('35', n_repr)
@@ -787,6 +795,23 @@ def test_get_metadata_unescape(self, vrt):
         self.assertEqual(meta1, {'key1': '" AAA " & > <', 'key2': "'BBB'"})
         self.assertEqual(meta2, meta0)
 
+    def test_reproject_pure_geolocation(self):
+        n0 = Nansat(self.test_file_gcps)
+        b0 = n0[1]
+        lon0, lat0 = n0.get_geolocation_grids()
+        d1 = Domain.from_lonlat(lon=lon0, lat=lat0)
+        d2 = Domain.from_lonlat(lon=lon0, lat=lat0, add_gcps=False)
+        d3 = Domain(NSR().wkt, '-te 27 70 31 72 -ts 500 500')
+
+        n1 = Nansat.from_domain(d1, b0)
+        n2 = Nansat.from_domain(d2, b0)
+
+        n1.reproject(d3)
+        n2.reproject(d3)
+
+        b1 = n1[1]
+        b2 = n2[1]
+        self.assertTrue(np.allclose(b1,b2))
 
 if __name__ == "__main__":
     unittest.main()

nansat/tests/test_nsr.py

@@ -2,7 +2,7 @@
 # Name:         test_nansat.py
 # Purpose:      Test the nansat module
 #
-# Author:       Morten Wergeland Hansen, Asuka Yamakawa
+# Author:       Morten Wergeland Hansen, Anton Korosov, Asuka Yamakawa
 # Modified: Morten Wergeland Hansen
 #
 # Created:  18.06.2014
@@ -50,6 +50,13 @@ def test_init_from_proj4(self):
         self.assertEqual(nsr.Validate(), 0)
         self.assertTrue('longlat' in nsr.ExportToProj4())
 
+    def test_init_from_proj4_unicode(self):
+        nsr = NSR(u'+proj=longlat')
+
+        self.assertEqual(type(nsr), NSR)
+        self.assertEqual(nsr.Validate(), 0)
+        self.assertTrue('longlat' in nsr.ExportToProj4())
+
     def test_init_from_wkt(self):
         nsr = NSR(osr.SRS_WKT_WGS84)