From 56996cc888ef301d73e418b0099caab1f5fc2d84 Mon Sep 17 00:00:00 2001 From: "Morten W. Hansen" Date: Fri, 23 Oct 2020 16:09:13 +0200 Subject: [PATCH 01/57] test that history metadata is added correctly --- nansat/tests/test_exporter.py | 9 +++++++++ 1 file changed, 9 insertions(+) diff --git a/nansat/tests/test_exporter.py b/nansat/tests/test_exporter.py index 40067a4f2..b1437436d 100644 --- a/nansat/tests/test_exporter.py +++ b/nansat/tests/test_exporter.py @@ -67,6 +67,15 @@ def test_special_characters_in_exported_metadata(self): '"'})) self.assertIsInstance(dd, dict) + def test_history_metadata(self): + orig = Nansat(self.test_file_gcps, mapper=self.default_mapper) + hh = datetime.datetime.now().replace(tzinfo=datetime.timezone.utc).isoformat() + orig.vrt.dataset.SetMetadataItem('history', hh) + orig.export(self.tmp_filename) + copy = Nansat(self.tmp_filename, mapper=self.default_mapper) + history = copy.get_metadata('history') + self.assertEqual(history, hh) + def test_time_coverage_metadata_of_exported_equals_original(self): orig = Nansat(self.test_file_gcps, mapper=self.default_mapper) orig.set_metadata('time_coverage_start', '2010-01-02T08:49:02.347809') From 7e5a1f6d1709a672190968f5c1540268b065b917 Mon Sep 17 00:00:00 2001 From: "Morten W. Hansen" Date: Fri, 11 Jun 2021 15:11:34 +0200 Subject: [PATCH 02/57] updated baseURLs to s1 thredds --- nansat/mappers/mapper_opendap_sentinel1.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/nansat/mappers/mapper_opendap_sentinel1.py b/nansat/mappers/mapper_opendap_sentinel1.py index a38a4df5d..9a2a68402 100644 --- a/nansat/mappers/mapper_opendap_sentinel1.py +++ b/nansat/mappers/mapper_opendap_sentinel1.py @@ -32,6 +32,8 @@ class Mapper(Opendap, Sentinel1): baseURLs = [ + 'https://nbstds.met.no/thredds/dodsC/NBS/S1A', + 'https://nbstds.met.no/thredds/dodsC/NBS/S1B', 'http://nbstds.met.no/thredds/dodsC/NBS/S1A', 'http://nbstds.met.no/thredds/dodsC/NBS/S1B', ] From fceb4b5c30d3097bef664afc86b80eb0f9baeca3 Mon Sep 17 00:00:00 2001 From: "Morten W. Hansen" Date: Wed, 26 Jan 2022 15:08:34 +0000 Subject: [PATCH 03/57] removed print line which hampers readability - the line was printed every time a nansat object was opened --- nansat/mappers/mapper_globcolour_l3m.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/nansat/mappers/mapper_globcolour_l3m.py b/nansat/mappers/mapper_globcolour_l3m.py index 207565795..826e43f5d 100644 --- a/nansat/mappers/mapper_globcolour_l3m.py +++ b/nansat/mappers/mapper_globcolour_l3m.py @@ -25,7 +25,7 @@ def __init__(self, filename, gdalDataset, gdalMetadata, **kwargs): ''' GLOBCOLOR L3M VRT ''' try: - print("=>%s<=" % gdalMetadata['NC_GLOBAL#title']) + x = "=>%s<=" % gdalMetadata['NC_GLOBAL#title'] except (TypeError, KeyError): raise WrongMapperError From 6d715b9905d8888f1a19c70ab0f8e10ddc6a6aa9 Mon Sep 17 00:00:00 2001 From: "Morten W. Hansen" Date: Wed, 26 Jan 2022 15:09:08 +0000 Subject: [PATCH 04/57] added exception handling of value error --- nansat/mappers/mapper_netcdf_cf.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/nansat/mappers/mapper_netcdf_cf.py b/nansat/mappers/mapper_netcdf_cf.py index acaee716e..ba87e99fd 100644 --- a/nansat/mappers/mapper_netcdf_cf.py +++ b/nansat/mappers/mapper_netcdf_cf.py @@ -317,7 +317,7 @@ def _band_dict(self, subfilename, band_num, subds, band=None, band_metadata=None try: band_metadata['time_iso_8601'] = self._time_count_to_np_datetime64( band_metadata[timecountname]) - except KeyError as e: + except (ValueError, KeyError) as e: # No timing information available for this band - it is # probably a constant, such as land area fraction or similar. # Then we don't need time for this band... From 1a53f4bf6e3dbe10b507b02d233d2d284eca80e2 Mon Sep 17 00:00:00 2001 From: "Morten W. Hansen" Date: Wed, 26 Jan 2022 17:28:37 +0000 Subject: [PATCH 05/57] use != instead of is not. Mappername default changed to None would be better. --- nansat/nansat.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/nansat/nansat.py b/nansat/nansat.py index 2c5c5306c..29e1e354a 100644 --- a/nansat/nansat.py +++ b/nansat/nansat.py @@ -1122,7 +1122,7 @@ def _get_mapper(self, mappername, **kwargs): tmp_vrt = None # TODO: There seems to be code repetition in this if-test - should be avoided... - if mappername is not '': + if mappername != '': # If a specific mapper is requested, we test only this one. # get the module name mappername = 'mapper_' + mappername.replace('mapper_', '').replace('.py', '').lower() From d90e2e3a157bdb77bed077da7bcfe521079294dd Mon Sep 17 00:00:00 2001 From: "Morten W. Hansen" Date: Thu, 27 Jan 2022 08:45:55 +0000 Subject: [PATCH 06/57] This is not an error. Might be a warning but it is still annoying, so I've changed it to a debug message. --- nansat/nansat.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/nansat/nansat.py b/nansat/nansat.py index 29e1e354a..47a35e47f 100644 --- a/nansat/nansat.py +++ b/nansat/nansat.py @@ -1067,7 +1067,7 @@ def _get_dataset_metadata(self): try: gdal_dataset = gdal.Open(self.filename) except RuntimeError: - self.logger.error('GDAL could not open %s, trying to read with Nansat mappers...' + self.logger.debug('GDAL could not open %s, trying to read with Nansat mappers...' % self.filename) if gdal_dataset is not None: # get metadata from the GDAL dataset From 21e4334125438d2e524d6dd18796f2bf6831662b Mon Sep 17 00:00:00 2001 From: "Morten W. Hansen" Date: Fri, 23 Dec 2022 11:04:59 +0100 Subject: [PATCH 07/57] #525: added function to export with xarray, some cleaning, and tests --- nansat/exporter.py | 74 +++++++++++++++++++++++++++++- nansat/mappers/mapper_netcdf_cf.py | 1 + nansat/nansat.py | 7 --- nansat/tests/test_exporter.py | 55 +++++++++++++++++++++- 4 files changed, 128 insertions(+), 9 deletions(-) diff --git a/nansat/exporter.py b/nansat/exporter.py index a0ac0fde0..4446e6127 100644 --- a/nansat/exporter.py +++ b/nansat/exporter.py @@ -18,6 +18,7 @@ import tempfile import datetime import warnings +import importlib from nansat.utils import gdal import numpy as np @@ -30,6 +31,11 @@ from nansat.exceptions import NansatGDALError +try: + import xarray as xr +except: + warnings.warn("'xarray' needs to be installed for Exporter.xr_export to work.") + class Exporter(object): """Abstract class for export functions """ @@ -37,7 +43,73 @@ class Exporter(object): DEFAULT_SOURCE = 'satellite remote sensing' UNWANTED_METADATA = ['dataType', 'SourceFilename', 'SourceBand', '_Unsigned', 'FillValue', - 'time', '_FillValue', 'type', 'scale', 'offset'] + '_FillValue', 'type', 'scale', 'offset', 'NETCDF_VARNAME'] + + def xr_export(self, filename, bands=None, encoding=None): + """Export Nansat object into netCDF using xarray. Longitude + and latitude arrays will be used as coordinates (of the 2D + dataset). Note that ACDD and CF metadata should be added to + the Nansat object before exporting. This method only removes + the UNWANTED_METADATA and creates the NetCDF file with the + metadata provided in the Nansat object. + + Parameters + ---------- + filename : str + output file name + bands: list (default=None) + Specify band numbers to export. + If None, all bands are exported. + encoding: dict + Nested dictionary with variable names as keys and + dictionaries of variable specific encodings as values, + e.g., {"my_variable": { + "dtype": "int16", + "scale_factor": 0.1, + "zlib": True, + "_FillValue": 9999, + }, ...} + The h5netcdf engine supports both the NetCDF4-style + compression encoding parameters + {"zlib": True, "complevel": 9} and the h5py ones + {"compression": "gzip", "compression_opts": 9}. This + allows using any compression plugin installed in the HDF5 + library, e.g. LZF. + + """ + xr_installed = importlib.util.find_spec("xarray") + if xr_installed is None: + raise ModuleNotFoundError("Please install 'xarray'") + datavars = {} + lon, lat = self.get_geolocation_grids() + for band in self.bands().keys(): + band_metadata = self.get_metadata(band_id=band) + # Clean up the use metadata + pop_keys = [] + for key in band_metadata.keys(): + if key in self.UNWANTED_METADATA: + pop_keys.append(key) + for key in pop_keys: + band_metadata.pop(key) + if bands is not None: + if (not band in bands) and (not band_metadata["name"] in bands): + continue + datavars[band_metadata["name"]] = xr.DataArray( + self[band], + dims = ["x", "y"], + attrs = band_metadata) + ds = xr.Dataset( + datavars, + coords = {"longitude": (["x", "y"], lon), "latitude": (["x", "y"], lat)}, + attrs = self.get_metadata()) + if encoding is None: + # Having a coordinate variable with a fill value applied + # violates the cf standard + encoding = { + "longitude": {"_FillValue": None}, + "latitude": {"_FillValue": None} + } + ds.to_netcdf(filename, encoding=encoding) def export(self, filename='', bands=None, rm_metadata=None, add_geolocation=True, driver='netCDF', options=None, hardcopy=False): diff --git a/nansat/mappers/mapper_netcdf_cf.py b/nansat/mappers/mapper_netcdf_cf.py index 68990d397..5ed0f43c8 100644 --- a/nansat/mappers/mapper_netcdf_cf.py +++ b/nansat/mappers/mapper_netcdf_cf.py @@ -277,6 +277,7 @@ def get_band_number(): subds = gdal.Open(fn) band = subds.GetRasterBand(Context.band_number) band_metadata = self._clean_band_metadata(band) + band_metadata['_FillValue'] = band.GetNoDataValue() return self._band_dict(fn, Context.band_number, subds, band=band, band_metadata=band_metadata) diff --git a/nansat/nansat.py b/nansat/nansat.py index 47a35e47f..7fabc6dd3 100644 --- a/nansat/nansat.py +++ b/nansat/nansat.py @@ -101,9 +101,6 @@ class Nansat(Domain, Exporter): """ - FILL_VALUE = 9.96921e+36 - ALT_FILL_VALUE = -10000. - # instance attributes logger = None filename = None @@ -265,10 +262,6 @@ def _fill_with_nan(self, band, band_data): """Fill input array with fill value taen from input band metadata""" fill_value = float(band.GetMetadata()['_FillValue']) band_data[band_data == fill_value] = np.nan - # quick hack to avoid problem with wrong _FillValue - see issue - # #123 - if fill_value == self.FILL_VALUE: - band_data[band_data == self.ALT_FILL_VALUE] = np.nan return band_data diff --git a/nansat/tests/test_exporter.py b/nansat/tests/test_exporter.py index a01222182..02ec036dc 100644 --- a/nansat/tests/test_exporter.py +++ b/nansat/tests/test_exporter.py @@ -37,7 +37,7 @@ from netCDF4 import Dataset -from nansat import Nansat, Domain, NSR +from nansat import Nansat, Domain, NSR, exporter from nansat.utils import gdal from nansat.tests.nansat_test_base import NansatTestBase @@ -46,6 +46,59 @@ class ExporterTest(NansatTestBase): + @patch('nansat.exporter.importlib.util.find_spec') + def test_xr_export__raises_module_not_found_error(self, mock_find_spec): + """Test that an error is raised if xarray is not installed.""" + mock_find_spec.return_value = None + n = Nansat(self.test_file_arctic) + with self.assertRaises(ModuleNotFoundError) as e: + n.xr_export('test.nc') + self.assertEqual(str(e.exception), "Please install 'xarray'") + + @patch.object(exporter.xr.Dataset, 'to_netcdf') + def test_xr_export__default(self, mock_to_netcdf): + """Test that the to_netcdf function is called for default + usage of xr_export.""" + n = Nansat(self.test_file_arctic) + n.xr_export('test.nc') + mock_to_netcdf.assert_called_once_with('test.nc', + encoding={"longitude": {"_FillValue": None}, "latitude": {"_FillValue": None}}) + + def test_xr_export__one_band(self): + """Test that only a given band is exported.""" + n = Nansat(self.test_file_arctic) + fd, tmp_ncfile = tempfile.mkstemp(suffix='.nc') + n.xr_export(tmp_ncfile, bands=['Bootstrap']) + ds = Dataset(tmp_ncfile) + self.assertEqual(list(ds.variables.keys()), ['Bootstrap', 'longitude', 'latitude']) + os.close(fd) + os.unlink(tmp_ncfile) + + def test_xr_export__with_specific_encoding_and_nan_values(self): + """Test that a band with nan-values is masked as expected, + and with the FillValue specified by the user.""" + n = Nansat(self.test_file_arctic) + xx = n['Bootstrap'].astype(float) + xx = np.ma.masked_where( + xx == float(n.get_metadata(band_id='Bootstrap', key='_FillValue')), xx) + xx.data[xx.mask] = np.nan + yy = xx.data + yy[2,:] = np.nan + n.add_band(yy, parameters={'name': 'test_band_with_nans'}) + encoding = { + "longitude": {"_FillValue": None}, + "latitude": {"_FillValue": None}, + "test_band_with_nans": {"_FillValue": 999}, + } + + fd, tmp_ncfile = tempfile.mkstemp(suffix='.nc') + n.xr_export(tmp_ncfile, encoding=encoding) + ds = Dataset(tmp_ncfile) + # Check that original invalid/masked data equals 999 + self.assertEqual(ds.variables['test_band_with_nans'][:].data[xx.mask][0], 999.) + # Check that data set to np.nan equals 999 + self.assertEqual(ds.variables['test_band_with_nans'][:].data[2,0], 999.) + def test_geolocation_of_exportedNC_vs_original(self): """ Lon/lat in original and exported file should coincide """ orig = Nansat(self.test_file_gcps, mapper=self.default_mapper) From 0ca6093e862c7f1d8cfe6c6d479b136b7dcd57f0 Mon Sep 17 00:00:00 2001 From: "Morten W. Hansen" Date: Fri, 23 Dec 2022 11:09:17 +0100 Subject: [PATCH 08/57] #525: ipdb lines were not meant to be committed.. --- nansat/tests/test_exporter.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/nansat/tests/test_exporter.py b/nansat/tests/test_exporter.py index 02ec036dc..ca6335027 100644 --- a/nansat/tests/test_exporter.py +++ b/nansat/tests/test_exporter.py @@ -78,6 +78,8 @@ def test_xr_export__with_specific_encoding_and_nan_values(self): """Test that a band with nan-values is masked as expected, and with the FillValue specified by the user.""" n = Nansat(self.test_file_arctic) + import ipdb + ipdb.set_trace() xx = n['Bootstrap'].astype(float) xx = np.ma.masked_where( xx == float(n.get_metadata(band_id='Bootstrap', key='_FillValue')), xx) From d1be8e9ac575f55d59fc5403af3028d1e5945b5c Mon Sep 17 00:00:00 2001 From: "Morten W. Hansen" Date: Fri, 23 Dec 2022 11:11:43 +0100 Subject: [PATCH 09/57] #525: removed ipdb lines --- nansat/tests/test_exporter.py | 2 -- 1 file changed, 2 deletions(-) diff --git a/nansat/tests/test_exporter.py b/nansat/tests/test_exporter.py index ca6335027..02ec036dc 100644 --- a/nansat/tests/test_exporter.py +++ b/nansat/tests/test_exporter.py @@ -78,8 +78,6 @@ def test_xr_export__with_specific_encoding_and_nan_values(self): """Test that a band with nan-values is masked as expected, and with the FillValue specified by the user.""" n = Nansat(self.test_file_arctic) - import ipdb - ipdb.set_trace() xx = n['Bootstrap'].astype(float) xx = np.ma.masked_where( xx == float(n.get_metadata(band_id='Bootstrap', key='_FillValue')), xx) From 110789cfe38a77fb909bff26902cff818116f830 Mon Sep 17 00:00:00 2001 From: "Morten W. Hansen" Date: Sun, 25 Dec 2022 13:45:12 +0100 Subject: [PATCH 10/57] #525: remover xarray based export function and modified the export function instead --- nansat/exporter.py | 136 ++++++++++++++++------------------ nansat/tests/test_exporter.py | 77 +++++++++---------- 2 files changed, 98 insertions(+), 115 deletions(-) diff --git a/nansat/exporter.py b/nansat/exporter.py index 4446e6127..e7ffe317b 100644 --- a/nansat/exporter.py +++ b/nansat/exporter.py @@ -45,74 +45,8 @@ class Exporter(object): UNWANTED_METADATA = ['dataType', 'SourceFilename', 'SourceBand', '_Unsigned', 'FillValue', '_FillValue', 'type', 'scale', 'offset', 'NETCDF_VARNAME'] - def xr_export(self, filename, bands=None, encoding=None): - """Export Nansat object into netCDF using xarray. Longitude - and latitude arrays will be used as coordinates (of the 2D - dataset). Note that ACDD and CF metadata should be added to - the Nansat object before exporting. This method only removes - the UNWANTED_METADATA and creates the NetCDF file with the - metadata provided in the Nansat object. - - Parameters - ---------- - filename : str - output file name - bands: list (default=None) - Specify band numbers to export. - If None, all bands are exported. - encoding: dict - Nested dictionary with variable names as keys and - dictionaries of variable specific encodings as values, - e.g., {"my_variable": { - "dtype": "int16", - "scale_factor": 0.1, - "zlib": True, - "_FillValue": 9999, - }, ...} - The h5netcdf engine supports both the NetCDF4-style - compression encoding parameters - {"zlib": True, "complevel": 9} and the h5py ones - {"compression": "gzip", "compression_opts": 9}. This - allows using any compression plugin installed in the HDF5 - library, e.g. LZF. - - """ - xr_installed = importlib.util.find_spec("xarray") - if xr_installed is None: - raise ModuleNotFoundError("Please install 'xarray'") - datavars = {} - lon, lat = self.get_geolocation_grids() - for band in self.bands().keys(): - band_metadata = self.get_metadata(band_id=band) - # Clean up the use metadata - pop_keys = [] - for key in band_metadata.keys(): - if key in self.UNWANTED_METADATA: - pop_keys.append(key) - for key in pop_keys: - band_metadata.pop(key) - if bands is not None: - if (not band in bands) and (not band_metadata["name"] in bands): - continue - datavars[band_metadata["name"]] = xr.DataArray( - self[band], - dims = ["x", "y"], - attrs = band_metadata) - ds = xr.Dataset( - datavars, - coords = {"longitude": (["x", "y"], lon), "latitude": (["x", "y"], lat)}, - attrs = self.get_metadata()) - if encoding is None: - # Having a coordinate variable with a fill value applied - # violates the cf standard - encoding = { - "longitude": {"_FillValue": None}, - "latitude": {"_FillValue": None} - } - ds.to_netcdf(filename, encoding=encoding) - def export(self, filename='', bands=None, rm_metadata=None, add_geolocation=True, - driver='netCDF', options=None, hardcopy=False): + driver='netCDF', options='FORMAT=NC4', hardcopy=False): """Export Nansat object into netCDF or GTiff file Parameters @@ -129,7 +63,8 @@ def export(self, filename='', bands=None, rm_metadata=None, add_geolocation=True add geolocation array datasets to exported file? driver : str Name of GDAL driver (format) - options : str or list + options : str or list (default: 'FORMAT=NC4' for NetCDF 4 + file format) GDAL export options in format of: 'OPT=VAL', or ['OPT1=VAL1', 'OP2='VAL2'] See also http://www.gdal.org/frmt_netcdf.html @@ -142,11 +77,6 @@ def export(self, filename='', bands=None, rm_metadata=None, add_geolocation=True Notes ------ - If number of bands is more than one, serial numbers are added at the end of each band name. - It is possible to fix it by changing line.4605 in GDAL/frmts/netcdf/netcdfdataset.cpp : - 'if( nBands > 1 ) sprintf(szBandName,"%s%d",tmpMetadata,iBand);' - --> 'if( nBands > 1 ) sprintf(szBandName,"%s",tmpMetadata);' - CreateCopy fails in case the band name has special characters, e.g. the slash in 'HH/VV'. @@ -188,14 +118,72 @@ def export(self, filename='', bands=None, rm_metadata=None, add_geolocation=True add_gcps = export_vrt.prepare_export_netcdf() # Create output file using GDAL - dataset = gdal.GetDriverByName(driver).CreateCopy(filename, export_vrt.dataset, options=options) + dataset = gdal.GetDriverByName(driver).CreateCopy(filename, export_vrt.dataset, + options=options) del dataset # add GCPs into netCDF file as separate float variables if add_gcps: Exporter._add_gcps(filename, export_vrt.dataset.GetGCPs()) + if driver=='netCDF': + # 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.rename_attributes(filename) + self.logger.debug('Export - OK!') + def rename_attributes(self, filename): + """ Rename global attributes to get rid of the "GDAL_"-string + added by gdal. + """ + GDAL = "GDAL_" + del_attrs = [] + rename_attrs = [] + # Open new file to edit attribute names + ds = Dataset(filename, 'r+') + for attr in ds.ncattrs(): + if GDAL in attr: + if attr.replace(GDAL, "") in ds.ncattrs(): + # Mark for deletion + del_attrs.append(attr) + else: + # Mark for renaming + rename_attrs.append(attr) + + # Delete repeated attributes.. + for attr in del_attrs: + ds.delncattr(attr) + # Rename attributes: + for attr in rename_attrs: + ds.renameAttribute(attr, attr.replace(GDAL, "")) + ds.close() + + def rename_variables(self, filename): + """ Rename variable names to reflect the name attribute of + the variable's metadata. + + Parameters + ---------- + filename : str + NetCDF file name + """ + # Open new file to edit variable names + ds = Dataset(filename, 'r+') + + # Decide which variables to rename + rename_vars = [] + for var in ds.variables.keys(): + if ('name' in ds.variables[var].ncattrs()) and ( + var != ds.variables[var].getncattr('name')): + rename_vars.append(var) + + # Rename selected variables + for var in rename_vars: + ds.renameVariable(var, ds.variables[var].getncattr('name')) + + ds.close() + def export2thredds(self, filename, bands=None, diff --git a/nansat/tests/test_exporter.py b/nansat/tests/test_exporter.py index 02ec036dc..d715acc11 100644 --- a/nansat/tests/test_exporter.py +++ b/nansat/tests/test_exporter.py @@ -46,37 +46,11 @@ class ExporterTest(NansatTestBase): - @patch('nansat.exporter.importlib.util.find_spec') - def test_xr_export__raises_module_not_found_error(self, mock_find_spec): - """Test that an error is raised if xarray is not installed.""" - mock_find_spec.return_value = None - n = Nansat(self.test_file_arctic) - with self.assertRaises(ModuleNotFoundError) as e: - n.xr_export('test.nc') - self.assertEqual(str(e.exception), "Please install 'xarray'") - - @patch.object(exporter.xr.Dataset, 'to_netcdf') - def test_xr_export__default(self, mock_to_netcdf): - """Test that the to_netcdf function is called for default - usage of xr_export.""" - n = Nansat(self.test_file_arctic) - n.xr_export('test.nc') - mock_to_netcdf.assert_called_once_with('test.nc', - encoding={"longitude": {"_FillValue": None}, "latitude": {"_FillValue": None}}) - - def test_xr_export__one_band(self): - """Test that only a given band is exported.""" - n = Nansat(self.test_file_arctic) - fd, tmp_ncfile = tempfile.mkstemp(suffix='.nc') - n.xr_export(tmp_ncfile, bands=['Bootstrap']) - ds = Dataset(tmp_ncfile) - self.assertEqual(list(ds.variables.keys()), ['Bootstrap', 'longitude', 'latitude']) - os.close(fd) - os.unlink(tmp_ncfile) - - def test_xr_export__with_specific_encoding_and_nan_values(self): + def test_export__with_nan_values(self): """Test that a band with nan-values is masked as expected, - and with the FillValue specified by the user.""" + that global attribute names are the same as the Nansat + metadata, and that variable names are the same as the Nansat + band names without appended numbers.""" n = Nansat(self.test_file_arctic) xx = n['Bootstrap'].astype(float) xx = np.ma.masked_where( @@ -85,19 +59,36 @@ def test_xr_export__with_specific_encoding_and_nan_values(self): yy = xx.data yy[2,:] = np.nan n.add_band(yy, parameters={'name': 'test_band_with_nans'}) - encoding = { - "longitude": {"_FillValue": None}, - "latitude": {"_FillValue": None}, - "test_band_with_nans": {"_FillValue": 999}, - } + + zz = np.zeros(yy.shape) + zz[np.isnan(yy)] = 1 + sumnan = zz.sum() + # temp file for exported netcdf fd, tmp_ncfile = tempfile.mkstemp(suffix='.nc') - n.xr_export(tmp_ncfile, encoding=encoding) + + # export with nansat + n.export(tmp_ncfile) + ds = Dataset(tmp_ncfile) - # Check that original invalid/masked data equals 999 - self.assertEqual(ds.variables['test_band_with_nans'][:].data[xx.mask][0], 999.) - # Check that data set to np.nan equals 999 - self.assertEqual(ds.variables['test_band_with_nans'][:].data[2,0], 999.) + # Check that the number of elements equal to np.nan is the same + xx = ds.variables['test_band_with_nans'][:] + yy = np.zeros(xx.shape) + yy[np.isnan(xx)] = 1 + self.assertEqual(yy.sum(), sumnan) + + # Check that the global metadata attribute names are the same + orig_metadata = list(n.get_metadata().keys()) + ncattrs = ds.ncattrs() + for attr in orig_metadata: + self.assertIn(attr, ncattrs) + + # Check that variable names don't contain band numbers + self.assertEqual(list(ds.variables.keys()), ['polar_stereographic', 'x', 'y', + 'UMass_AES', 'Bootstrap', 'Bristol', 'test_band_with_nans']) + + os.close(fd) + os.unlink(tmp_ncfile) def test_geolocation_of_exportedNC_vs_original(self): """ Lon/lat in original and exported file should coincide """ @@ -120,6 +111,8 @@ def test_special_characters_in_exported_metadata(self): '"'})) self.assertIsInstance(dd, dict) + # Date format in the history could be version dependent + # - probably fine to skip this test.. def test_history_metadata(self): orig = Nansat(self.test_file_gcps, mapper=self.default_mapper) hh = datetime.datetime.now().replace(tzinfo=datetime.timezone.utc).isoformat() @@ -471,7 +464,9 @@ def test_example1(self): def test_example2(self): n = Nansat(self.tmp_ncfile) - res = n.export2thredds(self.filename_exported, {'x_wind_10m': {'description': 'example'}}) + res = n.export2thredds( + self.filename_exported, + {'x_wind_10m': {'description': 'example'}}) self.assertEqual(res, None) def test_example3(self): From 2c453635f4e71f7abb08bac7fd080949421fdd6d Mon Sep 17 00:00:00 2001 From: "Morten W. Hansen" Date: Sun, 25 Dec 2022 13:49:04 +0100 Subject: [PATCH 11/57] #525: removed unnecessary test --- nansat/tests/test_exporter.py | 11 ----------- 1 file changed, 11 deletions(-) diff --git a/nansat/tests/test_exporter.py b/nansat/tests/test_exporter.py index d715acc11..78f42893c 100644 --- a/nansat/tests/test_exporter.py +++ b/nansat/tests/test_exporter.py @@ -111,17 +111,6 @@ def test_special_characters_in_exported_metadata(self): '"'})) self.assertIsInstance(dd, dict) - # Date format in the history could be version dependent - # - probably fine to skip this test.. - def test_history_metadata(self): - orig = Nansat(self.test_file_gcps, mapper=self.default_mapper) - hh = datetime.datetime.now().replace(tzinfo=datetime.timezone.utc).isoformat() - orig.vrt.dataset.SetMetadataItem('history', hh) - orig.export(self.tmp_filename) - copy = Nansat(self.tmp_filename, mapper=self.default_mapper) - history = copy.get_metadata('history') - self.assertEqual(history, hh) - def test_time_coverage_metadata_of_exported_equals_original(self): orig = Nansat(self.test_file_gcps, mapper=self.default_mapper) orig.set_metadata('time_coverage_start', '2010-01-02T08:49:02.347809') From d4d64951b59abc81a1989c901e03b2447389e515 Mon Sep 17 00:00:00 2001 From: "Morten W. Hansen" Date: Sun, 25 Dec 2022 15:48:12 +0100 Subject: [PATCH 12/57] #525: cleaned test code and added one test to demonstrate issue when fill values are not correctly picked up --- nansat/tests/test_exporter.py | 41 ++++++++++++++++++++++------------- 1 file changed, 26 insertions(+), 15 deletions(-) diff --git a/nansat/tests/test_exporter.py b/nansat/tests/test_exporter.py index 78f42893c..1ec085b88 100644 --- a/nansat/tests/test_exporter.py +++ b/nansat/tests/test_exporter.py @@ -52,17 +52,17 @@ def test_export__with_nan_values(self): metadata, and that variable names are the same as the Nansat band names without appended numbers.""" n = Nansat(self.test_file_arctic) - xx = n['Bootstrap'].astype(float) - xx = np.ma.masked_where( - xx == float(n.get_metadata(band_id='Bootstrap', key='_FillValue')), xx) - xx.data[xx.mask] = np.nan - yy = xx.data - yy[2,:] = np.nan - n.add_band(yy, parameters={'name': 'test_band_with_nans'}) + aa = n['Bootstrap'].astype(float) + aa = np.ma.masked_where( + aa == float(n.get_metadata(band_id='Bootstrap', key='_FillValue')), aa) + aa.data[aa.mask] = np.nan + bb = aa.data.copy() + bb[2,:] = np.nan + n.add_band(bb, parameters={'name': 'test_band_with_nans'}) - zz = np.zeros(yy.shape) - zz[np.isnan(yy)] = 1 - sumnan = zz.sum() + cc = np.zeros(bb.shape) + cc[np.isnan(bb)] = 1 + sumnan = cc.sum() # temp file for exported netcdf fd, tmp_ncfile = tempfile.mkstemp(suffix='.nc') @@ -71,11 +71,12 @@ def test_export__with_nan_values(self): n.export(tmp_ncfile) ds = Dataset(tmp_ncfile) - # Check that the number of elements equal to np.nan is the same - xx = ds.variables['test_band_with_nans'][:] - yy = np.zeros(xx.shape) - yy[np.isnan(xx)] = 1 - self.assertEqual(yy.sum(), sumnan) + # Check that the number of elements in bb and the new file + # equal to np.nan is the same when opened with netCDF4.Dataset: + dd = ds.variables['test_band_with_nans'][:] + ee = np.zeros(dd.shape) + ee[np.isnan(dd)] = 1 + self.assertEqual(ee.sum(), sumnan) # Check that the global metadata attribute names are the same orig_metadata = list(n.get_metadata().keys()) @@ -87,6 +88,16 @@ def test_export__with_nan_values(self): self.assertEqual(list(ds.variables.keys()), ['polar_stereographic', 'x', 'y', 'UMass_AES', 'Bootstrap', 'Bristol', 'test_band_with_nans']) + # Open the tmp file using nansat, and check that the Bootstrap + # band from the tmp file is the same as the one in the original + n = Nansat(tmp_ncfile) + ff = n['Bootstrap'] + ff = np.ma.masked_where( + ff == float(n.get_metadata(band_id='Bootstrap', key='_FillValue')), ff) + # This only works when the fill value (-10000) in arctic.nc is + # treated correctly (with python3.9 and GDAL3.4.1) + self.assertEqual(ff.mask[0,239], aa.mask[0,239]) + os.close(fd) os.unlink(tmp_ncfile) From f276d0ee05ce0126c9b681a1861cf52c1117f6e6 Mon Sep 17 00:00:00 2001 From: "Morten W. Hansen" Date: Mon, 26 Dec 2022 15:07:50 +0100 Subject: [PATCH 13/57] #525: static methods --- nansat/exporter.py | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/nansat/exporter.py b/nansat/exporter.py index e7ffe317b..810349ae2 100644 --- a/nansat/exporter.py +++ b/nansat/exporter.py @@ -133,7 +133,8 @@ def export(self, filename='', bands=None, rm_metadata=None, add_geolocation=True self.logger.debug('Export - OK!') - def rename_attributes(self, filename): + @staticmethod + def rename_attributes(filename): """ Rename global attributes to get rid of the "GDAL_"-string added by gdal. """ @@ -159,7 +160,8 @@ def rename_attributes(self, filename): ds.renameAttribute(attr, attr.replace(GDAL, "")) ds.close() - def rename_variables(self, filename): + @staticmethod + def rename_variables(filename): """ Rename variable names to reflect the name attribute of the variable's metadata. From a671dc8f6a75c98a08cf88db3a8eb600dee553a9 Mon Sep 17 00:00:00 2001 From: "Morten W. Hansen" Date: Tue, 27 Dec 2022 10:04:22 +0100 Subject: [PATCH 14/57] #525: adjusted netcdf-cf mapper to account for new attribute names without 'GDAL_' --- nansat/mappers/mapper_netcdf_cf.py | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/nansat/mappers/mapper_netcdf_cf.py b/nansat/mappers/mapper_netcdf_cf.py index 5ed0f43c8..d899e1264 100644 --- a/nansat/mappers/mapper_netcdf_cf.py +++ b/nansat/mappers/mapper_netcdf_cf.py @@ -41,8 +41,10 @@ def __init__(self, filename, gdal_dataset, gdal_metadata, *args, **kwargs): raise WrongMapperError if 'NC_GLOBAL#GDAL_NANSAT_GCPY_000' in list(gdal_metadata.keys()) or \ - 'NC_GLOBAL#GDAL_NANSAT_GCPProjection' in list(gdal_metadata.keys()): - # Probably Nansat generated netcdf of swath data - see issue #192 + 'NC_GLOBAL#GDAL_NANSAT_GCPProjection' in list(gdal_metadata.keys()) or \ + 'NC_GLOBAL#NANSAT_GCPY_000' in list(gdal_metadata.keys()) or \ + 'NC_GLOBAL#NANSAT_GCPProjection' in list(gdal_metadata.keys()): + # Nansat generated netcdf of swath data is not standard CF raise WrongMapperError metadata = VRT._remove_strings_in_metadata_keys(gdal_metadata, From 54f17ab0009964d41fc9ad37c37d5c93610b82fa Mon Sep 17 00:00:00 2001 From: "Morten W. Hansen" Date: Wed, 4 Jan 2023 22:16:12 +0000 Subject: [PATCH 15/57] #525: gdal adds its own Conventions attribute with value 'CF-1.5' which is obviously wrong. This is now deleted. Same goes for the history. That should also be handled in the processing software rather then by the gdal netcdf driver. --- nansat/exporter.py | 15 ++++++++++----- 1 file changed, 10 insertions(+), 5 deletions(-) diff --git a/nansat/exporter.py b/nansat/exporter.py index 810349ae2..47a87d667 100644 --- a/nansat/exporter.py +++ b/nansat/exporter.py @@ -143,14 +143,19 @@ def rename_attributes(filename): rename_attrs = [] # Open new file to edit attribute names ds = Dataset(filename, 'r+') + """ The netcdf driver adds the Conventions attribute with + value CF-1.5. This may be wrong, so it is better to use the + Conventions metadata from the Nansat object. Other attributes + added by gdal that are already present in Nansat, should also + be deleted.""" for attr in ds.ncattrs(): if GDAL in attr: if attr.replace(GDAL, "") in ds.ncattrs(): - # Mark for deletion - del_attrs.append(attr) - else: - # Mark for renaming - rename_attrs.append(attr) + # Mark the attribute created by the netcdf driver + # for deletion - ref above comment + del_attrs.append(attr.replace(GDAL, "")) + # Mark for renaming + rename_attrs.append(attr) # Delete repeated attributes.. for attr in del_attrs: From cfc5a285c79e6bac42cf6f6acafb4c3efb569fc5 Mon Sep 17 00:00:00 2001 From: "Morten W. Hansen" Date: Thu, 5 Jan 2023 13:32:34 +0000 Subject: [PATCH 16/57] #525: history is cut by gdal CreateCopy. This needs to be overridden. --- nansat/exporter.py | 21 ++++++++++++++------- 1 file changed, 14 insertions(+), 7 deletions(-) diff --git a/nansat/exporter.py b/nansat/exporter.py index 47a87d667..28b126f2e 100644 --- a/nansat/exporter.py +++ b/nansat/exporter.py @@ -129,14 +129,18 @@ 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.rename_attributes(filename) + self.correct_attributes(filename, history=self.vrt.dataset.GetMetadata()['history']) self.logger.debug('Export - OK!') @staticmethod - def rename_attributes(filename): + def correct_attributes(filename, history=None): """ Rename global attributes to get rid of the "GDAL_"-string - added by gdal. + added by gdal, remove attributes added by gdal that are + already present in the Nansat object, and correct the history + attribute (the latter may be reduced in length because + gdal.GetDriverByName(driver).CreateCopy limits the string + lenght to 161 characters). """ GDAL = "GDAL_" del_attrs = [] @@ -144,10 +148,10 @@ def rename_attributes(filename): # Open new file to edit attribute names ds = Dataset(filename, 'r+') """ The netcdf driver adds the Conventions attribute with - value CF-1.5. This may be wrong, so it is better to use the - Conventions metadata from the Nansat object. Other attributes - added by gdal that are already present in Nansat, should also - be deleted.""" + value CF-1.5. This in most cases wrong, so it is better to use + the Conventions metadata from the Nansat object. Other + attributes added by gdal that are already present in Nansat, + should also be deleted.""" for attr in ds.ncattrs(): if GDAL in attr: if attr.replace(GDAL, "") in ds.ncattrs(): @@ -163,6 +167,9 @@ def rename_attributes(filename): # Rename attributes: for attr in rename_attrs: ds.renameAttribute(attr, attr.replace(GDAL, "")) + # Correct the history + if history is not None: + ds.history = history ds.close() @staticmethod From 5562633c29f2f3edd423c668e60835f4556e55cd Mon Sep 17 00:00:00 2001 From: "Morten W. Hansen" Date: Thu, 1 Jun 2023 13:23:02 +0000 Subject: [PATCH 17/57] get proj4 string from grid mapping variable --- nansat/mappers/mapper_opendap_arome.py | 30 ++++++++++++++++++++------ 1 file changed, 23 insertions(+), 7 deletions(-) diff --git a/nansat/mappers/mapper_opendap_arome.py b/nansat/mappers/mapper_opendap_arome.py index 0d86fb13a..800270ee3 100644 --- a/nansat/mappers/mapper_opendap_arome.py +++ b/nansat/mappers/mapper_opendap_arome.py @@ -5,17 +5,19 @@ # 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 pyproj import CRS from nansat.mappers.mapper_arome import Mapper as MapperArome from nansat.mappers.opendap import Opendap from nansat.exceptions import WrongMapperError from nansat.nsr import NSR -import pythesint as pti -import os -from datetime import datetime -from netCDF4 import Dataset -import numpy as np -import json class Mapper(Opendap, MapperArome): @@ -37,10 +39,24 @@ def __init__(self, filename, gdal_dataset, gdal_metadata, date=None, timestamp = date if date else self.get_date(filename) ds = Dataset(filename) + shape=[] + for key in ds.dimensions.keys(): + shape.append(ds.dimensions[key].size) + + for var in ds.variables.keys(): + if ds[var].shape==tuple(shape): + break + try: - self.srcDSProjection = NSR(ds.variables['projection_lambert'].proj4) + grid_mapping = ds.variables[ds.variables[var].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) + + self.srcDSProjection = NSR(crs.to_proj4()) self.create_vrt(filename, gdal_dataset, gdal_metadata, timestamp, ds, bands, cachedir) From f797b687442b912850bffc5e66244e662f3962e8 Mon Sep 17 00:00:00 2001 From: "Morten W. Hansen" Date: Tue, 6 Jun 2023 11:57:36 +0000 Subject: [PATCH 18/57] specify bands in opendap arome mapper since files are too big --- nansat/mappers/mapper_opendap_arome.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/nansat/mappers/mapper_opendap_arome.py b/nansat/mappers/mapper_opendap_arome.py index 800270ee3..86e6e5f85 100644 --- a/nansat/mappers/mapper_opendap_arome.py +++ b/nansat/mappers/mapper_opendap_arome.py @@ -35,6 +35,9 @@ class Mapper(Opendap, MapperArome): def __init__(self, filename, gdal_dataset, gdal_metadata, date=None, ds=None, bands=None, cachedir=None, *args, **kwargs): + if bands is None: + bands = ['x_wind_10m', 'y_wind_10m'] + self.test_mapper(filename) timestamp = date if date else self.get_date(filename) ds = Dataset(filename) From eb5e6fbf8cccfb4c8bf1b2fb521044ef16d4b676 Mon Sep 17 00:00:00 2001 From: "Morten W. Hansen" Date: Tue, 6 Jun 2023 14:09:33 +0000 Subject: [PATCH 19/57] update to new netcdf attributes --- nansat/mappers/sentinel1.py | 12 ++++++++---- 1 file changed, 8 insertions(+), 4 deletions(-) diff --git a/nansat/mappers/sentinel1.py b/nansat/mappers/sentinel1.py index f64b355de..4c14732ee 100644 --- a/nansat/mappers/sentinel1.py +++ b/nansat/mappers/sentinel1.py @@ -69,10 +69,14 @@ def set_gcmd_dif_keywords(self): self.dataset.SetMetadataItem('instrument', json.dumps(mm)) self.dataset.SetMetadataItem('platform', json.dumps(ee)) - self.dataset.SetMetadataItem('time_coverage_start', - self.dataset.GetMetadataItem('ACQUISITION_START_TIME')) - self.dataset.SetMetadataItem('time_coverage_end', - self.dataset.GetMetadataItem('ACQUISITION_STOP_TIME')) + tstart = self.dataset.GetMetadataItem('time_coverage_start') if \ + self.dataset.GetMetadataItem('ACQUISITION_START_TIME') is None else \ + self.dataset.GetMetadataItem('ACQUISITION_START_TIME') + tend = self.dataset.GetMetadataItem('time_coverage_end') if \ + self.dataset.GetMetadataItem('ACQUISITION_STOP_TIME') is None else \ + self.dataset.GetMetadataItem('ACQUISITION_STOP_TIME') + self.dataset.SetMetadataItem('time_coverage_start', tstart) + self.dataset.SetMetadataItem('time_coverage_end', tend) def get_gcps(self, flip_gcp_line=False): """ Get Ground Control Points for the dataset. From aebc93c5198951a5e542af1d7349d30fc6bd8816 Mon Sep 17 00:00:00 2001 From: "Morten W. Hansen" Date: Tue, 13 Jun 2023 10:11:15 +0000 Subject: [PATCH 20/57] working mapper for opendap arome arctic --- nansat/mappers/mapper_arome.py | 21 ++- nansat/mappers/mapper_netcdf_cf.py | 194 ++++++++++++++++++------- nansat/mappers/mapper_opendap_arome.py | 177 +++++++++++++++------- nansat/mappers/opendap.py | 5 +- 4 files changed, 272 insertions(+), 125 deletions(-) diff --git a/nansat/mappers/mapper_arome.py b/nansat/mappers/mapper_arome.py index e3a4569f0..c7e4af60e 100644 --- a/nansat/mappers/mapper_arome.py +++ b/nansat/mappers/mapper_arome.py @@ -8,23 +8,22 @@ class Mapper(NetcdfCF): - def __init__(self, *args, **kwargs): + def __init__(self, filename, gdal_dataset, gdal_metadata, *args, **kwargs): - mm = args[2] # metadata - if not mm: + if not gdal_metadata: raise WrongMapperError - if 'NC_GLOBAL#source' not in list(mm.keys()): + if 'source' not in list(gdal_metadata.keys()): raise WrongMapperError - if not 'arome' in mm['NC_GLOBAL#source'].lower() and \ - not 'meps' in mm['NC_GLOBAL#source'].lower(): + if not 'arome' in gdal_metadata['title'].lower() and \ + not 'meps' in gdal_metadata['title'].lower(): raise WrongMapperError - super(Mapper, self).__init__(*args, **kwargs) + super(Mapper, self).__init__(filename, gdal_dataset, gdal_metadata, *args, **kwargs) - self.dataset.SetMetadataItem('time_coverage_start', - (parse(mm['NC_GLOBAL#min_time'], ignoretz=True, fuzzy=True).isoformat())) - self.dataset.SetMetadataItem('time_coverage_end', - (parse(mm['NC_GLOBAL#max_time'], ignoretz=True, fuzzy=True).isoformat())) + #self.dataset.SetMetadataItem('time_coverage_start', parse( + # gdal_metadata['NC_GLOBAL#min_time'], ignoretz=True, fuzzy=True).isoformat()) + #self.dataset.SetMetadataItem('time_coverage_end', parse( + # gdal_metadata['NC_GLOBAL#max_time'], ignoretz=True, fuzzy=True).isoformat())) # Get dictionary describing the instrument and platform according to # the GCMD keywords diff --git a/nansat/mappers/mapper_netcdf_cf.py b/nansat/mappers/mapper_netcdf_cf.py index d899e1264..8128507ee 100644 --- a/nansat/mappers/mapper_netcdf_cf.py +++ b/nansat/mappers/mapper_netcdf_cf.py @@ -4,13 +4,17 @@ http://cfconventions.org/compliance-checker.html ''' -import warnings, os, datetime -import numpy as np import collections -from nansat.utils import gdal +import datetime +import os +import warnings + +import numpy as np from dateutil.parser import parse from netCDF4 import Dataset +from osgeo import gdal +from pyproj import CRS from nansat.vrt import VRT from nansat.nsr import NSR @@ -28,7 +32,6 @@ class ContinueI(Exception): class Mapper(VRT): """ """ - input_filename = '' def __init__(self, filename, gdal_dataset, gdal_metadata, *args, **kwargs): @@ -124,9 +127,10 @@ def _time_reference(self, ds=None): rt = parse(times.units, fuzzy=True) # This sets timezone to local # Remove timezone information from epoch, which defaults to # utc (otherwise the timezone should be given in the dataset) + units = times.units epoch = datetime.datetime(rt.year, rt.month, rt.day, rt.hour, rt.minute, rt.second) - return epoch, times.units + return epoch, units def _timevarname(self, ds=None): if not ds: @@ -171,9 +175,10 @@ def _time_count_to_np_datetime64(self, time_count, time_reference=None): raise Exception('Check time units..') return tt - 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. + def _band_list(self, gdal_dataset, gdal_metadata, netcdf_dim=None, bands=None, *args, + **kwargs): + ''' Create list of dictionaries mapping source and destination + metadata of bands that should be added to the Nansat object. Parameters ---------- @@ -182,50 +187,48 @@ def _band_list(self, gdal_dataset, gdal_metadata, netcdf_dim={}, bands=[], *args gdal_metadata : dict Dictionary of global metadata netcdf_dim : dict - Dictionary of desired slice of a multi-dimensional array. Since - gdal only returns 2D bands, a multi-dimensional array (x,y,z) is - split into z bands accompanied with metadata information about the - position of the slice along the z-axis. The (key, value) pairs represent dimension name and - the desired value in that dimension (not the index), respectively. E.g., for a height - dimension of [10, 20, 30] m, netcdf_dim = {'height': 20} if you want to extract the - data from 20 m height. + Dictionary of desired slice of a multi-dimensional array. + Since gdal only returns 2D bands, a multi-dimensional + array (x,y,z) is split into z bands accompanied with + metadata information about the position of the slice along + the z-axis. The (key, value) pairs represent dimension + name and the desired value in that dimension (not the + index), respectively. E.g., for a height dimension of + [10, 20, 30] m, netcdf_dim = {'height': 20} if you want to + extract the data from 20 m height. bands : list - List of desired bands following NetCDF-CF standard names. NOTE: - some datasets have other bands as well, i.e., of data not yet - implemented in CF. We may at some point generalize this to provide - a dict with key name and value, where the key is, e.g., - "standard_name" or "metno_name", etc. + List of desired bands following NetCDF-CF standard names. + NOTE: some datasets have other bands as well, i.e., of + data not yet implemented in CF. We may at some point + generalize this to provide a dict with key name and value, + where the key is, e.g., "standard_name" or "metno_name", + etc. ''' + if netcdf_dim is None: + netcdf_dim = {} + if bands is None: + ds = Dataset(self.input_filename) + bands = ds.variables.keys() metadictlist = [] for fn in self._get_sub_filenames(gdal_dataset): + band = fn.split(':')[-1] + if band not in bands: + continue + # Don't process geolocation variables/bands if ('GEOLOCATION_X_DATASET' in fn or 'longitude' in fn or 'GEOLOCATION_Y_DATASET' in fn or 'latitude' in fn): continue try: - metadictlist.append(self._get_band_from_subfile(fn, netcdf_dim=netcdf_dim, bands=bands)) + metadictlist.append( + self._get_band_metadata(fn, netcdf_dim=netcdf_dim)) except ContinueI: continue - return metadictlist - - def _get_band_from_subfile(self, fn, netcdf_dim={}, bands=[]): - nc_ds = Dataset(self.input_filename) - band_name = fn.split(':')[-1] - if bands: - variable = nc_ds.variables[band_name] - if 'standard_name' not in variable.ncattrs() or not variable.standard_name in bands: - raise ContinueI - # TODO: consider to allow band name in addition or instead of standard_name in the band - # list kwarg - sub_band = nc_ds.variables[band_name] - dimension_names = [b.name for b in sub_band.get_dims()] - dimension_names.reverse() - dim_sizes = {} - for dim in sub_band.get_dims(): - dim_sizes[dim.name] = dim.size - - # Pop spatial dimensions (longitude and latitude, or x and y) + @staticmethod + def _pop_spatial_dimensions(dimension_names): + """ Pop spatial dimensions (longitude and latitude, or x and y) + """ for allowed in ALLOWED_SPATIAL_DIMENSIONS_X: try: ind_dim_x = [i for i, s in enumerate(dimension_names) if allowed in s.lower()][0] @@ -238,15 +241,29 @@ def _get_band_from_subfile(self, fn, netcdf_dim={}, bands=[]): dimension_names.pop(ind_dim_y) except IndexError: continue + + def _get_index_of_dimensions(self, dimension_names, netcdf_dim, dim_sizes): + """ Nansat only works with 2D data. For for datasets with + higher dimensions, we need to pick the correct slice. The + function returns a dictionary with the index of the wanted + data slice, as specified in the netcdf_dim dictionary. + """ + ds = Dataset(self.input_filename) index4key = collections.OrderedDict() for key in dimension_names: if key in netcdf_dim.keys(): val = netcdf_dim[key] - if key == 'time' and type(val) == np.datetime64: + if key == 'time': + if type(val) != datetime.datetime and type(val) != np.datetime64: + raise ValueError + if type(val) == datetime.datetime: + val = np.datetime64(val).astype('M8[s]') # Get band number from given timestamp index = int(np.argmin(np.abs(self.times() - val))) else: - index = int(np.argmin(np.abs(nc_ds.variables[key][:] - val))) + import pdb + pdb.set_trace() + index = int(np.argmin(np.abs(ds.variables[key][:] - val))) index4key[key] = { 'index': index, 'size': dim_sizes[key], @@ -256,32 +273,75 @@ def _get_band_from_subfile(self, fn, netcdf_dim={}, bands=[]): 'index': 0, 'size': dim_sizes[key], } - - # Works in Python 2 and 3 + return index4key + + @staticmethod + def _calculate_band_number(index4dimension, dimension_names): + """ Uses index and size values in provided dict to calculate + gdal band number. The index4dimension dict contains the index + in the dimensions given by dimension_names, and the size of + the dimensions. + """ class Context: band_number = 1 multiplier = 1 - #band_number = 1 # Only works in python 3 - #multiplier = 1 # Only works in Python 3 + def get_band_number(): - #nonlocal band_number # Only works in Python 3 - #nonlocal multiplier # Only works in Python 3 try: name_dim0 = dimension_names.pop(0) except: return - Context.band_number += index4key[name_dim0]['index']*Context.multiplier - Context.multiplier *= index4key[name_dim0]['size'] + Context.band_number += index4dimension[name_dim0]['index']*Context.multiplier + Context.multiplier *= index4dimension[name_dim0]['size'] get_band_number() get_band_number() + return Context.band_number + + def _get_dimension_info(self, band): + """ Get band list from a netCDF4.Dataset. See docs of + _band_list. + """ + ds = Dataset(self.input_filename) + var = ds.variables[band] + dimension_names = [b.name for b in var.get_dims()] + dimension_names.reverse() + dim_sizes = collections.OrderedDict() + for dim in var.get_dims(): + dim_sizes[dim.name] = dim.size + + return dimension_names, dim_sizes + + def _get_band_number(self, band, netcdf_dim={}): + """ Get gdal band number from multidimensional dataset. + """ + if netcdf_dim is None: + netcdf_dim = {} + + 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) + band_number = self._calculate_band_number(index, dimension_names) + + return band_number + + def _get_band_metadata(self, fn, netcdf_dim=None): + """ Gets band metadata using gdal + """ + if netcdf_dim is None: + netcdf_dim = {} + import pdb + pdb.set_trace() #check fn from file and from url + band = fn.split(':')[-1] + band_number = self._get_band_number(band, netcdf_dim) + subds = gdal.Open(fn) - band = subds.GetRasterBand(Context.band_number) + band = subds.GetRasterBand(band_number) band_metadata = self._clean_band_metadata(band) band_metadata['_FillValue'] = band.GetNoDataValue() - return self._band_dict(fn, Context.band_number, subds, band=band, + return self._band_dict(fn, band_number, subds, band=band, band_metadata=band_metadata) def _clean_band_metadata(self, band, remove = ['_Unsigned', 'ScaleRatio', @@ -394,15 +454,37 @@ def _create_empty(self, gdal_dataset, gdal_metadata): 'spatial reference WKT, assuming a regular longitude/latitude grid') self._create_empty_with_nansat_spatial_reference_WKT(gdal_dataset, gdal_metadata) - def _create_empty_from_projection_variable(self, gdal_dataset, gdal_metadata, - projection_variable='projection_lambert'): + def _create_empty_from_projection_variable(self, gdal_dataset, gdal_metadata): + + # TODO: only open dataset once or close them.. ds = Dataset(self.input_filename) - subdataset = gdal.Open(self._get_sub_filenames(gdal_dataset)[0]) + + shape=[] + for key in ds.dimensions.keys(): + shape.append(ds.dimensions[key].size) + + for var in ds.variables.keys(): + if ds[var].shape==tuple(shape): + break + try: + grid_mapping = ds.variables[ds.variables[var].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) + + try: + subdataset = gdal.Open(self._get_sub_filenames(gdal_dataset)[0]) + except IndexError: + subdataset = gdal_dataset self._init_from_dataset_params( x_size = subdataset.RasterXSize, y_size = subdataset.RasterYSize, geo_transform = subdataset.GetGeoTransform(), - projection = NSR(ds.variables[projection_variable].proj4).wkt, + projection = NSR(crs.to_proj4()).wkt, metadata = gdal_metadata) def _create_empty_from_subdatasets(self, gdal_dataset, metadata): diff --git a/nansat/mappers/mapper_opendap_arome.py b/nansat/mappers/mapper_opendap_arome.py index 86e6e5f85..01b4acb9c 100644 --- a/nansat/mappers/mapper_opendap_arome.py +++ b/nansat/mappers/mapper_opendap_arome.py @@ -12,15 +12,19 @@ from datetime import datetime from netCDF4 import Dataset +from osgeo import gdal from pyproj import CRS -from nansat.mappers.mapper_arome import Mapper as MapperArome -from nansat.mappers.opendap import Opendap 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(Opendap, MapperArome): +class Mapper(MapperNCCF): baseURLs = ['http://thredds.met.no/thredds/catalog/arome25/catalog.html', 'https://thredds.met.no/thredds/dodsC/aromearcticarchive', @@ -28,90 +32,151 @@ class Mapper(Opendap, MapperArome): 'https://thredds.met.no/thredds/dodsC/meps25epsarchive', 'http://thredds.met.no/thredds/dodsC/meps25epsarchive'] timeVarName = 'time' - xName = 'x' - yName = 'y' - timeCalendarStart = '1970-01-01' - def __init__(self, filename, gdal_dataset, gdal_metadata, date=None, + def __init__(self, filename, gdal_dataset, gdal_metadata, netcdf_dim=None, ds=None, bands=None, cachedir=None, *args, **kwargs): - if bands is None: - bands = ['x_wind_10m', 'y_wind_10m'] + 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 - self.test_mapper(filename) - timestamp = date if date else self.get_date(filename) - ds = Dataset(filename) + metadata = {} + for attr in ds.ncattrs(): + metadata[attr] = ds.getncattr(attr) - shape=[] - for key in ds.dimensions.keys(): - shape.append(ds.dimensions[key].size) + if not 'arome' in metadata['title'].lower() and \ + not 'meps' in metadata['title'].lower(): + raise WrongMapperError - for var in ds.variables.keys(): - if ds[var].shape==tuple(shape): - break + xsize = ds.dimensions['x'].size + ysize = ds.dimensions['y'].size + # Pick 10 meter height dimension only + height_dim = 'height7' + if height_dim not in ds.dimensions.keys(): + raise WrongMapperError + if ds.dimensions[height_dim].size != 1: + raise WrongMapperError + if ds.variables['height7'][0].data != 10: + raise WrongMapperError + varnames = [] + for var in ds.variables: + var_dimensions = ds.variables[var].dimensions + if var_dimensions == ('time', 'height7', 'y', 'x'): + varnames.append(var) + + # Projection try: - grid_mapping = ds.variables[ds.variables[var].grid_mapping] + 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()) - self.srcDSProjection = 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.create_vrt(filename, gdal_dataset, gdal_metadata, timestamp, ds, bands, cachedir) + self._init_from_dataset_params(xsize, ysize, gtrans, nsr.wkt) + + meta_dict = [] + if bands is None: + bands = varnames + for band in bands: + 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(attr), 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)) - md_item = 'Data Center' - if not self.dataset.GetMetadataItem(md_item): - self.dataset.SetMetadataItem(md_item, 'NO/MET') - md_item = 'Entry Title' - if not self.dataset.GetMetadataItem(md_item): - self.dataset.SetMetadataItem(md_item, str(ds.getncattr('title'))) - md_item = 'summary' - if not self.dataset.GetMetadataItem(md_item): - summary = """ - AROME_Arctic is a convection-permitting atmosphere model covering parts of the Barents - Sea and the Nordic Arctic. It has horizontal resolution of 2.5 km and 65 vertical - levels. AROME_Arctic runs for 66 hours four times a day (00,06,12,18) with three-hourly - cycling for data assimilation. Boundary data is from ECMWF. Model code based on HARMONIE - cy40h1.1 - """ - self.dataset.SetMetadataItem(md_item, str(summary)) + #md_item = 'Data Center' + #if not self.dataset.GetMetadataItem(md_item): + # self.dataset.SetMetadataItem(md_item, 'NO/MET') + #md_item = 'Entry Title' + #if not self.dataset.GetMetadataItem(md_item): + # self.dataset.SetMetadataItem(md_item, str(ds.getncattr('title'))) + #md_item = 'summary' + #if not self.dataset.GetMetadataItem(md_item): + # summary = """ + # AROME_Arctic is a convection-permitting atmosphere model covering parts of the Barents + # Sea and the Nordic Arctic. It has horizontal resolution of 2.5 km and 65 vertical + # levels. AROME_Arctic runs for 66 hours four times a day (00,06,12,18) with three-hourly + # cycling for data assimilation. Boundary data is from ECMWF. Model code based on HARMONIE + # cy40h1.1 + # """ + # self.dataset.SetMetadataItem(md_item, str(summary)) + + def get_band_metadata_dict(self, fn, ncvar): + gds = gdal.Open(fn) + meta_item = { + 'src': {'SourceFilename': fn, 'SourceBand': 1}, + 'dst': {'name': ncvar.name, 'dataType': 6} + } + + for attr_key in ncvar.ncattrs(): + attr_val = ncvar.getncattr(attr_key) + if attr_key in ['scale', 'scale_factor']: + meta_item['src']['ScaleRatio'] = attr_val + elif attr_key in ['offset', 'add_offset']: + meta_item['src']['ScaleOffset'] = attr_val + else: + meta_item['dst'][attr_key] = attr_val + + return meta_item - @staticmethod - def get_date(filename): - """Extract date and time parameters from filename and return - it as a formatted string - Parameters - ---------- - filename: str - nn - - Returns - ------- - str, YYYY-mm-ddThh:MMZ + @staticmethod + def _get_sub_filename(url, var, dim_sizes, index): + """ Opendap driver refers to subdatasets differently than the + standard way in vrt.py + """ + shape = [] + for item in dim_sizes.items(): + if item[0] in index.keys(): + shape.append(index[item[0]]['index']) + else: + shape.append(item[0]) + # assemble dimensions string + gd_shape = ''.join(['[%s]' % dimsize for dimsize in shape]) + return '{url}?{var}.{var}{shape}'.format(url=url, var=var, shape=gd_shape) - Examples - -------- - >>> Mapper.get_date('/path/to/arome_arctic_full_2_5km_20171030T21Z.nc') - '2017-10-30T21:00Z' + @staticmethod + def get_date(filename): + """Extract date and time parameters from filename and return + it as a datetime object. """ _, filename = os.path.split(filename) - t = datetime.strptime(filename.split('_')[-1], '%Y%m%dT%HZ.nc') - return datetime.strftime(t, '%Y-%m-%dT%H:%MZ') + 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.timeCalendarStart).astype('M8[s]') + [(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 88814e5af..c7707bb32 100644 --- a/nansat/mappers/opendap.py +++ b/nansat/mappers/opendap.py @@ -157,8 +157,9 @@ def get_metaitem(self, url, var_name, var_dimensions): # assemble dimensions string dims = ''.join(['[%s]' % dim for dim in var_dimensions]) sfname = '{url}?{var}.{var}{shape}'.format(url=url, var=var_name, shape=dims) - # For Sentinel-1, the source filename is not at the same format. Simple solution is to check - # if this is correct witha try-except but that may be too time consuming. Expecting + # For Sentinel-1, the source filename is not at the same + # format. Simple solution is to check if this is correct with + # a try-except but that may be too time consuming. Expecting # discussion... try: ds = gdal.Open(sfname) From 7876ece2ba34af3222652e2925787f08798afe14 Mon Sep 17 00:00:00 2001 From: "Morten W. Hansen" Date: Wed, 14 Jun 2023 13:31:58 +0000 Subject: [PATCH 21/57] opendap-arome now works --- nansat/domain.py | 4 ++-- nansat/mappers/mapper_netcdf_cf.py | 2 -- nansat/mappers/mapper_opendap_arome.py | 10 ++++++---- nansat/mappers/opendap.py | 6 +++++- 4 files changed, 13 insertions(+), 9 deletions(-) diff --git a/nansat/domain.py b/nansat/domain.py index b38ff0dd6..5ac51b9c7 100644 --- a/nansat/domain.py +++ b/nansat/domain.py @@ -526,9 +526,9 @@ def _create_extent_dict(extent_str): raise ValueError(' must contains exactly 2 parameters ' '("-te" or "-lle") and ("-ts" or "-tr")') key, extent_dict = Domain._add_to_dict(extent_dict, option) - if key is 'te' or key is 'lle': + if key == 'te' or key == 'lle': Domain._validate_te_lle(extent_dict[key]) - elif key is 'ts' or key is 'tr': + elif key == 'ts' or key == 'tr': Domain._validate_ts_tr(extent_dict[key]) return extent_dict diff --git a/nansat/mappers/mapper_netcdf_cf.py b/nansat/mappers/mapper_netcdf_cf.py index 8128507ee..4227bfb98 100644 --- a/nansat/mappers/mapper_netcdf_cf.py +++ b/nansat/mappers/mapper_netcdf_cf.py @@ -331,8 +331,6 @@ def _get_band_metadata(self, fn, netcdf_dim=None): """ if netcdf_dim is None: netcdf_dim = {} - import pdb - pdb.set_trace() #check fn from file and from url band = fn.split(':')[-1] band_number = self._get_band_number(band, netcdf_dim) diff --git a/nansat/mappers/mapper_opendap_arome.py b/nansat/mappers/mapper_opendap_arome.py index 01b4acb9c..8575fb63e 100644 --- a/nansat/mappers/mapper_opendap_arome.py +++ b/nansat/mappers/mapper_opendap_arome.py @@ -24,7 +24,7 @@ from nansat.mappers.opendap import Opendap -class Mapper(MapperNCCF): +class Mapper(MapperNCCF, Opendap): baseURLs = ['http://thredds.met.no/thredds/catalog/arome25/catalog.html', 'https://thredds.met.no/thredds/dodsC/aromearcticarchive', @@ -51,7 +51,7 @@ def __init__(self, filename, gdal_dataset, gdal_metadata, netcdf_dim=None, metadata = {} for attr in ds.ncattrs(): - metadata[attr] = ds.getncattr(attr) + metadata[attr] = self._fix_encoding(ds.getncattr(attr)) if not 'arome' in metadata['title'].lower() and \ not 'meps' in metadata['title'].lower(): @@ -97,6 +97,8 @@ def __init__(self, filename, gdal_dataset, gdal_metadata, netcdf_dim=None, 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) @@ -107,7 +109,7 @@ def __init__(self, filename, gdal_dataset, gdal_metadata, netcdf_dim=None, # Copy metadata for key in metadata.keys(): - self.dataset.SetMetadataItem(str(attr), str(metadata[key])) + self.dataset.SetMetadataItem(str(key), str(metadata[key])) mm = pti.get_gcmd_instrument('Computer') ee = pti.get_gcmd_platform('ecmwfifs') @@ -139,7 +141,7 @@ def get_band_metadata_dict(self, fn, ncvar): } for attr_key in ncvar.ncattrs(): - attr_val = ncvar.getncattr(attr_key) + attr_val = self._fix_encoding(ncvar.getncattr(attr_key)) if attr_key in ['scale', 'scale_factor']: meta_item['src']['ScaleRatio'] = attr_val elif attr_key in ['offset', 'add_offset']: diff --git a/nansat/mappers/opendap.py b/nansat/mappers/opendap.py index c7707bb32..418c21868 100644 --- a/nansat/mappers/opendap.py +++ b/nansat/mappers/opendap.py @@ -211,7 +211,11 @@ def _fix_encoding(var): >>> Opendap._fix_encoding('asnes') 'asnes' """ - return str(var.encode('ascii', 'ignore').decode()) + if isinstance(var, str): + retval = str(var.encode('ascii', 'ignore').decode()) + else: + retval = var + return retval def create_vrt(self, filename, gdalDataset, gdalMetadata, date, ds, bands, cachedir): """ Create VRT From d1bffc64b3f1af5dfb66423c4bbcd036c47806c7 Mon Sep 17 00:00:00 2001 From: "Morten W. Hansen" Date: Wed, 21 Jun 2023 11:29:48 +0200 Subject: [PATCH 22/57] update box --- Vagrantfile | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Vagrantfile b/Vagrantfile index aacc03103..33cfd808c 100644 --- a/Vagrantfile +++ b/Vagrantfile @@ -6,8 +6,8 @@ VAGRANTFILE_API_VERSION = "2" Vagrant.configure(VAGRANTFILE_API_VERSION) do |config| - config.vm.box = "ubuntu/trusty64" - config.vm.box_url = "https://atlas.hashicorp.com/ubuntu/trusty64" + config.vm.box = "ubuntu/bionic64" + #config.vm.box_url = "https://atlas.hashicorp.com/ubuntu/bionic64" config.vm.define "nansat", primary: true do |nansat| end From 1a62f3dcf32f6484a4cdda386ea37218a28b25a7 Mon Sep 17 00:00:00 2001 From: "Morten W. Hansen" Date: Mon, 3 Jul 2023 14:00:24 +0000 Subject: [PATCH 23/57] some error handling --- nansat/mappers/mapper_netcdf_cf.py | 8 +++----- nansat/mappers/mapper_opendap_arome.py | 4 ++-- 2 files changed, 5 insertions(+), 7 deletions(-) diff --git a/nansat/mappers/mapper_netcdf_cf.py b/nansat/mappers/mapper_netcdf_cf.py index 4227bfb98..c96061143 100644 --- a/nansat/mappers/mapper_netcdf_cf.py +++ b/nansat/mappers/mapper_netcdf_cf.py @@ -443,7 +443,7 @@ def _band_dict(self, subfilename, band_num, subds, band=None, band_metadata=None def _create_empty(self, gdal_dataset, gdal_metadata): try: self._create_empty_from_projection_variable(gdal_dataset, gdal_metadata) - except KeyError: + except (KeyError, AttributeError): try: self._create_empty_from_subdatasets(gdal_dataset, gdal_metadata) except NansatMissingProjectionError: @@ -464,10 +464,8 @@ def _create_empty_from_projection_variable(self, gdal_dataset, gdal_metadata): for var in ds.variables.keys(): if ds[var].shape==tuple(shape): break - try: - grid_mapping = ds.variables[ds.variables[var].grid_mapping] - except KeyError: - raise WrongMapperError + + grid_mapping = ds.variables[ds.variables[var].grid_mapping] grid_mapping_dict = {} for index in grid_mapping.ncattrs(): diff --git a/nansat/mappers/mapper_opendap_arome.py b/nansat/mappers/mapper_opendap_arome.py index 8575fb63e..056f70443 100644 --- a/nansat/mappers/mapper_opendap_arome.py +++ b/nansat/mappers/mapper_opendap_arome.py @@ -53,8 +53,8 @@ def __init__(self, filename, gdal_dataset, gdal_metadata, netcdf_dim=None, for attr in ds.ncattrs(): metadata[attr] = self._fix_encoding(ds.getncattr(attr)) - if not 'arome' in metadata['title'].lower() and \ - not 'meps' in metadata['title'].lower(): + if 'title' not in metadata.keys() or ('arome' not in metadata['title'].lower() and 'meps' \ + not in metadata['title'].lower()): raise WrongMapperError xsize = ds.dimensions['x'].size From 3f488fb5758b1bb2e6aed32c555975cb88ff0d64 Mon Sep 17 00:00:00 2001 From: "Morten W. Hansen" Date: Tue, 27 Feb 2024 22:14:35 +0100 Subject: [PATCH 24/57] start meps mapper --- nansat/mappers/mapper_meps.py | 84 +++++++++++++++++++++++++++++++++++ 1 file changed, 84 insertions(+) create mode 100644 nansat/mappers/mapper_meps.py diff --git a/nansat/mappers/mapper_meps.py b/nansat/mappers/mapper_meps.py new file mode 100644 index 000000000..96da492c5 --- /dev/null +++ b/nansat/mappers/mapper_meps.py @@ -0,0 +1,84 @@ +from pyproj import CRS +from netCDF4 import Dataset + +from nansat.nsr import NSR + +import requests + +import xml.etree.ElementTree as ET + +from dateutil.parser import parse + +import json +import pythesint as pti + +from nansat.mappers.mapper_netcdf_cf import Mapper as NetcdfCF +from nansat.exceptions import WrongMapperError + +class Mapper(NetcdfCF): + + def __init__(self, filename, gdal_dataset, gdal_metadata, file_num=0, *args, **kwargs): + + if not filename.endswith(".ncml"): + raise WrongMapperError + + ds = Dataset(filename) + for attr in ds.ncattrs(): + gdal_metadata[attr] = ds.getncattr(attr) + + if not 'meps' in gdal_metadata['title'].lower(): + raise WrongMapperError + + self.input_filename = filename + + self._create_empty_from_projection_variable(gdal_dataset, gdal_metadata) + + import ipdb + ipdb.set_trace() + # Add bands with metadata and corresponding values to the empty VRT + self.create_bands(self._band_list(gdal_dataset, metadata, *args, **kwargs)) + + #self.dataset.SetMetadataItem('time_coverage_start', parse( + # gdal_metadata['NC_GLOBAL#min_time'], ignoretz=True, fuzzy=True).isoformat()) + #self.dataset.SetMetadataItem('time_coverage_end', parse( + # gdal_metadata['NC_GLOBAL#max_time'], ignoretz=True, fuzzy=True).isoformat())) + + # Get dictionary describing the instrument and platform according to + # the GCMD keywords + 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)) + + + def _create_empty_from_projection_variable(self, gdal_dataset, gdal_metadata): + + # TODO: only open dataset once or close them.. + ds = Dataset(self.input_filename) + + shape=[] + for key in ds.dimensions.keys(): + shape.append(ds.dimensions[key].size) + + for var in ds.variables.keys(): + if ds[var].shape==tuple(shape): + break + + grid_mapping = ds.variables[ds.variables[var].grid_mapping] + + grid_mapping_dict = {} + for index in grid_mapping.ncattrs(): + grid_mapping_dict[index] = grid_mapping.getncattr(index) + crs = CRS.from_cf(grid_mapping_dict) + + xx = ds.variables["x"][0:2] + yy = ds.variables["y"][0:2] + geotr = (xx[0], xx[1]-xx[0], 0, yy[0], 0, yy[1]-yy[0]) + + self._init_from_dataset_params( + x_size = ds.dimensions["x"].size, + y_size = ds.dimensions["y"].size, + geo_transform = geotr, + projection = NSR(crs.to_proj4()).wkt, + metadata = gdal_metadata) From 0e8447fc4a951b592234648a4254df75e6146735 Mon Sep 17 00:00:00 2001 From: "Morten W. Hansen" Date: Wed, 28 Feb 2024 11:49:35 +0100 Subject: [PATCH 25/57] Should work now but does not.. --- nansat/mappers/mapper_meps.py | 35 +++++++++++++++-------------------- 1 file changed, 15 insertions(+), 20 deletions(-) diff --git a/nansat/mappers/mapper_meps.py b/nansat/mappers/mapper_meps.py index 96da492c5..05c281025 100644 --- a/nansat/mappers/mapper_meps.py +++ b/nansat/mappers/mapper_meps.py @@ -1,17 +1,13 @@ +import os +import json + +import pythesint as pti + +from osgeo import gdal from pyproj import CRS from netCDF4 import Dataset from nansat.nsr import NSR - -import requests - -import xml.etree.ElementTree as ET - -from dateutil.parser import parse - -import json -import pythesint as pti - from nansat.mappers.mapper_netcdf_cf import Mapper as NetcdfCF from nansat.exceptions import WrongMapperError @@ -29,19 +25,13 @@ def __init__(self, filename, gdal_dataset, gdal_metadata, file_num=0, *args, **k if not 'meps' in gdal_metadata['title'].lower(): raise WrongMapperError - self.input_filename = filename - - self._create_empty_from_projection_variable(gdal_dataset, gdal_metadata) - + url = self._get_odap_url(filename, file_num) + gdal_dataset = gdal.Open(url) + metadata = gdal_dataset.GetMetadata() import ipdb ipdb.set_trace() - # Add bands with metadata and corresponding values to the empty VRT - self.create_bands(self._band_list(gdal_dataset, metadata, *args, **kwargs)) - #self.dataset.SetMetadataItem('time_coverage_start', parse( - # gdal_metadata['NC_GLOBAL#min_time'], ignoretz=True, fuzzy=True).isoformat()) - #self.dataset.SetMetadataItem('time_coverage_end', parse( - # gdal_metadata['NC_GLOBAL#max_time'], ignoretz=True, fuzzy=True).isoformat())) + super(Mapper, self).__init__(url, gdal_dataset, metadata, *args, **kwargs) # Get dictionary describing the instrument and platform according to # the GCMD keywords @@ -52,6 +42,11 @@ def __init__(self, filename, gdal_dataset, gdal_metadata, file_num=0, *args, **k self.dataset.SetMetadataItem('platform', json.dumps(ee)) + def _get_odap_url(self, fn, file_num): + url = os.path.split(fn)[0] + "/member_%02d" % int(os.path.basename(fn)[8:11]) + "/meps_" + os.path.basename(fn).split("_")[2] + "_%02d_" % file_num + os.path.basename(fn).split("_")[3][:-2] + return url + + def _create_empty_from_projection_variable(self, gdal_dataset, gdal_metadata): # TODO: only open dataset once or close them.. From 4e71e74509892742cfcf83c44b43ed7e007f68b3 Mon Sep 17 00:00:00 2001 From: "Morten W. Hansen" Date: Wed, 28 Feb 2024 18:03:24 +0100 Subject: [PATCH 26/57] update --- nansat/mappers/mapper_arome.py | 3 +-- nansat/mappers/{mapper_meps.py => mapper_meps_ncml.py} | 2 +- 2 files changed, 2 insertions(+), 3 deletions(-) rename nansat/mappers/{mapper_meps.py => mapper_meps_ncml.py} (96%) diff --git a/nansat/mappers/mapper_arome.py b/nansat/mappers/mapper_arome.py index c7e4af60e..e02cc728f 100644 --- a/nansat/mappers/mapper_arome.py +++ b/nansat/mappers/mapper_arome.py @@ -14,8 +14,7 @@ def __init__(self, filename, gdal_dataset, gdal_metadata, *args, **kwargs): raise WrongMapperError if 'source' not in list(gdal_metadata.keys()): raise WrongMapperError - if not 'arome' in gdal_metadata['title'].lower() and \ - not 'meps' in gdal_metadata['title'].lower(): + if not 'arome' in gdal_metadata['title'].lower(): raise WrongMapperError super(Mapper, self).__init__(filename, gdal_dataset, gdal_metadata, *args, **kwargs) diff --git a/nansat/mappers/mapper_meps.py b/nansat/mappers/mapper_meps_ncml.py similarity index 96% rename from nansat/mappers/mapper_meps.py rename to nansat/mappers/mapper_meps_ncml.py index 05c281025..9cb6c1369 100644 --- a/nansat/mappers/mapper_meps.py +++ b/nansat/mappers/mapper_meps_ncml.py @@ -31,7 +31,7 @@ def __init__(self, filename, gdal_dataset, gdal_metadata, file_num=0, *args, **k import ipdb ipdb.set_trace() - super(Mapper, self).__init__(url, gdal_dataset, metadata, *args, **kwargs) + super(Mapper, self).__init__(url, gdal_dataset, gdal_metadata, *args, **kwargs) # Get dictionary describing the instrument and platform according to # the GCMD keywords From 6cdea78ea4c97762bdda3d3666dd6a7c8a08d804 Mon Sep 17 00:00:00 2001 From: "Morten W. Hansen" Date: Wed, 28 Feb 2024 18:51:03 +0100 Subject: [PATCH 27/57] reorganise --- nansat/mappers/mapper_meps.py | 36 +++++++++++++++ nansat/mappers/mapper_meps_ncml.py | 71 ++++-------------------------- 2 files changed, 44 insertions(+), 63 deletions(-) create mode 100644 nansat/mappers/mapper_meps.py diff --git a/nansat/mappers/mapper_meps.py b/nansat/mappers/mapper_meps.py new file mode 100644 index 000000000..9c7e00a5e --- /dev/null +++ b/nansat/mappers/mapper_meps.py @@ -0,0 +1,36 @@ +import json +import pythesint as pti + +from osgeo import gdal +from netCDF4 import Dataset + +from nansat.exceptions import WrongMapperError +from nansat.mappers.mapper_netcdf_cf import Mapper as NetcdfCF + + +class Mapper(NetcdfCF): + + def __init__(self, url, gdal_dataset, gdal_metadata, file_num=0, *args, **kwargs): + + if not url.endswith(".nc"): + raise WrongMapperError + + ds = Dataset(url) + for attr in ds.ncattrs(): + gdal_metadata[attr] = ds.getncattr(attr) + + if not 'meps' in gdal_metadata['title'].lower(): + raise WrongMapperError + + gdal_dataset = gdal.Open(url) + metadata = gdal_dataset.GetMetadata() + + super(Mapper, self).__init__(url, gdal_dataset, gdal_metadata, *args, **kwargs) + + # Get dictionary describing the instrument and platform according to + # the GCMD keywords + 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)) diff --git a/nansat/mappers/mapper_meps_ncml.py b/nansat/mappers/mapper_meps_ncml.py index 9cb6c1369..95bcea9aa 100644 --- a/nansat/mappers/mapper_meps_ncml.py +++ b/nansat/mappers/mapper_meps_ncml.py @@ -1,79 +1,24 @@ import os -import json -import pythesint as pti - -from osgeo import gdal -from pyproj import CRS -from netCDF4 import Dataset - -from nansat.nsr import NSR -from nansat.mappers.mapper_netcdf_cf import Mapper as NetcdfCF +from nansat.mappers.mapper_meps import Mapper as NCMapper from nansat.exceptions import WrongMapperError -class Mapper(NetcdfCF): +class Mapper(NCMapper): - def __init__(self, filename, gdal_dataset, gdal_metadata, file_num=0, *args, **kwargs): - - if not filename.endswith(".ncml"): - raise WrongMapperError - ds = Dataset(filename) - for attr in ds.ncattrs(): - gdal_metadata[attr] = ds.getncattr(attr) + def __init__(self, ncml_url, gdal_dataset, gdal_metadata, file_num=0, *args, **kwargs): - if not 'meps' in gdal_metadata['title'].lower(): + if not ncml_url.endswith(".ncml"): raise WrongMapperError - url = self._get_odap_url(filename, file_num) - gdal_dataset = gdal.Open(url) - metadata = gdal_dataset.GetMetadata() - import ipdb - ipdb.set_trace() + url = self._get_odap_url(ncml_url, file_num) + # Dette fungerer ikke fordi gdal ikke finner gds.GetSubDatasets() = [] + # Men ncml-filen inneholder lustre-stier, så vi kan evt bruke de direkte, + # og lese som netcdf på ppi - men dette blir en workaround for MET. Opendap burde fungere... super(Mapper, self).__init__(url, gdal_dataset, gdal_metadata, *args, **kwargs) - # Get dictionary describing the instrument and platform according to - # the GCMD keywords - 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)) - def _get_odap_url(self, fn, file_num): url = os.path.split(fn)[0] + "/member_%02d" % int(os.path.basename(fn)[8:11]) + "/meps_" + os.path.basename(fn).split("_")[2] + "_%02d_" % file_num + os.path.basename(fn).split("_")[3][:-2] return url - - - def _create_empty_from_projection_variable(self, gdal_dataset, gdal_metadata): - - # TODO: only open dataset once or close them.. - ds = Dataset(self.input_filename) - - shape=[] - for key in ds.dimensions.keys(): - shape.append(ds.dimensions[key].size) - - for var in ds.variables.keys(): - if ds[var].shape==tuple(shape): - break - - grid_mapping = ds.variables[ds.variables[var].grid_mapping] - - grid_mapping_dict = {} - for index in grid_mapping.ncattrs(): - grid_mapping_dict[index] = grid_mapping.getncattr(index) - crs = CRS.from_cf(grid_mapping_dict) - - xx = ds.variables["x"][0:2] - yy = ds.variables["y"][0:2] - geotr = (xx[0], xx[1]-xx[0], 0, yy[0], 0, yy[1]-yy[0]) - - self._init_from_dataset_params( - x_size = ds.dimensions["x"].size, - y_size = ds.dimensions["y"].size, - geo_transform = geotr, - projection = NSR(crs.to_proj4()).wkt, - metadata = gdal_metadata) From 90e09f9ce62aa4dd81514cc18d91ed6d263709b5 Mon Sep 17 00:00:00 2001 From: "Morten W. Hansen" Date: Fri, 8 Mar 2024 16:40:17 +0100 Subject: [PATCH 28/57] MEPS mapper works --- nansat/mappers/mapper_meps.py | 76 +++++++++++++++++++++++--- nansat/mappers/mapper_meps_ncml.py | 14 +++-- nansat/mappers/mapper_netcdf_cf.py | 2 - nansat/mappers/mapper_opendap_arome.py | 60 ++------------------ nansat/mappers/opendap.py | 36 ++++++++++++ 5 files changed, 117 insertions(+), 71 deletions(-) diff --git a/nansat/mappers/mapper_meps.py b/nansat/mappers/mapper_meps.py index 9c7e00a5e..841ae1787 100644 --- a/nansat/mappers/mapper_meps.py +++ b/nansat/mappers/mapper_meps.py @@ -2,30 +2,90 @@ import pythesint as pti from osgeo import gdal +from pyproj import CRS from netCDF4 import Dataset +from nansat.nsr import NSR from nansat.exceptions import WrongMapperError from nansat.mappers.mapper_netcdf_cf import Mapper as NetcdfCF +from nansat.mappers.opendap import Opendap -class Mapper(NetcdfCF): +class Mapper(NetcdfCF, Opendap): - def __init__(self, url, gdal_dataset, gdal_metadata, file_num=0, *args, **kwargs): + def __init__(self, url, gdal_dataset, gdal_metadata, file_num=0, bands=None, *args, **kwargs): if not url.endswith(".nc"): raise WrongMapperError - ds = Dataset(url) + try: + ds = Dataset(url) + except OSError: + raise WrongMapperError + + metadata = {} for attr in ds.ncattrs(): - gdal_metadata[attr] = ds.getncattr(attr) + metadata[attr] = ds.getncattr(attr) + + if 'title' not in metadata.keys() or 'meps' not in metadata['title'].lower(): + raise WrongMapperError - if not 'meps' in gdal_metadata['title'].lower(): + self.input_filename = url + + xsize = ds.dimensions['x'].size + ysize = ds.dimensions['y'].size + + # Pick 10 meter height dimension only + height_dim = 'height6' + if height_dim not in ds.dimensions.keys(): + raise WrongMapperError + if ds.dimensions[height_dim].size != 1: + raise WrongMapperError + if ds.variables[height_dim][0].data != 10: raise WrongMapperError - gdal_dataset = gdal.Open(url) - metadata = gdal_dataset.GetMetadata() + varnames = [] + for var in ds.variables: + var_dimensions = ds.variables[var].dimensions + if var_dimensions == ('time', height_dim, '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, {}, dim_sizes) + fn = self._get_sub_filename(url, band, dim_sizes, index) + meta_dict.append(self.get_band_metadata_dict(fn, ds.variables[band])) + + self.create_bands(meta_dict) - super(Mapper, self).__init__(url, gdal_dataset, gdal_metadata, *args, **kwargs) + # Copy metadata + for key in metadata.keys(): + self.dataset.SetMetadataItem(str(key), str(metadata[key])) # Get dictionary describing the instrument and platform according to # the GCMD keywords diff --git a/nansat/mappers/mapper_meps_ncml.py b/nansat/mappers/mapper_meps_ncml.py index 95bcea9aa..4fa8602bf 100644 --- a/nansat/mappers/mapper_meps_ncml.py +++ b/nansat/mappers/mapper_meps_ncml.py @@ -12,13 +12,17 @@ def __init__(self, ncml_url, gdal_dataset, gdal_metadata, file_num=0, *args, **k raise WrongMapperError url = self._get_odap_url(ncml_url, file_num) - # Dette fungerer ikke fordi gdal ikke finner gds.GetSubDatasets() = [] - # Men ncml-filen inneholder lustre-stier, så vi kan evt bruke de direkte, - # og lese som netcdf på ppi - men dette blir en workaround for MET. Opendap burde fungere... super(Mapper, self).__init__(url, gdal_dataset, gdal_metadata, *args, **kwargs) - def _get_odap_url(self, fn, file_num): - url = os.path.split(fn)[0] + "/member_%02d" % int(os.path.basename(fn)[8:11]) + "/meps_" + os.path.basename(fn).split("_")[2] + "_%02d_" % file_num + os.path.basename(fn).split("_")[3][:-2] + def _get_odap_url(self, fn, file_num=0): + """ Get the opendap url to file number 'file_num'. The + default file number is 0, and yields the forecast time. + """ + url = ( + "" + os.path.split(fn)[0] + "/member_%02d" + "/meps_" + os.path.basename(fn).split("_")[2] + + "_%02d_" + os.path.basename(fn).split("_")[3][:-2] + ) % (int(os.path.basename(fn)[8:11]), file_num) return url diff --git a/nansat/mappers/mapper_netcdf_cf.py b/nansat/mappers/mapper_netcdf_cf.py index c96061143..e481209bb 100644 --- a/nansat/mappers/mapper_netcdf_cf.py +++ b/nansat/mappers/mapper_netcdf_cf.py @@ -261,8 +261,6 @@ def _get_index_of_dimensions(self, dimension_names, netcdf_dim, dim_sizes): # Get band number from given timestamp index = int(np.argmin(np.abs(self.times() - val))) else: - import pdb - pdb.set_trace() index = int(np.argmin(np.abs(ds.variables[key][:] - val))) index4key[key] = { 'index': index, diff --git a/nansat/mappers/mapper_opendap_arome.py b/nansat/mappers/mapper_opendap_arome.py index 056f70443..2175685f8 100644 --- a/nansat/mappers/mapper_opendap_arome.py +++ b/nansat/mappers/mapper_opendap_arome.py @@ -53,8 +53,7 @@ def __init__(self, filename, gdal_dataset, gdal_metadata, netcdf_dim=None, for attr in ds.ncattrs(): metadata[attr] = self._fix_encoding(ds.getncattr(attr)) - if 'title' not in metadata.keys() or ('arome' not in metadata['title'].lower() and 'meps' \ - not in metadata['title'].lower()): + if 'title' not in metadata.keys() or 'arome' not in metadata['title'].lower(): raise WrongMapperError xsize = ds.dimensions['x'].size @@ -66,12 +65,13 @@ def __init__(self, filename, gdal_dataset, gdal_metadata, netcdf_dim=None, raise WrongMapperError if ds.dimensions[height_dim].size != 1: raise WrongMapperError - if ds.variables['height7'][0].data != 10: + if ds.variables[height_dim][0].data != 10: raise WrongMapperError + varnames = [] for var in ds.variables: var_dimensions = ds.variables[var].dimensions - if var_dimensions == ('time', 'height7', 'y', 'x'): + if var_dimensions == ('time', height_dim, 'y', 'x'): varnames.append(var) # Projection @@ -116,58 +116,6 @@ def __init__(self, filename, gdal_dataset, gdal_metadata, netcdf_dim=None, self.dataset.SetMetadataItem('instrument', json.dumps(mm)) self.dataset.SetMetadataItem('platform', json.dumps(ee)) - #md_item = 'Data Center' - #if not self.dataset.GetMetadataItem(md_item): - # self.dataset.SetMetadataItem(md_item, 'NO/MET') - #md_item = 'Entry Title' - #if not self.dataset.GetMetadataItem(md_item): - # self.dataset.SetMetadataItem(md_item, str(ds.getncattr('title'))) - #md_item = 'summary' - #if not self.dataset.GetMetadataItem(md_item): - # summary = """ - # AROME_Arctic is a convection-permitting atmosphere model covering parts of the Barents - # Sea and the Nordic Arctic. It has horizontal resolution of 2.5 km and 65 vertical - # levels. AROME_Arctic runs for 66 hours four times a day (00,06,12,18) with three-hourly - # cycling for data assimilation. Boundary data is from ECMWF. Model code based on HARMONIE - # cy40h1.1 - # """ - # self.dataset.SetMetadataItem(md_item, str(summary)) - - def get_band_metadata_dict(self, fn, ncvar): - gds = gdal.Open(fn) - meta_item = { - 'src': {'SourceFilename': fn, 'SourceBand': 1}, - 'dst': {'name': ncvar.name, 'dataType': 6} - } - - for attr_key in ncvar.ncattrs(): - attr_val = self._fix_encoding(ncvar.getncattr(attr_key)) - if attr_key in ['scale', 'scale_factor']: - meta_item['src']['ScaleRatio'] = attr_val - elif attr_key in ['offset', 'add_offset']: - meta_item['src']['ScaleOffset'] = attr_val - else: - meta_item['dst'][attr_key] = attr_val - - return meta_item - - - - @staticmethod - def _get_sub_filename(url, var, dim_sizes, index): - """ Opendap driver refers to subdatasets differently than the - standard way in vrt.py - """ - shape = [] - for item in dim_sizes.items(): - if item[0] in index.keys(): - shape.append(index[item[0]]['index']) - else: - shape.append(item[0]) - # assemble dimensions string - gd_shape = ''.join(['[%s]' % dimsize for dimsize in shape]) - return '{url}?{var}.{var}{shape}'.format(url=url, var=var, shape=gd_shape) - @staticmethod def get_date(filename): """Extract date and time parameters from filename and return diff --git a/nansat/mappers/opendap.py b/nansat/mappers/opendap.py index 418c21868..83bcb812b 100644 --- a/nansat/mappers/opendap.py +++ b/nansat/mappers/opendap.py @@ -217,6 +217,42 @@ def _fix_encoding(var): retval = var return retval + + def get_band_metadata_dict(self, fn, ncvar): + gds = gdal.Open(fn) + meta_item = { + 'src': {'SourceFilename': fn, 'SourceBand': 1}, + 'dst': {'name': ncvar.name, 'dataType': 6} + } + + for attr_key in ncvar.ncattrs(): + attr_val = self._fix_encoding(ncvar.getncattr(attr_key)) + if attr_key in ['scale', 'scale_factor']: + meta_item['src']['ScaleRatio'] = attr_val + elif attr_key in ['offset', 'add_offset']: + meta_item['src']['ScaleOffset'] = attr_val + else: + meta_item['dst'][attr_key] = attr_val + + return meta_item + + + @staticmethod + def _get_sub_filename(url, var, dim_sizes, index): + """ Opendap driver refers to subdatasets differently than the + standard way in vrt.py + """ + shape = [] + for item in dim_sizes.items(): + if item[0] in index.keys(): + shape.append(index[item[0]]['index']) + else: + shape.append(item[0]) + # assemble dimensions string + gd_shape = ''.join(['[%s]' % dimsize for dimsize in shape]) + return '{url}?{var}.{var}{shape}'.format(url=url, var=var, shape=gd_shape) + + def create_vrt(self, filename, gdalDataset, gdalMetadata, date, ds, bands, cachedir): """ Create VRT From 2e7dedf42e5d565fa89b4a7dab50641a0c195762 Mon Sep 17 00:00:00 2001 From: "Morten W. Hansen" Date: Tue, 9 Apr 2024 12:28:33 +0200 Subject: [PATCH 29/57] Add MET Nordic mapper --- nansat/mappers/mapper_opendap_arome.py | 4 ---- 1 file changed, 4 deletions(-) diff --git a/nansat/mappers/mapper_opendap_arome.py b/nansat/mappers/mapper_opendap_arome.py index 2175685f8..69f9eef17 100644 --- a/nansat/mappers/mapper_opendap_arome.py +++ b/nansat/mappers/mapper_opendap_arome.py @@ -1,7 +1,3 @@ -# 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 From b26dedb0b7a40543dd685c552c3c5fe53635c21f Mon Sep 17 00:00:00 2001 From: Morten Wergeland Hansen Date: Mon, 15 Apr 2024 18:11:44 +0200 Subject: [PATCH 30/57] add source filename to metadata --- nansat/mappers/mapper_meps.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/nansat/mappers/mapper_meps.py b/nansat/mappers/mapper_meps.py index 841ae1787..a8b635b61 100644 --- a/nansat/mappers/mapper_meps.py +++ b/nansat/mappers/mapper_meps.py @@ -94,3 +94,6 @@ def __init__(self, url, gdal_dataset, gdal_metadata, file_num=0, bands=None, *ar 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) From f4306c556753edf0ddcde2ec050ab20f9d8e6d18 Mon Sep 17 00:00:00 2001 From: "Morten W. Hansen" Date: Wed, 17 Apr 2024 08:03:08 +0200 Subject: [PATCH 31/57] Now add MET Nordic mapper --- nansat/mappers/mapper_meps.py | 6 +- nansat/mappers/mapper_opendap_met_nordic.py | 116 ++++++++++++++++++++ 2 files changed, 121 insertions(+), 1 deletion(-) create mode 100644 nansat/mappers/mapper_opendap_met_nordic.py diff --git a/nansat/mappers/mapper_meps.py b/nansat/mappers/mapper_meps.py index 841ae1787..21c597416 100644 --- a/nansat/mappers/mapper_meps.py +++ b/nansat/mappers/mapper_meps.py @@ -25,7 +25,11 @@ def __init__(self, url, gdal_dataset, gdal_metadata, file_num=0, bands=None, *ar metadata = {} for attr in ds.ncattrs(): - metadata[attr] = ds.getncattr(attr) + content = ds.getncattr(attr) + content.replace("æ", "ae") + content.replace("ø", "oe") + content.replace("å", "aa") + metadata[attr] = content if 'title' not in metadata.keys() or 'meps' not in metadata['title'].lower(): raise WrongMapperError diff --git a/nansat/mappers/mapper_opendap_met_nordic.py b/nansat/mappers/mapper_opendap_met_nordic.py new file mode 100644 index 000000000..d17b2b967 --- /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 + From f43e9b573bb46090d3113f90c2be9432d6de3033 Mon Sep 17 00:00:00 2001 From: "Morten W. Hansen" Date: Wed, 17 Apr 2024 08:23:29 +0200 Subject: [PATCH 32/57] fix encoding of metadata --- nansat/mappers/mapper_meps.py | 9 ++++----- 1 file changed, 4 insertions(+), 5 deletions(-) diff --git a/nansat/mappers/mapper_meps.py b/nansat/mappers/mapper_meps.py index 7486b07f9..f1aed9f6b 100644 --- a/nansat/mappers/mapper_meps.py +++ b/nansat/mappers/mapper_meps.py @@ -13,7 +13,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 @@ -25,10 +26,8 @@ def __init__(self, url, gdal_dataset, gdal_metadata, file_num=0, bands=None, *ar metadata = {} for attr in ds.ncattrs(): - content = ds.getncattr(attr) - content.replace("æ", "ae") - content.replace("ø", "oe") - content.replace("å", "aa") + content = ds.getncattr(attr).replace("æ", "ae").replace("ø", "oe").replace("å", + "aa") metadata[attr] = content if 'title' not in metadata.keys() or 'meps' not in metadata['title'].lower(): From 251dba123603a4f387820646e4a1177d05f46f9b Mon Sep 17 00:00:00 2001 From: Morten Wergeland Hansen Date: Wed, 17 Apr 2024 14:15:39 +0200 Subject: [PATCH 33/57] reorganise and change apostrophes --- nansat/mappers/mapper_meps.py | 32 ++++++++++++++++---------------- 1 file changed, 16 insertions(+), 16 deletions(-) diff --git a/nansat/mappers/mapper_meps.py b/nansat/mappers/mapper_meps.py index f1aed9f6b..4c1e3463a 100644 --- a/nansat/mappers/mapper_meps.py +++ b/nansat/mappers/mapper_meps.py @@ -24,22 +24,22 @@ 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(): - content = ds.getncattr(attr).replace("æ", "ae").replace("ø", "oe").replace("å", - "aa") + content = ds.getncattr(attr) + content = content.replace("æ", "ae").replace("ø", "oe").replace("å", "aa") metadata[attr] = content - if 'title' not in metadata.keys() or 'meps' not in metadata['title'].lower(): - raise WrongMapperError - 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: @@ -50,7 +50,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 @@ -66,8 +66,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) @@ -92,11 +92,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) + self.dataset.SetMetadataItem("nc_file", self.input_filename) From a9556bb4dcdf1258b5a49fa85962e76a9346192c Mon Sep 17 00:00:00 2001 From: Morten Wergeland Hansen Date: Wed, 17 Apr 2024 16:01:50 +0200 Subject: [PATCH 34/57] Handle non-string attributes --- nansat/mappers/mapper_meps.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/nansat/mappers/mapper_meps.py b/nansat/mappers/mapper_meps.py index 4c1e3463a..399cea557 100644 --- a/nansat/mappers/mapper_meps.py +++ b/nansat/mappers/mapper_meps.py @@ -30,7 +30,8 @@ def __init__(self, url, gdal_dataset, gdal_metadata, file_num=0, bands=None, *ar metadata = {} for attr in ds.ncattrs(): content = ds.getncattr(attr) - content = content.replace("æ", "ae").replace("ø", "oe").replace("å", "aa") + if isinstance(content, str): + content = content.replace("æ", "ae").replace("ø", "oe").replace("å", "aa") metadata[attr] = content self.input_filename = url From cd929047e681f59b6eb4c162dafc1df6baa481f9 Mon Sep 17 00:00:00 2001 From: "Morten W. Hansen" Date: Mon, 29 Apr 2024 10:11:48 +0200 Subject: [PATCH 35/57] Remove warnings --- nansat/mappers/opendap.py | 2 -- 1 file changed, 2 deletions(-) diff --git a/nansat/mappers/opendap.py b/nansat/mappers/opendap.py index 83bcb812b..c0d604ecb 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}) From 731e5fb5b609d127fa5f8328042bc752d3beb72b Mon Sep 17 00:00:00 2001 From: "Morten W. Hansen" Date: Tue, 30 Apr 2024 16:51:43 +0200 Subject: [PATCH 36/57] Make it possible to avoid creating gcps --- nansat/exporter.py | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/nansat/exporter.py b/nansat/exporter.py index 2844e14d1..53e0048d1 100644 --- a/nansat/exporter.py +++ b/nansat/exporter.py @@ -46,7 +46,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 +112,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, From f410d5b476d269d4fd8dad6746ab1190b48fd3e4 Mon Sep 17 00:00:00 2001 From: "Morten W. Hansen" Date: Tue, 30 Apr 2024 20:42:08 +0200 Subject: [PATCH 37/57] Add time metadata to bands --- nansat/mappers/mapper_meps.py | 11 ++++++++++- nansat/mappers/mapper_opendap_arome.py | 16 ++++++++++++---- 2 files changed, 22 insertions(+), 5 deletions(-) diff --git a/nansat/mappers/mapper_meps.py b/nansat/mappers/mapper_meps.py index 399cea557..c756315c2 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 @@ -83,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) diff --git a/nansat/mappers/mapper_opendap_arome.py b/nansat/mappers/mapper_opendap_arome.py index 69f9eef17..6f45fdeba 100644 --- a/nansat/mappers/mapper_opendap_arome.py +++ b/nansat/mappers/mapper_opendap_arome.py @@ -1,12 +1,14 @@ # 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 @@ -99,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) @@ -118,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""" From 913c7aaa6ce126b3794c9b94fbd9e5aa01e6c4e6 Mon Sep 17 00:00:00 2001 From: Morten Wergeland Hansen Date: Fri, 3 May 2024 17:17:34 +0200 Subject: [PATCH 38/57] allow missing history, and set timezone to utc --- nansat/exporter.py | 9 +++++++-- 1 file changed, 7 insertions(+), 2 deletions(-) diff --git a/nansat/exporter.py b/nansat/exporter.py index 53e0048d1..25c50a872 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 @@ -130,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!') @@ -368,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) From 7097cfddc287661022cb00dc1ba3c2e57d4649ee Mon Sep 17 00:00:00 2001 From: "Morten W. Hansen" Date: Sat, 4 May 2024 15:01:44 +0200 Subject: [PATCH 39/57] find correct nc-file based on input time --- nansat/mappers/mapper_meps_ncml.py | 14 +++++++++++--- 1 file changed, 11 insertions(+), 3 deletions(-) diff --git a/nansat/mappers/mapper_meps_ncml.py b/nansat/mappers/mapper_meps_ncml.py index 4fa8602bf..686b081e9 100644 --- a/nansat/mappers/mapper_meps_ncml.py +++ b/nansat/mappers/mapper_meps_ncml.py @@ -1,17 +1,25 @@ import os +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) + ds = netCDF4.Dataset(ncml_url) + time = netcdf_dim["time"] + dt = time - np.datetime64( + datetime.datetime.fromisoformat(ds.time_coverage_start.replace("Z", "+00:00"))) + url = self._get_odap_url(ncml_url, np.round(dt)) super(Mapper, self).__init__(url, gdal_dataset, gdal_metadata, *args, **kwargs) From e222105a6ecf8766d418cb754e645c0b75087db2 Mon Sep 17 00:00:00 2001 From: Morten Wergeland Hansen Date: Mon, 13 May 2024 16:02:17 +0200 Subject: [PATCH 40/57] Assert correct gcp shape --- nansat/mappers/sentinel1.py | 32 +++++++++++++++++++++++--------- 1 file changed, 23 insertions(+), 9 deletions(-) diff --git a/nansat/mappers/sentinel1.py b/nansat/mappers/sentinel1.py index 4c14732ee..d763e2809 100644 --- a/nansat/mappers/sentinel1.py +++ b/nansat/mappers/sentinel1.py @@ -107,13 +107,30 @@ 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. + """ + if gcp_y*gcp_x != pixel.size: + gcp_y = gcp_y - pixel.size % gcp_y + gcp_x = gcp_x - pixel.size % gcp_x + assert gcp_y*gcp_x == pixel.size + + return gcp_y, 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 +145,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): From a3457f17e6430d3053c97fca4e8195ec7598e970 Mon Sep 17 00:00:00 2001 From: Morten Wergeland Hansen Date: Mon, 13 May 2024 17:16:16 +0200 Subject: [PATCH 41/57] better handle input dict --- nansat/mappers/mapper_meps_ncml.py | 12 ++++++++---- 1 file changed, 8 insertions(+), 4 deletions(-) diff --git a/nansat/mappers/mapper_meps_ncml.py b/nansat/mappers/mapper_meps_ncml.py index 686b081e9..d5f878c46 100644 --- a/nansat/mappers/mapper_meps_ncml.py +++ b/nansat/mappers/mapper_meps_ncml.py @@ -15,10 +15,14 @@ def __init__(self, ncml_url, gdal_dataset, gdal_metadata, netcdf_dim=None, *args if not ncml_url.endswith(".ncml"): raise WrongMapperError - ds = netCDF4.Dataset(ncml_url) - time = netcdf_dim["time"] - dt = time - np.datetime64( - datetime.datetime.fromisoformat(ds.time_coverage_start.replace("Z", "+00:00"))) + dt = 0 + if netcdf_dim is not None and "time" in netcdf_dim.keys(): + ds = netCDF4.Dataset(ncml_url) + time = netcdf_dim["time"] + dt = time - np.datetime64( + datetime.datetime.fromisoformat(ds.time_coverage_start.replace("Z", "+00:00"))) + import ipdb + ipdb.set_trace() url = self._get_odap_url(ncml_url, np.round(dt)) super(Mapper, self).__init__(url, gdal_dataset, gdal_metadata, *args, **kwargs) From cdb110daa385640e9f6fe14789640945017cd35d Mon Sep 17 00:00:00 2001 From: "Morten W. Hansen" Date: Mon, 13 May 2024 20:09:14 +0200 Subject: [PATCH 42/57] remove ipdb lines --- nansat/mappers/mapper_meps_ncml.py | 2 -- 1 file changed, 2 deletions(-) diff --git a/nansat/mappers/mapper_meps_ncml.py b/nansat/mappers/mapper_meps_ncml.py index d5f878c46..44ad98bbb 100644 --- a/nansat/mappers/mapper_meps_ncml.py +++ b/nansat/mappers/mapper_meps_ncml.py @@ -21,8 +21,6 @@ def __init__(self, ncml_url, gdal_dataset, gdal_metadata, netcdf_dim=None, *args time = netcdf_dim["time"] dt = time - np.datetime64( datetime.datetime.fromisoformat(ds.time_coverage_start.replace("Z", "+00:00"))) - import ipdb - ipdb.set_trace() url = self._get_odap_url(ncml_url, np.round(dt)) super(Mapper, self).__init__(url, gdal_dataset, gdal_metadata, *args, **kwargs) From 03e579a8fac6c53998aed0b18d12e967c895b3f1 Mon Sep 17 00:00:00 2001 From: Morten Wergeland Hansen Date: Tue, 14 May 2024 11:01:42 +0200 Subject: [PATCH 43/57] bug fix --- nansat/mappers/mapper_meps_ncml.py | 12 +++++++----- 1 file changed, 7 insertions(+), 5 deletions(-) diff --git a/nansat/mappers/mapper_meps_ncml.py b/nansat/mappers/mapper_meps_ncml.py index 44ad98bbb..c5fe86c74 100644 --- a/nansat/mappers/mapper_meps_ncml.py +++ b/nansat/mappers/mapper_meps_ncml.py @@ -1,4 +1,5 @@ import os +import pytz import netCDF4 import datetime @@ -15,13 +16,14 @@ def __init__(self, ncml_url, gdal_dataset, gdal_metadata, netcdf_dim=None, *args if not ncml_url.endswith(".ncml"): raise WrongMapperError - dt = 0 + 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"] - dt = time - np.datetime64( - datetime.datetime.fromisoformat(ds.time_coverage_start.replace("Z", "+00:00"))) - url = self._get_odap_url(ncml_url, np.round(dt)) + 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) From a559d6190993200ee1a1639eb90918d666475e14 Mon Sep 17 00:00:00 2001 From: Morten Wergeland Hansen Date: Tue, 14 May 2024 16:12:57 +0200 Subject: [PATCH 44/57] started --- nansat/mappers/sentinel1.py | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/nansat/mappers/sentinel1.py b/nansat/mappers/sentinel1.py index d763e2809..506b559b8 100644 --- a/nansat/mappers/sentinel1.py +++ b/nansat/mappers/sentinel1.py @@ -120,10 +120,15 @@ def get_gcp_shape(self): that the gcp_index dimension divided by gcp_x and gcp_y dimensions results in an integer. """ + # test modulo fra -/+ gcp_x/4 til remainder=0 + import ipdb + ipdb.set_trace() + if gcp_y*gcp_x != pixel.size: gcp_y = gcp_y - pixel.size % gcp_y gcp_x = gcp_x - pixel.size % gcp_x - assert gcp_y*gcp_x == pixel.size + if gcp_y*gcp_x != pixel.size: + raise ValueError("GCP dimension mismatch") return gcp_y, gcp_x From 2744857e2d9f9564c18c1a653730f0d43de43460 Mon Sep 17 00:00:00 2001 From: Morten Wergeland Hansen Date: Wed, 15 May 2024 10:39:04 +0200 Subject: [PATCH 45/57] Fix gcp shape by looping --- nansat/mappers/sentinel1.py | 23 +++++++++++++++++------ 1 file changed, 17 insertions(+), 6 deletions(-) diff --git a/nansat/mappers/sentinel1.py b/nansat/mappers/sentinel1.py index 506b559b8..4121e87ae 100644 --- a/nansat/mappers/sentinel1.py +++ b/nansat/mappers/sentinel1.py @@ -120,13 +120,24 @@ def get_gcp_shape(self): that the gcp_index dimension divided by gcp_x and gcp_y dimensions results in an integer. """ - # test modulo fra -/+ gcp_x/4 til remainder=0 - import ipdb - ipdb.set_trace() + def test(gcp_dim): + """ Test modulo from -/+ gcp_dim/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 = test(gcp_y) + gcp_x = test(gcp_x) - if gcp_y*gcp_x != pixel.size: - gcp_y = gcp_y - pixel.size % gcp_y - gcp_x = gcp_x - pixel.size % gcp_x if gcp_y*gcp_x != pixel.size: raise ValueError("GCP dimension mismatch") From d0d4f729ae7729dbb53ecbc1d98a4cb3dee9d62b Mon Sep 17 00:00:00 2001 From: Morten Wergeland Hansen Date: Thu, 16 May 2024 11:09:30 +0200 Subject: [PATCH 46/57] improve gcp reading --- nansat/mappers/sentinel1.py | 26 +++++++++++++++++++++----- 1 file changed, 21 insertions(+), 5 deletions(-) diff --git a/nansat/mappers/sentinel1.py b/nansat/mappers/sentinel1.py index 4121e87ae..2491b7585 100644 --- a/nansat/mappers/sentinel1.py +++ b/nansat/mappers/sentinel1.py @@ -120,8 +120,25 @@ def get_gcp_shape(self): that the gcp_index dimension divided by gcp_x and gcp_y dimensions results in an integer. """ - def test(gcp_dim): - """ Test modulo from -/+ gcp_dim/4 until remainder=0 + 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")): @@ -135,13 +152,12 @@ def test(gcp_dim): break return gcp_dim - gcp_y = test(gcp_y) - gcp_x = test(gcp_x) + gcp_y, gcp_x = test1(gcp_y, gcp_x) if gcp_y*gcp_x != pixel.size: raise ValueError("GCP dimension mismatch") - return gcp_y, gcp_x + return int(gcp_y), int(gcp_x) def add_incidence_angle_band(self): gcp_y, gcp_x = self.get_gcp_shape() From a7c92afb267520c4fe57c62cb2641d8eb2a8a2df Mon Sep 17 00:00:00 2001 From: Morten Wergeland Hansen Date: Thu, 16 May 2024 12:39:41 +0200 Subject: [PATCH 47/57] log gcp shape --- nansat/mappers/sentinel1.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/nansat/mappers/sentinel1.py b/nansat/mappers/sentinel1.py index 2491b7585..57e7fabab 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: @@ -154,6 +155,9 @@ def test2(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") From 09c607965b0a65a2ba670861addff238bbafb690 Mon Sep 17 00:00:00 2001 From: "Morten W. Hansen" Date: Fri, 7 Jun 2024 14:25:27 +0000 Subject: [PATCH 48/57] Init from lonlat --- nansat/mappers/mapper_netcdf_cf.py | 18 +++++++++++++++--- nansat/nsr.py | 2 +- 2 files changed, 16 insertions(+), 4 deletions(-) diff --git a/nansat/mappers/mapper_netcdf_cf.py b/nansat/mappers/mapper_netcdf_cf.py index e481209bb..2551ab7b7 100644 --- a/nansat/mappers/mapper_netcdf_cf.py +++ b/nansat/mappers/mapper_netcdf_cf.py @@ -497,9 +497,21 @@ def _create_empty_from_subdatasets(self, gdal_dataset, metadata): no_projection = False break if no_projection: - raise NansatMissingProjectionError - # Initialise VRT with subdataset containing projection - self._init_from_gdal_dataset(sub, metadata=metadata) + sub_datasets = gdal_dataset.GetSubDatasets() + filenames = [f[0] for f in sub_datasets] + for _, filename in enumerate(filenames): + if 'longitude' in filename: + xDatasetSource = filename + if 'latitude' in filename: + yDatasetSource = filename + lon_dataset = gdal.Open(xDatasetSource) + lon = lon_dataset.GetRasterBand(1).ReadAsArray() + lat_dataset = gdal.Open(yDatasetSource) + lat = lat_dataset.GetRasterBand(1).ReadAsArray() + self._init_from_lonlat(lon, lat, add_gcps=False) + else: + # Initialise VRT with subdataset containing projection + self._init_from_gdal_dataset(sub, metadata=metadata) def _create_empty_with_nansat_spatial_reference_WKT(self, gdal_dataset, gdal_metadata): """ In this case, gdal cannot find the projection of any diff --git a/nansat/nsr.py b/nansat/nsr.py index 46247bb29..5937ad744 100644 --- a/nansat/nsr.py +++ b/nansat/nsr.py @@ -60,7 +60,7 @@ def __init__(self, srs=0): osr.SpatialReference.__init__(self) # parse input parameters - if srs is 0: + if srs == 0: # generate default WGS84 SRS self.ImportFromEPSG(4326) elif isinstance(srs, str_types): From 7c5af871a5eae24d741d1253288392bfbe70fe4a Mon Sep 17 00:00:00 2001 From: "Morten W. Hansen" Date: Sun, 9 Jun 2024 11:49:41 +0000 Subject: [PATCH 49/57] Use lat/lon grids for initialization --- nansat/mappers/mapper_netcdf_cf.py | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/nansat/mappers/mapper_netcdf_cf.py b/nansat/mappers/mapper_netcdf_cf.py index 2551ab7b7..637fb9384 100644 --- a/nansat/mappers/mapper_netcdf_cf.py +++ b/nansat/mappers/mapper_netcdf_cf.py @@ -497,6 +497,7 @@ def _create_empty_from_subdatasets(self, gdal_dataset, metadata): no_projection = False break if no_projection: + #raise NansatMissingProjectionError sub_datasets = gdal_dataset.GetSubDatasets() filenames = [f[0] for f in sub_datasets] for _, filename in enumerate(filenames): @@ -504,6 +505,10 @@ def _create_empty_from_subdatasets(self, gdal_dataset, metadata): xDatasetSource = filename if 'latitude' in filename: yDatasetSource = filename + if not "xDatasetSource" in locals(): + raise NansatMissingProjectionError + if not "yDatasetSource" in locals(): + raise NansatMissingProjectionError lon_dataset = gdal.Open(xDatasetSource) lon = lon_dataset.GetRasterBand(1).ReadAsArray() lat_dataset = gdal.Open(yDatasetSource) From a6af99f550d402d501121a1d5234a659487fed74 Mon Sep 17 00:00:00 2001 From: "Morten W. Hansen" Date: Thu, 20 Jun 2024 14:19:03 +0200 Subject: [PATCH 50/57] bug fix? --- nansat/exporter.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/nansat/exporter.py b/nansat/exporter.py index 25c50a872..538849723 100644 --- a/nansat/exporter.py +++ b/nansat/exporter.py @@ -437,8 +437,8 @@ def _post_proc_thredds(self, fill_value = None if '_FillValue' in inp_var.ncattrs(): fill_value = inp_var._FillValue - if '_FillValue' in band_metadata[inp_var_name]: - fill_value = band_metadata['_FillValue'] + elif '_FillValue' in band_metadata[inp_var_name]: + fill_value = band_metadata[inp_var_name]['_FillValue'] dimensions = ('time', ) + inp_var.dimensions out_var = Exporter._copy_nc_var(inp_var, nc_out, inp_var_name, inp_var.dtype, dimensions, fill_value=fill_value, zlib=zlib) From ed63700a83a9b76e3d385e40a39ea44931040d46 Mon Sep 17 00:00:00 2001 From: Morten Wergeland Hansen Date: Tue, 25 Jun 2024 11:30:09 +0200 Subject: [PATCH 51/57] bug fix --- nansat/exporter.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/nansat/exporter.py b/nansat/exporter.py index 25c50a872..63138afe3 100644 --- a/nansat/exporter.py +++ b/nansat/exporter.py @@ -438,7 +438,7 @@ def _post_proc_thredds(self, if '_FillValue' in inp_var.ncattrs(): fill_value = inp_var._FillValue if '_FillValue' in band_metadata[inp_var_name]: - fill_value = band_metadata['_FillValue'] + fill_value = band_metadata[inp_var_name]['_FillValue'] dimensions = ('time', ) + inp_var.dimensions out_var = Exporter._copy_nc_var(inp_var, nc_out, inp_var_name, inp_var.dtype, dimensions, fill_value=fill_value, zlib=zlib) From 7e576c3ffa1e1f1c223ebe84d2af0d7206ec06c7 Mon Sep 17 00:00:00 2001 From: Morten Wergeland Hansen Date: Tue, 25 Jun 2024 15:04:24 +0200 Subject: [PATCH 52/57] Adjust gcp sizes --- nansat/mappers/sentinel1.py | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/nansat/mappers/sentinel1.py b/nansat/mappers/sentinel1.py index 57e7fabab..90c1d11f5 100644 --- a/nansat/mappers/sentinel1.py +++ b/nansat/mappers/sentinel1.py @@ -157,9 +157,13 @@ def test2(gcp_dim): logging.debug("GCPY size: %d" % gcp_y) logging.debug("GCPX size: %d" % gcp_x) + logging.debug("Pixel(s) size: %d" % pixel.size) if gcp_y*gcp_x != pixel.size: - raise ValueError("GCP dimension mismatch") + if gcp_y*gcp_x > pixel.size: + gcp_x, gcp_y = test1(gcp_x-1, gcp_y-1) + if gcp_y*gcp_x != pixel.size: + raise ValueError("GCP dimension mismatch") return int(gcp_y), int(gcp_x) From 3e98589c5a35f43dbe218e00868ddde86b5ccdef Mon Sep 17 00:00:00 2001 From: "Morten W. Hansen" Date: Fri, 30 Aug 2024 11:01:18 +0000 Subject: [PATCH 53/57] get rid of potential norwegian characters --- nansat/mappers/mapper_netcdf_cf.py | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/nansat/mappers/mapper_netcdf_cf.py b/nansat/mappers/mapper_netcdf_cf.py index 637fb9384..9b75c4f44 100644 --- a/nansat/mappers/mapper_netcdf_cf.py +++ b/nansat/mappers/mapper_netcdf_cf.py @@ -88,6 +88,11 @@ def __init__(self, filename, gdal_dataset, gdal_metadata, *args, **kwargs): # Set GCMD/DIF compatible metadata if available self._set_time_coverage_metadata(metadata) + # Add global metadata + metadata.pop("title_no", "") + metadata.pop("summary_no", "") + self.dataset.SetMetadata(metadata) + # Then add remaining GCMD/DIF compatible metadata in inheriting mappers def times(self): From 6b244c6d559a960d08c6eac3763f7c31bde661a7 Mon Sep 17 00:00:00 2001 From: Morten Wergeland Hansen Date: Fri, 25 Oct 2024 16:27:22 +0200 Subject: [PATCH 54/57] new mapper --- .../mapper_ncep_global_atm_model_best.py | 91 +++++++++++++++++++ nansat/mappers/mapper_opendap_arome.py | 3 + nansat/mappers/mapper_opendap_met_nordic.py | 3 + 3 files changed, 97 insertions(+) create mode 100644 nansat/mappers/mapper_ncep_global_atm_model_best.py diff --git a/nansat/mappers/mapper_ncep_global_atm_model_best.py b/nansat/mappers/mapper_ncep_global_atm_model_best.py new file mode 100644 index 000000000..56c7f6420 --- /dev/null +++ b/nansat/mappers/mapper_ncep_global_atm_model_best.py @@ -0,0 +1,91 @@ +import pytz +import netCDF4 +import datetime + +import numpy as np + +from dateutil.parser import parse + +from nansat.exceptions import WrongMapperError +from nansat.vrt import VRT + + +class Mapper(VRT): + + def __init__(self, filename, gdal_dataset, metadata, netcdf_dim=None, *args, **kwargs): + """Read NCEP_Global_Atmospheric_Model_best.ncd from + https://pae-paha.pacioos.hawaii.edu/thredds/dodsC/ncep_global/, + and get Nansat object with arrays as close as possible to the + given time. + + Time must be type datetime with timezone, e.g., + datetime.datetime(2024, 9, 10, 15, 0, 0, tzinfo=pytz.utc). + """ + + fn = ("https://pae-paha.pacioos.hawaii.edu/thredds/dodsC/" + "ncep_global/NCEP_Global_Atmospheric_Model_best.ncd") + if filename != fn: + raise WrongMapperError + + if netcdf_dim is None: + raise WrongMapperError + + ds = netCDF4.Dataset(filename) + metadata = {} + for key in ds.ncattrs(): + metadata[key] = str(ds.getncattr(key)) + + lon = ds["longitude"][:] + lat = ds["latitude"][:] + times = ds["time"][:] + + longitude, latitude = np.meshgrid(lon, lat) + + super(Mapper, self)._init_from_lonlat(longitude, latitude) + self.dataset.SetMetadata(metadata) + + t0 = parse(ds["time"].units.strip("hours since ")) + t1 = ds["time"][:] + dt = [datetime.timedelta(hours=t) for t in t1.data[:]] + times = [t0 + tt for tt in dt] + diff = np.array(times) - \ + netcdf_dim["time"].astype(datetime.datetime).replace(tzinfo=pytz.utc) + time_index = np.abs(diff).argmin() + + for var, val in ds.variables.items(): + if "standard_name" in ds[var].ncattrs(): + if val.standard_name == "longitude": + lon = var + if val.standard_name == "latitude": + lat = var + + for attr in ds.ncattrs(): + content = ds.getncattr(attr) + if type(content) is str: + metadata[attr] = content + + self.band_vrts = {} + metaDict = [] + + for key, val in ds.variables.items(): + data = ds[key] + if data.shape != (t1.shape[0], longitude.shape[0], longitude.shape[1]): + continue + self.band_vrts[key] = VRT.from_array( + data[time_index, :, :].filled(fill_value=np.nan)) + key_metadict = {} + for attr in ds[key].ncattrs(): + key_metadict[attr] = ds[key].getncattr(attr) + key_metadict["name"] = key + key_metadict["time"] = times[time_index].isoformat() + metaDict.append({ + 'src': { + 'SourceFilename': self.band_vrts[key].filename, + 'SourceBand': 1 + }, + 'dst': key_metadict, + }) + + self.create_bands(metaDict) + + #self.fix_global_metadata diff --git a/nansat/mappers/mapper_opendap_arome.py b/nansat/mappers/mapper_opendap_arome.py index 6f45fdeba..1a4583bd1 100644 --- a/nansat/mappers/mapper_opendap_arome.py +++ b/nansat/mappers/mapper_opendap_arome.py @@ -39,6 +39,9 @@ def __init__(self, filename, gdal_dataset, gdal_metadata, netcdf_dim=None, self.input_filename = filename + if "thredds.met.no" not in filename: + raise WrongMapperError + if gdal_dataset is None: gdal_dataset = gdal.Open(filename) diff --git a/nansat/mappers/mapper_opendap_met_nordic.py b/nansat/mappers/mapper_opendap_met_nordic.py index d17b2b967..f21682133 100644 --- a/nansat/mappers/mapper_opendap_met_nordic.py +++ b/nansat/mappers/mapper_opendap_met_nordic.py @@ -31,6 +31,9 @@ def __init__(self, filename, gdal_dataset, gdal_metadata, netcdf_dim=None, if netcdf_dim is None: netcdf_dim = {} + if "thredds.met.no" not in filename: + raise WrongMapperError + self.input_filename = filename if gdal_dataset is None: From b1d7edea65abf55398e918efbf64ae96b0a90c74 Mon Sep 17 00:00:00 2001 From: "Morten W. Hansen" Date: Tue, 29 Oct 2024 15:40:45 +0100 Subject: [PATCH 55/57] correct docstring --- nansat/mappers/mapper_ncep_global_atm_model_best.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/nansat/mappers/mapper_ncep_global_atm_model_best.py b/nansat/mappers/mapper_ncep_global_atm_model_best.py index 56c7f6420..9a3675ddc 100644 --- a/nansat/mappers/mapper_ncep_global_atm_model_best.py +++ b/nansat/mappers/mapper_ncep_global_atm_model_best.py @@ -18,8 +18,7 @@ def __init__(self, filename, gdal_dataset, metadata, netcdf_dim=None, *args, **k and get Nansat object with arrays as close as possible to the given time. - Time must be type datetime with timezone, e.g., - datetime.datetime(2024, 9, 10, 15, 0, 0, tzinfo=pytz.utc). + Time must be type datetime64. Timezone is assumed to be UTC. """ fn = ("https://pae-paha.pacioos.hawaii.edu/thredds/dodsC/" From 191914f33c38967745960dc4812dd9184b179828 Mon Sep 17 00:00:00 2001 From: "Morten W. Hansen" Date: Tue, 29 Oct 2024 19:46:32 +0100 Subject: [PATCH 56/57] fix docstring --- nansat/nansat.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/nansat/nansat.py b/nansat/nansat.py index 7fabc6dd3..afe9869ca 100644 --- a/nansat/nansat.py +++ b/nansat/nansat.py @@ -1416,7 +1416,7 @@ def crop_lonlat(self, lonlim, latlim): Examples -------- - >>> extent = n.crop(lonlim=[-10,10], latlim=[-20,20]) # crop for given lon/lat limits + >>> extent = n.crop([-10,10], [-20,20]) # crop for given lon/lat limits """ # lon/lat lists for four corners From 2d634ceb9fb046703099a659a30910adf4daccd5 Mon Sep 17 00:00:00 2001 From: "Morten W. Hansen" Date: Wed, 30 Oct 2024 14:17:25 +0100 Subject: [PATCH 57/57] reorganise arrays --- .../mapper_ncep_global_atm_model_best.py | 30 ++++++++++--------- 1 file changed, 16 insertions(+), 14 deletions(-) diff --git a/nansat/mappers/mapper_ncep_global_atm_model_best.py b/nansat/mappers/mapper_ncep_global_atm_model_best.py index 9a3675ddc..e94024e93 100644 --- a/nansat/mappers/mapper_ncep_global_atm_model_best.py +++ b/nansat/mappers/mapper_ncep_global_atm_model_best.py @@ -34,13 +34,20 @@ def __init__(self, filename, gdal_dataset, metadata, netcdf_dim=None, *args, **k for key in ds.ncattrs(): metadata[key] = str(ds.getncattr(key)) - lon = ds["longitude"][:] + lon = (ds["longitude"][:] + 180.) % 360 - 180. + #lon = ds["longitude"][:] lat = ds["latitude"][:] times = ds["time"][:] longitude, latitude = np.meshgrid(lon, lat) - - super(Mapper, self)._init_from_lonlat(longitude, latitude) + # Shift array to yield -180 to 180 deg longitude + ind = np.where(lon<0)[0] + xx = np.empty(longitude.shape) + xx[:, 0:ind.shape[0]] = longitude[:, ind.min():] + xx[:, ind.shape[0]:] = longitude[:, 0:ind.min()] + longitude = xx + + super(Mapper, self)._init_from_lonlat(longitude, latitude, add_gcps=True) self.dataset.SetMetadata(metadata) t0 = parse(ds["time"].units.strip("hours since ")) @@ -51,13 +58,6 @@ def __init__(self, filename, gdal_dataset, metadata, netcdf_dim=None, *args, **k netcdf_dim["time"].astype(datetime.datetime).replace(tzinfo=pytz.utc) time_index = np.abs(diff).argmin() - for var, val in ds.variables.items(): - if "standard_name" in ds[var].ncattrs(): - if val.standard_name == "longitude": - lon = var - if val.standard_name == "latitude": - lat = var - for attr in ds.ncattrs(): content = ds.getncattr(attr) if type(content) is str: @@ -67,11 +67,13 @@ def __init__(self, filename, gdal_dataset, metadata, netcdf_dim=None, *args, **k metaDict = [] for key, val in ds.variables.items(): - data = ds[key] - if data.shape != (t1.shape[0], longitude.shape[0], longitude.shape[1]): + if ds[key].shape != (t1.shape[0], longitude.shape[0], longitude.shape[1]): continue - self.band_vrts[key] = VRT.from_array( - data[time_index, :, :].filled(fill_value=np.nan)) + data = ds[key][time_index, :, :].filled(fill_value=np.nan) + xx = np.empty(longitude.shape) + xx[:, 0:ind.shape[0]] = data[:, ind.min():] + xx[:, ind.shape[0]:] = data[:, 0:ind.min()] + self.band_vrts[key] = VRT.from_array(xx) key_metadict = {} for attr in ds[key].ncattrs(): key_metadict[attr] = ds[key].getncattr(attr)