CHANGELOG.rst

@@ -19,12 +19,16 @@ CHANGELOG
   * Improved load_nes.py removing redundant code
   * Direct access to variable data. (`#77 <https://earth.bsc.es/gitlab/es/NES/-/issues/77>`_)
   * New functionalities for vertical extrapolation. (`#74 <https://earth.bsc.es/gitlab/es/NES/-/issues/74>`_)
+  * Removed cfunits and psutil dependencies.
+  * Updated the requirements and environment.yml for the Conda environment created in MN5 (`#78 <https://earth.bsc.es/gitlab/es/NES/-/issues/78>`_)
   * Bugfix:
 
     * Vertical interpolation for descendant level values (Pressure) (`#71 <https://earth.bsc.es/gitlab/es/NES/-/issues/71>`_)
     * Removed lat-lon dimension on the NetCDF projections that not need them (`#72 <https://earth.bsc.es/gitlab/es/NES/-/issues/72>`_)
     * Fixed the bug when creating the spatial bounds after selecting a region (`#68 <https://earth.bsc.es/gitlab/es/NES/-/issues/68>`_)
-
+    * Fixed the bug related to Shapely deprecated function TopologicalError(`#76 <https://earth.bsc.es/gitlab/es/NES/-/issues/76>`_)
+    * Fixed the bug related to NumPy deprecated np.object(`#76 <https://earth.bsc.es/gitlab/es/NES/-/issues/76>`_)
+    * Removed DeprecationWarnings from shapely and pyproj libraries, needed for the porting to MN5 (`#78 <https://earth.bsc.es/gitlab/es/NES/-/issues/78>`_)
 1.1.3
 ============
 

environment.yml

 ---
 
-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

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'],

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)

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])

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)))

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

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):

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'

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'

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():

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():

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():

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"

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()
 

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()
 

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()
 

requirements.txt

 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