diff --git a/CHANGELOG.rst b/CHANGELOG.rst index f94241b1be8f3d8bdc786a18fc2e9bdf8da154c3..4251c85b036bcb2d1e2d959a1a52521ac7d8190e 100644 --- a/CHANGELOG.rst +++ b/CHANGELOG.rst @@ -19,12 +19,16 @@ CHANGELOG * Improved load_nes.py removing redundant code * Direct access to variable data. (`#77 `_) * New functionalities for vertical extrapolation. (`#74 `_) + * Removed cfunits and psutil dependencies. + * Updated the requirements and environment.yml for the Conda environment created in MN5 (`#78 `_) * Bugfix: * Vertical interpolation for descendant level values (Pressure) (`#71 `_) * Removed lat-lon dimension on the NetCDF projections that not need them (`#72 `_) * Fixed the bug when creating the spatial bounds after selecting a region (`#68 `_) - + * Fixed the bug related to Shapely deprecated function TopologicalError(`#76 `_) + * Fixed the bug related to NumPy deprecated np.object(`#76 `_) + * Removed DeprecationWarnings from shapely and pyproj libraries, needed for the porting to MN5 (`#78 `_) 1.1.3 ============ diff --git a/environment.yml b/environment.yml index 810b716ad0272f873492d22d3f8b0cceb4b3f5b6..7283f91772879dfd9366460d15db933c43158dcf 100755 --- a/environment.yml +++ b/environment.yml @@ -1,14 +1,18 @@ --- -name: nes +name: NES_v1.1.3 channels: - conda-forge - - anaconda + - bioconda + - defaults dependencies: - - python = 3 - # Testing - - pytest - - pytest-cov - - pycodestyle \ No newline at end of file + - python >= 3.8 + - mpich + - mpi4py ~= 3.1.4 + - geopandas + - netcdf4=*=mpi* + - eccodes + - filelock + - configargparse \ No newline at end of file diff --git a/nes/load_nes.py b/nes/load_nes.py index 839b473660e4d88c62498788a64748222c0e12f6..b30086940f5d1af7f68fc66b0d9c790e10f4ad83 100644 --- a/nes/load_nes.py +++ b/nes/load_nes.py @@ -289,11 +289,11 @@ def concatenate_netcdfs(nessy_list, comm=None, info=False, parallel_method='Y', data = var_info[nessy_first.read_axis_limits['y_min']:nessy_first.read_axis_limits['y_max'], nessy_first.read_axis_limits['x_min']:nessy_first.read_axis_limits['x_max'], :] - data_aux = np.empty(shape=(data.shape[0], data.shape[1]), dtype=np.object) + data_aux = np.empty(shape=(data.shape[0], data.shape[1]), dtype=object) for lat_n in range(data.shape[0]): for lon_n in range(data.shape[1]): data_aux[lat_n, lon_n] = ''.join( - data[lat_n, lon_n].tostring().decode('ascii').replace('\x00', '')) + data[lat_n, lon_n].tobytes().decode('ascii').replace('\x00', '')) data = data_aux.reshape((1, 1, data_aux.shape[-2], data_aux.shape[-1])) else: data = var_info[nessy_first.read_axis_limits['t_min']:nessy_first.read_axis_limits['t_max'], diff --git a/nes/methods/cell_measures.py b/nes/methods/cell_measures.py index 3e561c554f604eb9acf9e5b8fdaa9ffdf708dc2e..5288a02ae50a533306f5bec7ca6b20c64b667ea9 100644 --- a/nes/methods/cell_measures.py +++ b/nes/methods/cell_measures.py @@ -81,10 +81,10 @@ def calculate_geometry_area(geometry_list, earth_radius_minor_axis=6356752.3142, for geom_ind in range(0, len(geometry_list)): # Calculate the area of each geometry in multipolygon and collection objects - if geometry_list[geom_ind].type in ['MultiPolygon', 'GeometryCollection']: + if geometry_list[geom_ind].geom_type in ['MultiPolygon', 'GeometryCollection']: multi_geom_area = 0 for multi_geom_ind in range(0, len(geometry_list[geom_ind].geoms)): - if geometry_list[geom_ind].geoms[multi_geom_ind].type == 'Point': + if geometry_list[geom_ind].geoms[multi_geom_ind].geom_type == 'Point': continue geometry_corner_lon, geometry_corner_lat = geometry_list[geom_ind].geoms[multi_geom_ind].exterior.coords.xy geometry_corner_lon = np.array(geometry_corner_lon) diff --git a/nes/methods/horizontal_interpolation.py b/nes/methods/horizontal_interpolation.py index 3887bd905b1dc47b8fac119a31fb9a331484b56e..006f5ffa780ee4330ed58c47199c73a72fecfe6e 100644 --- a/nes/methods/horizontal_interpolation.py +++ b/nes/methods/horizontal_interpolation.py @@ -15,7 +15,6 @@ from warnings import warn import copy import pyproj import gc -import psutil # CONSTANTS NEAREST_OPTS = ['NearestNeighbour', 'NearestNeighbours', 'nn', 'NN'] @@ -597,7 +596,7 @@ def create_area_conservative_weight_matrix(self, dst_nes, wm_path=None, flux=Fal src_grid['FID_src'] = src_grid.index src_grid = src_grid.reset_index() dst_grid = dst_grid.reset_index() - fid_src, fid_dst = dst_grid.sindex.query_bulk(src_grid.geometry, predicate='intersects') + fid_src, fid_dst = dst_grid.sindex.query(src_grid.geometry, predicate='intersects') # Calculate intersected areas and fractions intersection_df = pd.DataFrame(columns=["FID_src", "FID_dst"]) @@ -755,6 +754,8 @@ def lon_lat_to_cartesian_ecef(lon, lat): lla = pyproj.Proj(proj='latlong', ellps='WGS84', datum='WGS84') ecef = pyproj.Proj(proj='geocent', ellps='WGS84', datum='WGS84') - x, y, z = pyproj.transform(lla, ecef, lon, lat, np.zeros(lon.shape), radians=False) - + # x, y, z = pyproj.transform(lla, ecef, lon, lat, np.zeros(lon.shape), radians=False) + # Deprecated: https://pyproj4.github.io/pyproj/stable/gotchas.html#upgrading-to-pyproj-2-from-pyproj-1 + transformer = pyproj.Transformer.from_proj(lla, ecef) + x, y, z = transformer.transform(lon, lat, np.zeros(lon.shape), radians=False) return np.column_stack([x, y, z]) diff --git a/nes/methods/spatial_join.py b/nes/methods/spatial_join.py index 52c51ed7fece30164a173903f56ff7a43a36024f..d647ebb66a6753fd1ac9b1b7c6400dd4c375c343 100644 --- a/nes/methods/spatial_join.py +++ b/nes/methods/spatial_join.py @@ -7,7 +7,7 @@ from geopandas import GeoDataFrame import nes import numpy as np import pandas as pd -from shapely.geos import TopologicalError +from shapely.errors import TopologicalError def spatial_join(self, ext_shp, method=None, var_list=None, info=False, apply_bbox=True): @@ -247,8 +247,8 @@ def spatial_join_intersection(self, ext_shp, info=False): grid_shp = grid_shp.reset_index() # Get intersected areas - # inp, res = ext_shp.sindex.query_bulk(grid_shp.geometry, predicate='intersects') - inp, res = grid_shp.sindex.query_bulk(ext_shp.geometry, predicate='intersects') + # inp, res = ext_shp.sindex.query(grid_shp.geometry, predicate='intersects') + inp, res = grid_shp.sindex.query(ext_shp.geometry, predicate='intersects') if info: print('\t\tRank {0:03d}: {1} intersected areas found'.format(self.rank, len(inp))) diff --git a/nes/methods/vertical_interpolation.py b/nes/methods/vertical_interpolation.py index 7260fdbe191484129649b0feff7e3efcc4d95728..d1868e413043a90d566db81ea2eb3dba68565621 100644 --- a/nes/methods/vertical_interpolation.py +++ b/nes/methods/vertical_interpolation.py @@ -328,6 +328,6 @@ def interpolate_vertical(self, new_levels, new_src_vertical=None, kind='linear', # Remove original file information self.__ini_path = None self.dataset = None - self.netcdf = None + self.dataset = None return self diff --git a/nes/nc_projections/default_nes.py b/nes/nc_projections/default_nes.py index e3fb5ecae8d96932b311c2d251fdc17447deca41..0fb60ac77597190944160021bd84dcba81f7640f 100644 --- a/nes/nc_projections/default_nes.py +++ b/nes/nc_projections/default_nes.py @@ -8,8 +8,7 @@ import pandas as pd from datetime import timedelta from netCDF4 import Dataset, num2date, date2num, stringtochar from mpi4py import MPI -from cfunits import Units -from shapely.geos import TopologicalError +from shapely.errors import TopologicalError import geopandas as gpd from shapely.geometry import Polygon, Point from copy import deepcopy, copy @@ -2004,7 +2003,7 @@ class Nes(object): variables[var_name]['data'] = None variables[var_name]['dimensions'] = var_info.dimensions variables[var_name]['dtype'] = var_info.dtype - if variables[var_name]['dtype'] in [str, np.object]: + if variables[var_name]['dtype'] in [str, object]: if self.strlen is None: self.set_strlen() variables[var_name]['dtype'] = str @@ -2013,7 +2012,7 @@ class Nes(object): for attrname in var_info.ncattrs(): if attrname not in ['missing_value', '_FillValue']: value = getattr(var_info, attrname) - if value in ['unitless', '-']: + if str(value) in ['unitless', '-']: value = '' variables[var_name][attrname] = value else: @@ -2052,11 +2051,11 @@ class Nes(object): data = nc_var[self.read_axis_limits['y_min']:self.read_axis_limits['y_max'], self.read_axis_limits['x_min']:self.read_axis_limits['x_max'], :] - data_aux = np.empty(shape=(data.shape[0], data.shape[1]), dtype=np.object) + data_aux = np.empty(shape=(data.shape[0], data.shape[1]), dtype=object) for lat_n in range(data.shape[0]): for lon_n in range(data.shape[1]): data_aux[lat_n, lon_n] = ''.join( - data[lat_n, lon_n].tostring().decode('ascii').replace('\x00', '')) + data[lat_n, lon_n].tobytes().decode('ascii').replace('\x00', '')) data = data_aux.reshape((1, 1, data_aux.shape[-2], data_aux.shape[-1])) else: data = nc_var[self.read_axis_limits['t_min']:self.read_axis_limits['t_max'], @@ -2075,13 +2074,13 @@ class Nes(object): self.read_axis_limits['y_min']:self.read_axis_limits['y_max'], self.read_axis_limits['x_min']:self.read_axis_limits['x_max'], :] - data_aux = np.empty(shape=(data.shape[0], data.shape[1], data.shape[2], data.shape[3]), dtype=np.object) + data_aux = np.empty(shape=(data.shape[0], data.shape[1], data.shape[2], data.shape[3]), dtype=object) for time_n in range(data.shape[0]): for lev_n in range(data.shape[1]): for lat_n in range(data.shape[2]): for lon_n in range(data.shape[3]): data_aux[time_n, lev_n, lat_n, lon_n] = ''.join( - data[time_n, lev_n, lat_n, lon_n].tostring().decode('ascii').replace('\x00', '')) + data[time_n, lev_n, lat_n, lon_n].tobytes().decode('ascii').replace('\x00', '')) data = data_aux else: # data = nc_var[self.read_axis_limits['t_min']:self.read_axis_limits['t_max'], @@ -2452,7 +2451,7 @@ class Nes(object): lev = netcdf.createVariable('lev', np.float32, ('lev',), zlib=self.zip_lvl > 0, complevel=self.zip_lvl) if 'units' in self._lev.keys(): - lev.units = Units(self._lev['units'], formatted=True).units + lev.units = self._lev['units'] else: lev.units = '' if 'positive' in self._lev.keys(): @@ -2549,7 +2548,7 @@ class Nes(object): lev = netcdf.createVariable('lev', self._lev['data'].dtype, ('lev',), zlib=self.zip_lvl > 0, complevel=self.zip_lvl) if 'units' in self._lev.keys(): - lev.units = Units(self._lev['units'], formatted=True).units + lev.units = self._lev['units'] else: lev.units = '' if 'positive' in self._lev.keys(): @@ -2687,8 +2686,8 @@ class Nes(object): raise TypeError("It was not possible to cast the data to the input dtype.") else: var_dtype = var_dict['data'].dtype - if var_dtype is np.object: - raise TypeError("Data dtype is np.object. Define dtype explicitly as dictionary key 'dtype'") + if var_dtype is object: + raise TypeError("Data dtype is object. Define dtype explicitly as dictionary key 'dtype'") if var_dict['data'] is not None: @@ -3534,8 +3533,10 @@ class Nes(object): lla = pyproj.Proj(proj='latlong', ellps='WGS84', datum='WGS84') ecef = pyproj.Proj(proj='geocent', ellps='WGS84', datum='WGS84') - x, y, z = pyproj.transform(lla, ecef, lon, lat, np.zeros(lon.shape), radians=False) - + # x, y, z = pyproj.transform(lla, ecef, lon, lat, np.zeros(lon.shape), radians=False) + # Deprecated: https://pyproj4.github.io/pyproj/stable/gotchas.html#upgrading-to-pyproj-2-from-pyproj-1 + transformer = pyproj.Transformer.from_proj(lla, ecef) + x, y, z = transformer.transform(lon, lat, np.zeros(lon.shape), radians=False) return np.column_stack([x, y, z]) def add_4d_vertical_info(self, info_to_add): diff --git a/nes/nc_projections/lcc_nes.py b/nes/nc_projections/lcc_nes.py index c225d5b769956068cee4c84d8e33bdc9cf24cd4b..79417ac87f266e6f3cfe86267fdad144a2788cb0 100644 --- a/nes/nc_projections/lcc_nes.py +++ b/nes/nc_projections/lcc_nes.py @@ -4,7 +4,6 @@ import warnings import sys import numpy as np import pandas as pd -from cfunits import Units from pyproj import Proj from copy import deepcopy import geopandas as gpd @@ -279,7 +278,7 @@ class LCCNes(Nes): y = netcdf.createVariable('y', self._y['data'].dtype, ('y',)) y.long_name = 'y coordinate of projection' if 'units' in self._y.keys(): - y.units = Units(self._y['units'], formatted=True).units + y.units = self._y['units'] else: y.units = 'm' y.standard_name = 'projection_y_coordinate' @@ -291,7 +290,7 @@ class LCCNes(Nes): x = netcdf.createVariable('x', self._x['data'].dtype, ('x',)) x.long_name = 'x coordinate of projection' if 'units' in self._x.keys(): - x.units = Units(self._x['units'], formatted=True).units + x.units = self._x['units'] else: x.units = 'm' x.standard_name = 'projection_x_coordinate' diff --git a/nes/nc_projections/mercator_nes.py b/nes/nc_projections/mercator_nes.py index 3cd46182389d258bd853518e9e67b08f707a0f95..263ecc74693f348c1365c0712a8bf048fcbd25f3 100644 --- a/nes/nc_projections/mercator_nes.py +++ b/nes/nc_projections/mercator_nes.py @@ -4,7 +4,6 @@ import warnings import sys import numpy as np import pandas as pd -from cfunits import Units from pyproj import Proj from copy import deepcopy import geopandas as gpd @@ -264,7 +263,7 @@ class MercatorNes(Nes): y = netcdf.createVariable('y', self._y['data'].dtype, ('y',)) y.long_name = 'y coordinate of projection' if 'units' in self._y.keys(): - y.units = Units(self._y['units'], formatted=True).units + y.units = self._y['units'] else: y.units = 'm' y.standard_name = 'projection_y_coordinate' @@ -276,7 +275,7 @@ class MercatorNes(Nes): x = netcdf.createVariable('x', self._x['data'].dtype, ('x',)) x.long_name = 'x coordinate of projection' if 'units' in self._x.keys(): - x.units = Units(self._x['units'], formatted=True).units + x.units = self._x['units'] else: x.units = 'm' x.standard_name = 'projection_x_coordinate' diff --git a/nes/nc_projections/points_nes.py b/nes/nc_projections/points_nes.py index bca62c070b24f1d72614cb7c416b181161a57d97..0a762f86ba0dccd6581c7257c44740e69a65195a 100644 --- a/nes/nc_projections/points_nes.py +++ b/nes/nc_projections/points_nes.py @@ -310,7 +310,7 @@ class PointsNes(Nes): Portion of the variable data corresponding to the rank. """ - nc_var = self.netcdf.variables[var_name] + nc_var = self.dataset.variables[var_name] var_dims = nc_var.dimensions # Read data in 1 or 2 dimensions @@ -319,8 +319,8 @@ class PointsNes(Nes): elif len(var_dims) == 2: if 'strlen' in var_dims: data = nc_var[self.read_axis_limits['x_min']:self.read_axis_limits['x_max'], :] - data = np.array([''.join(i.tostring().decode('ascii').replace('\x00', '')) for i in data], - dtype=np.object) + data = np.array([''.join(i.tobytes().decode('ascii').replace('\x00', '')) for i in data], + dtype=object) else: data = nc_var[self.read_axis_limits['t_min']:self.read_axis_limits['t_max'], self.read_axis_limits['x_min']:self.read_axis_limits['x_max']] @@ -362,8 +362,8 @@ class PointsNes(Nes): raise e("It was not possible to cast the data to the input dtype.") else: var_dtype = var_dict['data'].dtype - if var_dtype is np.object: - raise TypeError("Data dtype is np.object. Define dtype explicitly as dictionary key 'dtype'") + if var_dtype is object: + raise TypeError("Data dtype is object. Define dtype explicitly as dictionary key 'dtype'") # Get dimensions when reading datasets if 'dimensions' in var_dict.keys(): diff --git a/nes/nc_projections/points_nes_ghost.py b/nes/nc_projections/points_nes_ghost.py index 2acdd81d99db809a0271310f51fcbe0005e43958..f91d565be7ee565fb3cdf9a8a590a71fc0760189 100644 --- a/nes/nc_projections/points_nes_ghost.py +++ b/nes/nc_projections/points_nes_ghost.py @@ -279,7 +279,7 @@ class PointsNesGHOST(PointsNes): Portion of the variable data corresponding to the rank. """ - nc_var = self.netcdf.variables[var_name] + nc_var = self.dataset.variables[var_name] var_dims = nc_var.dimensions # Read data in 1 or 2 dimensions @@ -330,8 +330,8 @@ class PointsNesGHOST(PointsNes): raise e("It was not possible to cast the data to the input dtype.") else: var_dtype = var_dict['data'].dtype - if var_dtype is np.object: - raise TypeError("Data dtype is np.object. Define dtype explicitly as dictionary key 'dtype'") + if var_dtype is object: + raise TypeError("Data dtype is object. Define dtype explicitly as dictionary key 'dtype'") # Get dimensions when reading datasets if 'dimensions' in var_dict.keys(): diff --git a/nes/nc_projections/points_nes_providentia.py b/nes/nc_projections/points_nes_providentia.py index d2643f5f147f15552cd420299b1de9483624d117..e35ff3d2789dc3a873bcba595ccd64884fb312c9 100644 --- a/nes/nc_projections/points_nes_providentia.py +++ b/nes/nc_projections/points_nes_providentia.py @@ -314,8 +314,7 @@ class PointsNesProvidentia(PointsNes): data: np.array Portion of the variable data corresponding to the rank. """ - - nc_var = self.netcdf.variables[var_name] + nc_var = self.dataset.variables[var_name] var_dims = nc_var.dimensions # Read data in 1, 2 or 3 dimensions @@ -367,8 +366,8 @@ class PointsNesProvidentia(PointsNes): raise e("It was not possible to cast the data to the input dtype.") else: var_dtype = var_dict['data'].dtype - if var_dtype is np.object: - raise TypeError("Data dtype is np.object. Define dtype explicitly as dictionary key 'dtype'") + if var_dtype is yobject: + raise TypeError("Data dtype is object. Define dtype explicitly as dictionary key 'dtype'") # Get dimensions when reading datasets if 'dimensions' in var_dict.keys(): diff --git a/nes/nc_projections/rotated_nes.py b/nes/nc_projections/rotated_nes.py index b6f86c1501d67910bc647a27a51a8c34b0d8b425..c29b1fd9ad869b5e4f350499d04649ea10f4fc91 100644 --- a/nes/nc_projections/rotated_nes.py +++ b/nes/nc_projections/rotated_nes.py @@ -5,7 +5,6 @@ import sys import numpy as np import pandas as pd import math -from cfunits import Units from pyproj import Proj from copy import deepcopy import geopandas as gpd @@ -260,7 +259,7 @@ class RotatedNes(Nes): rlat = netcdf.createVariable('rlat', self._rlat['data'].dtype, ('rlat',)) rlat.long_name = "latitude in rotated pole grid" if 'units' in self._rlat.keys(): - rlat.units = Units(self._rlat['units'], formatted=True).units + rlat.units = self._rlat['units'] else: rlat.units = 'degrees' rlat.standard_name = "grid_latitude" @@ -272,7 +271,7 @@ class RotatedNes(Nes): rlon = netcdf.createVariable('rlon', self._rlon['data'].dtype, ('rlon',)) rlon.long_name = "longitude in rotated pole grid" if 'units' in self._rlon.keys(): - rlon.units = Units(self._rlon['units'], formatted=True).units + rlon.units = self._rlon['units'] else: rlon.units = 'degrees' rlon.standard_name = "grid_longitude" diff --git a/nes/nes_formats/cmaq_format.py b/nes/nes_formats/cmaq_format.py index f715b765562da6a544fc1075138f01a1e8964785..e8fd3990ee7d31dddb81a2c97025dbd76f362d19 100644 --- a/nes/nes_formats/cmaq_format.py +++ b/nes/nes_formats/cmaq_format.py @@ -60,7 +60,7 @@ def to_netcdf_cmaq(self, path, chunking=False, keep_open=False): # Close NetCDF if keep_open: - self.netcdf = netcdf + self.dataset = netcdf else: netcdf.close() diff --git a/nes/nes_formats/monarch_format.py b/nes/nes_formats/monarch_format.py index a8f5e24753395200911efa47f3b3a9e2c139227f..34a22d10af47ac09313623379d42fb1bbc5778a1 100644 --- a/nes/nes_formats/monarch_format.py +++ b/nes/nes_formats/monarch_format.py @@ -73,7 +73,7 @@ def to_netcdf_monarch(self, path, chunking=False, keep_open=False): netcdf.setncattr('Conventions', 'CF-1.7') if keep_open: - self.netcdf = netcdf + self.dataset = netcdf else: netcdf.close() diff --git a/nes/nes_formats/wrf_chem_format.py b/nes/nes_formats/wrf_chem_format.py index f37a74d4618be1c83889b5dc1c3e0bd7fad19a77..d2a71cabe0d4f919de8cd5dde76c5f9af82e9418 100644 --- a/nes/nes_formats/wrf_chem_format.py +++ b/nes/nes_formats/wrf_chem_format.py @@ -68,7 +68,7 @@ def to_netcdf_wrf_chem(self, path, chunking=False, keep_open=False): # Close NetCDF if keep_open: - self.netcdf = netcdf + self.dataset = netcdf else: netcdf.close() diff --git a/requirements.txt b/requirements.txt index 7bdad0d8fa2b2e14b2c548c652ba5ef1d817d62f..c31235437f78ae07517150e6056c96006101cf56 100644 --- a/requirements.txt +++ b/requirements.txt @@ -1,5 +1,5 @@ pycodestyle>=2.10.0 -geopandas>=0.10.2 +geopandas>=1.0.0 rtree>=0.9.0 pandas>=1.3.5 netcdf4>=1.6.2 @@ -11,9 +11,7 @@ Shapely>=2.0.0 scipy>=1.7.3 filelock>=3.9.0 eccodes-python~=0.9.5 -cfunits>=3.3.5 mpi4py>=3.1.4 # rasterio>=1.1.3 sphinx>=7.2.6 -sphinx-rtd-theme==2.0.0 -psutil>=5.9.6 \ No newline at end of file +sphinx-rtd-theme==2.0.0 \ No newline at end of file