diff --git a/nansat/exporter.py b/nansat/exporter.py index 2844e14d..25c50a87 100644 --- a/nansat/exporter.py +++ b/nansat/exporter.py @@ -15,6 +15,7 @@ from __future__ import print_function, absolute_import, division import os +import pytz import tempfile import datetime import warnings @@ -46,7 +47,7 @@ class Exporter(object): '_FillValue', 'type', 'scale', 'offset', 'NETCDF_VARNAME'] def export(self, filename='', bands=None, rm_metadata=None, add_geolocation=True, - driver='netCDF', options='FORMAT=NC4', hardcopy=False): + driver='netCDF', options='FORMAT=NC4', hardcopy=False, add_gcps=True): """Export Nansat object into netCDF or GTiff file Parameters @@ -112,10 +113,11 @@ def export(self, filename='', bands=None, rm_metadata=None, add_geolocation=True if self.filename == filename or hardcopy: export_vrt.hardcopy_bands() - if driver == 'GTiff': - add_gcps = export_vrt.prepare_export_gtiff() - else: - add_gcps = export_vrt.prepare_export_netcdf() + if add_gcps: + if driver == 'GTiff': + add_gcps = export_vrt.prepare_export_gtiff() + else: + add_gcps = export_vrt.prepare_export_netcdf() # Create output file using GDAL dataset = gdal.GetDriverByName(driver).CreateCopy(filename, export_vrt.dataset, @@ -129,7 +131,11 @@ def export(self, filename='', bands=None, rm_metadata=None, add_geolocation=True # Rename variable names to get rid of the band numbers self.rename_variables(filename) # Rename attributes to get rid of "GDAL_" added by gdal - self.correct_attributes(filename, history=self.vrt.dataset.GetMetadata()['history']) + try: + history = self.vrt.dataset.GetMetadata()['history'] + except KeyError: + history = None + self.correct_attributes(filename, history=history) self.logger.debug('Export - OK!') @@ -367,7 +373,7 @@ def _create_dimensions(self, nc_inp, nc_out, time): nc_out.createDimension(dim_name, dim_shapes[dim_name]) # create value for time variable - td = time - datetime.datetime(1900, 1, 1) + td = time - datetime.datetime(1900, 1, 1).replace(tzinfo=pytz.timezone("utc")) days = td.days + (float(td.seconds) / 60.0 / 60.0 / 24.0) # add time dimension nc_out.createDimension('time', 1) diff --git a/nansat/mappers/mapper_meps.py b/nansat/mappers/mapper_meps.py index 841ae178..c756315c 100644 --- a/nansat/mappers/mapper_meps.py +++ b/nansat/mappers/mapper_meps.py @@ -1,4 +1,7 @@ import json +import pytz +import datetime + import pythesint as pti from osgeo import gdal @@ -13,7 +16,8 @@ class Mapper(NetcdfCF, Opendap): - def __init__(self, url, gdal_dataset, gdal_metadata, file_num=0, bands=None, *args, **kwargs): + def __init__(self, url, gdal_dataset, gdal_metadata, file_num=0, bands=None, *args, + **kwargs): if not url.endswith(".nc"): raise WrongMapperError @@ -23,20 +27,23 @@ def __init__(self, url, gdal_dataset, gdal_metadata, file_num=0, bands=None, *ar except OSError: raise WrongMapperError + if "title" not in ds.ncattrs() or "meps" not in ds.getncattr("title").lower(): + raise WrongMapperError + metadata = {} for attr in ds.ncattrs(): - metadata[attr] = ds.getncattr(attr) - - if 'title' not in metadata.keys() or 'meps' not in metadata['title'].lower(): - raise WrongMapperError + content = ds.getncattr(attr) + if isinstance(content, str): + content = content.replace("æ", "ae").replace("ø", "oe").replace("å", "aa") + metadata[attr] = content self.input_filename = url - xsize = ds.dimensions['x'].size - ysize = ds.dimensions['y'].size + xsize = ds.dimensions["x"].size + ysize = ds.dimensions["y"].size # Pick 10 meter height dimension only - height_dim = 'height6' + height_dim = "height6" if height_dim not in ds.dimensions.keys(): raise WrongMapperError if ds.dimensions[height_dim].size != 1: @@ -47,7 +54,7 @@ def __init__(self, url, gdal_dataset, gdal_metadata, file_num=0, bands=None, *ar varnames = [] for var in ds.variables: var_dimensions = ds.variables[var].dimensions - if var_dimensions == ('time', height_dim, 'y', 'x'): + if var_dimensions == ("time", height_dim, "y", "x"): varnames.append(var) # Projection @@ -63,8 +70,8 @@ def __init__(self, url, gdal_dataset, gdal_metadata, file_num=0, bands=None, *ar nsr = NSR(crs.to_proj4()) # Geotransform - xx = ds.variables['x'][0:2] - yy = ds.variables['y'][0:2] + xx = ds.variables["x"][0:2] + yy = ds.variables["y"][0:2] gtrans = xx[0], xx[1]-xx[0], 0, yy[0], 0, yy[1]-yy[0] self._init_from_dataset_params(xsize, ysize, gtrans, nsr.wkt) @@ -79,7 +86,13 @@ def __init__(self, url, gdal_dataset, gdal_metadata, file_num=0, bands=None, *ar self._pop_spatial_dimensions(dimension_names) index = self._get_index_of_dimensions(dimension_names, {}, dim_sizes) fn = self._get_sub_filename(url, band, dim_sizes, index) - meta_dict.append(self.get_band_metadata_dict(fn, ds.variables[band])) + band_metadata = self.get_band_metadata_dict(fn, ds.variables[band]) + # Add time stamp to band metadata + tt = datetime.datetime.fromisoformat(str(self.times()[index["time"]["index"]])) + if tt.tzinfo is None: + tt = pytz.utc.localize(tt) + band_metadata["dst"]["time"] = tt.isoformat() + meta_dict.append(band_metadata) self.create_bands(meta_dict) @@ -89,8 +102,11 @@ def __init__(self, url, gdal_dataset, gdal_metadata, file_num=0, bands=None, *ar # Get dictionary describing the instrument and platform according to # the GCMD keywords - mm = pti.get_gcmd_instrument('computer') - ee = pti.get_gcmd_platform('models') + mm = pti.get_gcmd_instrument("computer") + ee = pti.get_gcmd_platform("models") + + self.dataset.SetMetadataItem("instrument", json.dumps(mm)) + self.dataset.SetMetadataItem("platform", json.dumps(ee)) - self.dataset.SetMetadataItem('instrument', json.dumps(mm)) - self.dataset.SetMetadataItem('platform', json.dumps(ee)) + # Set input filename + self.dataset.SetMetadataItem("nc_file", self.input_filename) diff --git a/nansat/mappers/mapper_meps_ncml.py b/nansat/mappers/mapper_meps_ncml.py index 4fa8602b..c5fe86c7 100644 --- a/nansat/mappers/mapper_meps_ncml.py +++ b/nansat/mappers/mapper_meps_ncml.py @@ -1,17 +1,29 @@ import os +import pytz +import netCDF4 +import datetime + +import numpy as np -from nansat.mappers.mapper_meps import Mapper as NCMapper from nansat.exceptions import WrongMapperError +from nansat.mappers.mapper_meps import Mapper as NCMapper class Mapper(NCMapper): - def __init__(self, ncml_url, gdal_dataset, gdal_metadata, file_num=0, *args, **kwargs): + def __init__(self, ncml_url, gdal_dataset, gdal_metadata, netcdf_dim=None, *args, **kwargs): if not ncml_url.endswith(".ncml"): raise WrongMapperError - url = self._get_odap_url(ncml_url, file_num) + dt = datetime.timedelta(0) + if netcdf_dim is not None and "time" in netcdf_dim.keys(): + ds = netCDF4.Dataset(ncml_url) + time = netcdf_dim["time"].astype(datetime.datetime).replace( + tzinfo=pytz.timezone("utc")) + dt = time - datetime.datetime.fromisoformat(ds.time_coverage_start.replace( + "Z", "+00:00")) + url = self._get_odap_url(ncml_url, np.round(dt.total_seconds()/3600)) super(Mapper, self).__init__(url, gdal_dataset, gdal_metadata, *args, **kwargs) diff --git a/nansat/mappers/mapper_opendap_arome.py b/nansat/mappers/mapper_opendap_arome.py index 2175685f..6f45fdeb 100644 --- a/nansat/mappers/mapper_opendap_arome.py +++ b/nansat/mappers/mapper_opendap_arome.py @@ -1,16 +1,14 @@ -# Name: mapper_arome.py -# Purpose: Nansat mapping for AROME-Arctic and MEPS (MetCoOp Ensemble -# Prediction System) data provided by MET.NO -# Author: Artem Moiseev # Licence: This file is part of NANSAT. You can redistribute it or modify # under the terms of GNU General Public License, v.3 # http://www.gnu.org/licenses/gpl-3.0.html -import json import os +import json +import pytz +import datetime + import numpy as np import pythesint as pti -from datetime import datetime from netCDF4 import Dataset from osgeo import gdal from pyproj import CRS @@ -103,7 +101,13 @@ def __init__(self, filename, gdal_dataset, gdal_metadata, netcdf_dim=None, self._pop_spatial_dimensions(dimension_names) index = self._get_index_of_dimensions(dimension_names, netcdf_dim, dim_sizes) fn = self._get_sub_filename(filename, band, dim_sizes, index) - meta_dict.append(self.get_band_metadata_dict(fn, ds.variables[band])) + band_metadata = self.get_band_metadata_dict(fn, ds.variables[band]) + # Add time stamp to band metadata + tt = datetime.datetime.fromisoformat(str(self.times()[index["time"]["index"]])) + if tt.tzinfo is None: + tt = pytz.utc.localize(tt) + band_metadata["dst"]["time"] = tt.isoformat() + meta_dict.append(band_metadata) self.create_bands(meta_dict) @@ -122,7 +126,7 @@ def get_date(filename): it as a datetime object. """ _, filename = os.path.split(filename) - return datetime.strptime(filename.split('_')[-1], '%Y%m%dT%HZ.nc') + return datetime.datetime.strptime(filename.split('_')[-1], '%Y%m%dT%HZ.nc') def convert_dstime_datetimes(self, ds_time): """Convert time variable to np.datetime64""" diff --git a/nansat/mappers/mapper_opendap_met_nordic.py b/nansat/mappers/mapper_opendap_met_nordic.py new file mode 100644 index 00000000..d17b2b96 --- /dev/null +++ b/nansat/mappers/mapper_opendap_met_nordic.py @@ -0,0 +1,116 @@ +# Licence: This file is part of NANSAT. You can redistribute it or modify +# under the terms of GNU General Public License, v.3 +# http://www.gnu.org/licenses/gpl-3.0.html +import json +import os +import numpy as np +import pythesint as pti + +from datetime import datetime +from netCDF4 import Dataset +from osgeo import gdal +from pyproj import CRS + +from nansat.exceptions import WrongMapperError +from nansat.nsr import NSR +from nansat.vrt import VRT + +from nansat.mappers.mapper_netcdf_cf import Mapper as MapperNCCF +#from nansat.mappers.mapper_arome import Mapper as MapperArome +from nansat.mappers.opendap import Opendap + + +class Mapper(MapperNCCF, Opendap): + + baseURLs = ['https://thredds.met.no/thredds/dodsC/metpparchivev3',] + timeVarName = 'time' + + def __init__(self, filename, gdal_dataset, gdal_metadata, netcdf_dim=None, + ds=None, bands=None, cachedir=None, *args, **kwargs): + + if netcdf_dim is None: + netcdf_dim = {} + + self.input_filename = filename + + if gdal_dataset is None: + gdal_dataset = gdal.Open(filename) + + try: + ds = Dataset(filename) + except OSError: + raise WrongMapperError + + metadata = {} + for attr in ds.ncattrs(): + metadata[attr] = self._fix_encoding(ds.getncattr(attr)) + + if 'title' not in metadata.keys() or 'met nordic' not in metadata['title'].lower(): + raise WrongMapperError + + xsize = ds.dimensions['x'].size + ysize = ds.dimensions['y'].size + + varnames = [] + for var in ds.variables: + var_dimensions = ds.variables[var].dimensions + if var_dimensions == ('time', 'y', 'x'): + varnames.append(var) + + # Projection + try: + grid_mapping = ds.variables[ds.variables[varnames[0]].grid_mapping] + except KeyError: + raise WrongMapperError + + grid_mapping_dict = {} + for index in grid_mapping.ncattrs(): + grid_mapping_dict[index] = grid_mapping.getncattr(index) + crs = CRS.from_cf(grid_mapping_dict) + nsr = NSR(crs.to_proj4()) + + # Geotransform + xx = ds.variables['x'][0:2] + yy = ds.variables['y'][0:2] + gtrans = xx[0], xx[1]-xx[0], 0, yy[0], 0, yy[1]-yy[0] + + self._init_from_dataset_params(xsize, ysize, gtrans, nsr.wkt) + + meta_dict = [] + if bands is None: + bands = varnames + for band in bands: + if band not in ds.variables.keys(): + continue + dimension_names, dim_sizes = self._get_dimension_info(band) + self._pop_spatial_dimensions(dimension_names) + index = self._get_index_of_dimensions(dimension_names, netcdf_dim, dim_sizes) + fn = self._get_sub_filename(filename, band, dim_sizes, index) + meta_dict.append(self.get_band_metadata_dict(fn, ds.variables[band])) + + self.create_bands(meta_dict) + + # Copy metadata + for key in metadata.keys(): + self.dataset.SetMetadataItem(str(key), str(metadata[key])) + + mm = pti.get_gcmd_instrument('Computer') + ee = pti.get_gcmd_platform('ecmwfifs') + self.dataset.SetMetadataItem('instrument', json.dumps(mm)) + self.dataset.SetMetadataItem('platform', json.dumps(ee)) + + @staticmethod + def get_date(filename): + """Extract date and time parameters from filename and return + it as a datetime object. + """ + _, filename = os.path.split(filename) + return datetime.strptime(filename.split('_')[-1], '%Y%m%dT%HZ.nc') + + def convert_dstime_datetimes(self, ds_time): + """Convert time variable to np.datetime64""" + ds_datetimes = np.array( + [(np.datetime64(self.epoch).astype('M8[s]') + + np.timedelta64(int(sec), 's').astype('m8[s]')) for sec in ds_time]).astype('M8[s]') + return ds_datetimes + diff --git a/nansat/mappers/opendap.py b/nansat/mappers/opendap.py index 83bcb812..c0d604ec 100644 --- a/nansat/mappers/opendap.py +++ b/nansat/mappers/opendap.py @@ -109,9 +109,7 @@ def get_dataset_time(self): if os.path.exists(cachefile): ds_time = np.load(cachefile)[self.timeVarName] else: - warnings.warn('Time consuming loading time from OpenDAP...') ds_time = self.ds.variables[self.timeVarName][:] - warnings.warn('Loading time - OK!') if os.path.exists(cachefile): np.savez(cachefile, **{self.timeVarName: ds_time}) diff --git a/nansat/mappers/sentinel1.py b/nansat/mappers/sentinel1.py index 4c14732e..57e7faba 100644 --- a/nansat/mappers/sentinel1.py +++ b/nansat/mappers/sentinel1.py @@ -1,4 +1,5 @@ import json +import logging import numpy as np from dateutil.parser import parse try: @@ -107,13 +108,65 @@ def get_gcps(self, flip_gcp_line=False): return gcps - def add_incidence_angle_band(self): + def get_gcp_shape(self): + """Get GCP shape from GCP pixel and line data. Returns the + GCP y and x dimensions. + """ # Get GCP variables pixel = self.ds['GCP_pixel_'+self.ds.polarisation[:2]][:].data line = self.ds['GCP_line_'+self.ds.polarisation[:2]][:].data + gcp_y = np.unique(line[:].data).shape[0] + gcp_x = np.unique(pixel[:].data).shape[0] + """ Uniqueness is not guaranteed at the correct grid - assert + that the gcp_index dimension divided by gcp_x and gcp_y + dimensions results in an integer. + """ + def test1(gcp_y, gcp_x): + """ Check if one of them is correct, then modify the + other. If that does not work, use test2 function. + """ + if pixel.size % gcp_x != 0: + if pixel.size % gcp_y == 0: + gcp_x = pixel.size/gcp_y + if pixel.size % gcp_y != 0: + if pixel.size % gcp_x == 0: + gcp_y = pixel.size/gcp_x + + if gcp_y*gcp_x != pixel.size: + gcp_x = test2(gcp_x) + gcp_y = test2(gcp_y) + + return gcp_y, gcp_x + + def test2(gcp_dim): + """Test modulo from -/+gcp_*/4 until remainder=0 + """ + if pixel.size % gcp_dim != 0: + for i in range(np.round(gcp_dim/4).astype("int")): + test_dim = gcp_dim - i + if pixel.size % test_dim == 0: + gcp_dim = test_dim + break + test_dim = gcp_dim + i + if pixel.size % test_dim == 0: + gcp_dim = test_dim + break + return gcp_dim + + gcp_y, gcp_x = test1(gcp_y, gcp_x) + + logging.debug("GCPY size: %d" % gcp_y) + logging.debug("GCPX size: %d" % gcp_x) + + if gcp_y*gcp_x != pixel.size: + raise ValueError("GCP dimension mismatch") + + return int(gcp_y), int(gcp_x) + + def add_incidence_angle_band(self): + gcp_y, gcp_x = self.get_gcp_shape() inci = self.ds['GCP_incidenceAngle_'+self.ds.polarisation[:2]][:].data - inci = inci.reshape(np.unique(line[:].data).shape[0], - np.unique(pixel[:].data).shape[0]) + inci = inci.reshape(gcp_y, gcp_x) # Add incidence angle band inciVRT = VRT.from_array(inci) @@ -128,14 +181,11 @@ def add_incidence_angle_band(self): def get_full_size_GCPs(self): # Get GCP variables - pixel = self.ds['GCP_pixel_' + self.ds.polarisation[:2]][:].data - line = self.ds['GCP_line_' + self.ds.polarisation[:2]][:].data + gcp_y, gcp_x = self.get_gcp_shape() lon = self.ds['GCP_longitude_' + self.ds.polarisation[:2]][:].data + lon = lon.reshape(gcp_y, gcp_x) lat = self.ds['GCP_latitude_' + self.ds.polarisation[:2]][:].data - lon = lon.reshape(np.unique(line[:].data).shape[0], - np.unique(pixel[:].data).shape[0]) - lat = lat.reshape(np.unique(line[:].data).shape[0], - np.unique(pixel[:].data).shape[0]) + lat = lat.reshape(gcp_y, gcp_x) return lon, lat def add_look_direction_band(self):