Merge fd4f4e1 into 72751e8

.gitignore

@@ -19,3 +19,5 @@ docs/_build/
 .idea
 .settings
 .vagrant
+provisioning/roles/andrewrothstein*
+provisioning/site.retry

nansat/mappers/mapper_ascat.py

@@ -0,0 +1,65 @@
+import re
+import json
+import numpy as np
+from datetime import datetime
+import pythesint as pti
+
+from netCDF4 import Dataset
+
+import gdal
+
+from nansat.vrt import VRT
+from nansat.geolocation import Geolocation
+from nansat.nsr import NSR
+from nansat.domain import Domain
+from nansat.mappers.scatterometers import Mapper as ScatterometryMapper
+from nansat.exceptions import WrongMapperError
+
+class Mapper(ScatterometryMapper):
+    """ Nansat mapper for ASCAT """
+
+    def __init__(self, filename, gdal_dataset, metadata, quartile=0, *args, **kwargs):
+
+        if not 'ascat' in metadata.get('NC_GLOBAL#source', '').lower():
+            raise WrongMapperError
+
+        super(Mapper, self).__init__(filename, gdal_dataset, metadata, quartile=quartile, *args, **kwargs)
+
+        lat = self.dataset.GetRasterBand(self._latitude_band_number(gdal_dataset)).ReadAsArray()
+        lon = self.dataset.GetRasterBand(self._longitude_band_number(gdal_dataset)).ReadAsArray()
+        lon = ScatterometryMapper.shift_longitudes(lon)
+        self.set_gcps(lon, lat, gdal_dataset)
+
+        # Get dictionary describing the instrument and platform according to
+        # the GCMD keywords
+        ii = pti.get_gcmd_instrument('ascat')
+        pp = pti.get_gcmd_platform(metadata['NC_GLOBAL#source'].split(' ')[0])
+        provider = pti.get_gcmd_provider(re.split('[^a-zA-Z]',
+            metadata['NC_GLOBAL#institution'])[0])
+
+        # TODO: Validate that the found instrument and platform are indeed what
+        # we want....
+
+        self.dataset.SetMetadataItem('instrument', json.dumps(ii))
+        self.dataset.SetMetadataItem('platform', json.dumps(pp))
+        self.dataset.SetMetadataItem('data_center', json.dumps(provider))
+        self.dataset.SetMetadataItem('entry_title', metadata['NC_GLOBAL#title'])
+        self.dataset.SetMetadataItem('ISO_topic_category',
+                json.dumps(pti.get_iso19115_topic_category('Oceans')))
+
+    def times(self):
+        """ Get times from time variable
+        """
+        ds = Dataset(self.input_filename)
+
+        # Get datetime object of epoch and time_units string
+        time_units = self._time_units(ds=ds)
+
+        # Get all times - slight difference from NetCDF-CF mappers times method...
+        times = ds.variables[self._timevarname(ds=ds)][:,0]
+
+        # Create numpy array of np.datetime64 times (provide epoch to save time)
+        tt = np.array([self._time_count_to_np_datetime64(tn,
+            time_units=time_units) for tn in times])
+
+        return tt

nansat/mappers/mapper_netcdf_cf.py

@@ -142,7 +142,7 @@ def _time_count_to_np_datetime64(self, time_count, time_units=None):
             raise Exception('Check time units..')
         return tt
 
-    def _band_list(self, gdal_dataset, gdal_metadata, netcdf_dim={}, bands=[]):
+    def _band_list(self, gdal_dataset, gdal_metadata, netcdf_dim={}, bands=[], *args, **kwargs):
         ''' Create list of dictionaries mapping source and destination metadata
         of bands that should be added to the Nansat object.
 
@@ -250,8 +250,8 @@ def _band_dict(self, subfilename, band_num, subds, band=None, band_metadata=None
             try:
                 band = subds.GetRasterBand(band_num)
             except RuntimeError as e:
-                if 'illegal band' in e.message.lower():
-                    warnings.warn('Skipping band due to GDAL error: %s' %e.message)
+                if 'illegal band' in str(e).lower():
+                    warnings.warn('Skipping band due to GDAL error: %s' %str(e))
                     return {}
                 else:
                     raise
@@ -272,7 +272,7 @@ def _band_dict(self, subfilename, band_num, subds, band=None, band_metadata=None
                 # Then we don't need time for this band...
                 warnings.warn(
                         '%s: %s - %s Continuing without time metadata for band %s'
-                        %(e.__repr__().split('(')[0], e.message, e.__doc__,
+                        %(e.__repr__().split('(')[0], str(e), e.__doc__,
                             band_metadata['NETCDF_VARNAME']))
 
         # Generate source metadata

nansat/mappers/mapper_quikscat.py

@@ -0,0 +1,47 @@
+import re
+import json
+import numpy as np
+from datetime import datetime
+import pythesint as pti
+
+import gdal
+
+from nansat.vrt import VRT
+from nansat.geolocation import Geolocation
+from nansat.nsr import NSR
+from nansat.domain import Domain
+#from nansat.mappers.mapper_netcdf_cf import Mapper as NetcdfCF
+from nansat.mappers.scatterometers import Mapper as ScatterometryMapper
+from nansat.exceptions import WrongMapperError
+
+#class Mapper(NetcdfCF):
+class Mapper(ScatterometryMapper):
+    """ Nansat mapper for QuikScat """
+
+    def __init__(self, filename, gdal_dataset, metadata, quartile=0, *args, **kwargs):
+
+        if not 'quikscat' in metadata.get('NC_GLOBAL#source', '').lower():
+            raise WrongMapperError
+
+        super(Mapper, self).__init__(filename, gdal_dataset, metadata, quartile=quartile, *args, **kwargs)
+
+        lat = self.dataset.GetRasterBand(self._latitude_band_number(gdal_dataset)).ReadAsArray()
+        lon = self.dataset.GetRasterBand(self._longitude_band_number(gdal_dataset)).ReadAsArray()
+        lon = ScatterometryMapper.shift_longitudes(lon)
+        self.set_gcps(lon, lat, gdal_dataset)
+
+        # Get dictionary describing the instrument and platform according to
+        # the GCMD keywords
+        mm = pti.get_gcmd_instrument('seawinds')
+        ee = pti.get_gcmd_platform('quikscat')
+        provider = metadata['NC_GLOBAL#institution']
+        if provider.lower()=='jpl':
+            provider = 'NASA/JPL/QUIKSCAT'
+        provider = pti.get_gcmd_provider(provider)
+
+        self.dataset.SetMetadataItem('instrument', json.dumps(mm))
+        self.dataset.SetMetadataItem('platform', json.dumps(ee))
+        self.dataset.SetMetadataItem('data_center', json.dumps(provider))
+        self.dataset.SetMetadataItem('entry_title', metadata['NC_GLOBAL#title'])
+        self.dataset.SetMetadataItem('ISO_topic_category',
+                json.dumps(pti.get_iso19115_topic_category('Oceans')))

nansat/mappers/scatterometers.py

@@ -0,0 +1,104 @@
+import json
+import numpy as np
+from datetime import datetime
+import pythesint as pti
+
+import gdal
+
+from nansat.vrt import VRT
+from nansat.geolocation import Geolocation
+from nansat.nsr import NSR
+from nansat.domain import Domain
+from nansat.mappers.mapper_netcdf_cf import Mapper as NetcdfCF
+from nansat.exceptions import WrongMapperError
+
+class Mapper(NetcdfCF):
+    """ Nansat mapper for scatterometers """
+
+    def __init__(self, filename, gdal_dataset, metadata, quartile=0, *args, **kwargs):
+
+        super(Mapper, self).__init__(filename, gdal_dataset, metadata, *args, **kwargs)
+
+        intervals = [0,1,2,3]
+        if not quartile in intervals:
+            raise ValueError('quartile must be one of [0,1,2,3]')
+
+        y_size = self.dataset.RasterYSize/4
+        y_offset = [y_size*qq for qq in intervals][quartile]
+
+        # Crop
+        self.set_offset_size('y', y_offset, y_size)
+
+        # Create band of times
+        # TODO: resolve nansat issue #263 (https://github.com/nansencenter/nansat/issues/263)
+        #import ipdb; ipdb.set_trace()
+        tt = self.times()[int(y_offset) : int(y_offset + y_size)]
+        self.dataset.SetMetadataItem('time_coverage_start', tt[0].astype(datetime).isoformat())
+        self.dataset.SetMetadataItem('time_coverage_end', tt[-1].astype(datetime).isoformat())
+        time_stamps = (tt - tt[0]) / np.timedelta64(1, 's')
+        self.band_vrts['time'] = VRT.from_array(
+                np.tile(time_stamps, (self.dataset.RasterXSize, 1)).transpose()
+            )
+        self.create_band(
+                src = {
+                    'SourceFilename': self.band_vrts['time'].filename,
+                    'SourceBand': 1,
+                },
+                dst = {
+                    'name': 'timestamp',
+                    'time_coverage_start': tt[0].astype(datetime).isoformat(),
+                    'units': 'seconds since time_coverage_start',
+                }
+            )
+
+        # Set projection to wkt
+        #self.dataset.SetProjection(NSR().wkt)
+
+    def set_gcps(self, lon, lat, gdal_dataset):
+        """ Set gcps """
+        self.band_vrts['new_lon_VRT'] = VRT.from_array(lon)
+        self.dataset.SetGCPs(VRT._lonlat2gcps(lon, lat, n_gcps=400), NSR().wkt)
+
+        # Add geolocation from correct longitudes and latitudes
+        self._add_geolocation(
+                Geolocation(self.band_vrts['new_lon_VRT'], self, x_band=1, y_band=self._latitude_band_number(gdal_dataset))
+            )
+
+    def _geoloc_band_number(self, gdal_dataset, sub_filename_index, long_name):
+        """ Return the band number associated to a specific geolocation sub-file and long_name
+        """
+        band_index = [ii for ii, ll in enumerate(self._get_sub_filenames(gdal_dataset)) if
+                ':'+sub_filename_index in ll]
+        if band_index:
+            band_num = band_index[0] + 1
+        else:
+            raise WrongMapperError
+        # Check that it is actually latitudes
+        band = gdal.Open(
+                self._get_sub_filenames(gdal_dataset)[band_index[0]]
+            )
+        if not band.GetRasterBand(1).GetMetadata()['long_name'] == long_name:
+            raise ValueError('Cannot find %s band'%long_name)
+        return band_num
+
+    def _latitude_band_number(self, gdal_dataset):
+        return self._geoloc_band_number(gdal_dataset, 'lat', 'latitude')
+
+    def _longitude_band_number(self, gdal_dataset):
+        return self._geoloc_band_number(gdal_dataset, 'lon', 'longitude')
+
+    @staticmethod
+    def shift_longitudes(lon):
+        """ Apply correction of longitudes (they are defined on 0:360 degrees but also contain
+        egative values)
+
+        TODO: consider making this core to nansat - different ways of defining longitudes (-180:180
+        og 0:360 degrees) often cause problems...
+        """
+        return np.mod(lon+180., 360.) - 180.
+
+    def _create_empty(self, gdal_dataset, gdal_metadata):
+        lat = gdal.Open(
+                self._get_sub_filenames(gdal_dataset)[self._latitude_band_number(gdal_dataset)]
+            )
+        super(NetcdfCF, self).__init__(lat.RasterXSize, lat.RasterYSize, metadata=gdal_metadata)

nansat/nansat.py

@@ -1099,6 +1099,7 @@ def _get_mapper(self, mappername, **kwargs):
         gdal_dataset, metadata = self._get_dataset_metadata()
         tmp_vrt = None
 
+        # TODO: There seems to be code repetition in this if-test - should be avoided...
         if mappername is not '':
             # If a specific mapper is requested, we test only this one.
             # get the module name