From 4f674897ebdf1961d7c2a4d7ce3f700f1881a686 Mon Sep 17 00:00:00 2001 From: Alba Vilanova Cortezon Date: Wed, 31 Aug 2022 12:44:17 +0200 Subject: [PATCH] #3 Implement to_shapefile, update tutorials and correct bounds --- nes/nc_projections/default_nes.py | 51 +++++ nes/nc_projections/lcc_nes.py | 168 +++++++++----- nes/nc_projections/mercator_nes.py | 112 +++++++--- nes/nc_projections/points_nes.py | 48 ++-- nes/nc_projections/rotated_nes.py | 71 ++++-- tests/scalability_test_nord3v2.bash | 3 +- tests/test_bash_mn4.cmd | 3 +- tests/test_bash_nord3v2.cmd | 2 + .../1.5.Read_Write_Mercator.ipynb | 113 ++-------- tutorials/2.Creation/2.2.Create_Rotated.ipynb | 16 +- tutorials/2.Creation/2.5.Create-LCC.ipynb | 16 +- .../2.Creation/2.6.Create_Mercator.ipynb | 62 +++--- .../4.3.Providentia_Interpolation.ipynb | 42 ++-- .../5.Others/5.2.Add_Coordinates_Bounds.ipynb | 178 +++++++-------- tutorials/5.Others/5.3.To_Shapefile.ipynb | 210 ++++++++++++++++++ tutorials/Jupyter_bash_nord3v2.cmd | 2 + 16 files changed, 713 insertions(+), 384 deletions(-) create mode 100644 tutorials/5.Others/5.3.To_Shapefile.ipynb diff --git a/nes/nc_projections/default_nes.py b/nes/nc_projections/default_nes.py index 81ab079..aca209c 100644 --- a/nes/nc_projections/default_nes.py +++ b/nes/nc_projections/default_nes.py @@ -10,6 +10,8 @@ from netCDF4 import Dataset, num2date, date2num from mpi4py import MPI from cfunits import Units from numpy.ma.core import MaskError +import geopandas as gpd +from shapely.geometry import Polygon from copy import deepcopy import datetime from ..interpolation import vertical_interpolation @@ -2061,6 +2063,55 @@ class Nes(object): return None + def to_shapefile(self, path): + """ + Write output file with shapefile format. + + Parameters + ---------- + path : str + Path to the output file. + """ + + if self._lat_bnds is None or self._lon_bnds is None: + self.create_spatial_bounds() + + # Reshape arrays to create geometry + aux_shape = (self.lat_bnds.shape[0], self.lon_bnds.shape[0], 4) + lon_bnds_aux = np.empty(aux_shape) + lon_bnds_aux[:, :, 0] = self.lon_bnds[np.newaxis, :, 0] + lon_bnds_aux[:, :, 1] = self.lon_bnds[np.newaxis, :, 1] + lon_bnds_aux[:, :, 2] = self.lon_bnds[np.newaxis, :, 1] + lon_bnds_aux[:, :, 3] = self.lon_bnds[np.newaxis, :, 0] + + lon_bnds = lon_bnds_aux + del lon_bnds_aux + + lat_bnds_aux = np.empty(aux_shape) + lat_bnds_aux[:, :, 0] = self.lat_bnds[:, np.newaxis, 0] + lat_bnds_aux[:, :, 1] = self.lat_bnds[:, np.newaxis, 0] + lat_bnds_aux[:, :, 2] = self.lat_bnds[:, np.newaxis, 1] + lat_bnds_aux[:, :, 3] = self.lat_bnds[:, np.newaxis, 1] + + lat_bnds = lat_bnds_aux + del lat_bnds_aux + + aux_b_lats = lat_bnds.reshape((lat_bnds.shape[0] * lat_bnds.shape[1], lat_bnds.shape[2])) + aux_b_lons = lon_bnds.reshape((lon_bnds.shape[0] * lon_bnds.shape[1], lon_bnds.shape[2])) + + # Create dataframe cointaining all polygons + geometry = [] + for i in range(aux_b_lons.shape[0]): + geometry.append(Polygon([(aux_b_lons[i, 0], aux_b_lats[i, 0]), + (aux_b_lons[i, 1], aux_b_lats[i, 1]), + (aux_b_lons[i, 2], aux_b_lats[i, 2]), + (aux_b_lons[i, 3], aux_b_lats[i, 3]), + (aux_b_lons[i, 0], aux_b_lats[i, 0])])) + gdf = gpd.GeoDataFrame(index=range(aux_b_lons.shape[0]), geometry=geometry, crs="EPSG:4326") + gdf.to_file(path) + + return None + def __gather_data_py_object(self): """ Gather all the variable data into the MPI rank 0 to perform a serial write. diff --git a/nes/nc_projections/lcc_nes.py b/nes/nc_projections/lcc_nes.py index 309767d..f6605a6 100644 --- a/nes/nc_projections/lcc_nes.py +++ b/nes/nc_projections/lcc_nes.py @@ -4,6 +4,8 @@ import numpy as np from cfunits import Units from pyproj import Proj from copy import deepcopy +import geopandas as gpd +from shapely.geometry import Polygon from .default_nes import Nes @@ -127,24 +129,26 @@ class LCCNes(Nes): Returns ------- - projection : dict + projection_data : dict Dictionary with the projection data. """ if create_nes: - projection = {'data': None, - 'dimensions': (), - 'grid_mapping_name': 'lambert_conformal_conic', - 'standard_parallel': [kwargs['lat_1'], kwargs['lat_2']], - 'longitude_of_central_meridian': kwargs['lon_0'], - 'latitude_of_projection_origin': kwargs['lat_0'], - } + projection_data = {'data': None, + 'dimensions': (), + 'grid_mapping_name': 'lambert_conformal_conic', + 'standard_parallel': [kwargs['lat_1'], kwargs['lat_2']], + 'longitude_of_central_meridian': kwargs['lon_0'], + 'latitude_of_projection_origin': kwargs['lat_0'], + } else: - projection = self.variables['Lambert_conformal'] + projection_data = self.variables['Lambert_conformal'] + projection_data['standard_parallel'] = [projection_data['standard_parallel'].split(', ')[0], + projection_data['standard_parallel'].split(', ')[1]] self.free_vars('Lambert_conformal') - return projection + return projection_data def _create_dimensions(self, netcdf): """ @@ -230,23 +234,23 @@ class LCCNes(Nes): x = np.array([self._x['data']] * len(self._y['data'])) y = np.array([self._y['data']] * len(self._x['data'])).T - projection = Proj( - proj='lcc', - ellps='WGS84', - R=6356752.3142, # WGS84_SEMIMINOR_AXIS (as defined in Cartopy source code) - lat_1=kwargs['lat_1'], - lat_2=kwargs['lat_2'], - lon_0=kwargs['lon_0'], - lat_0=kwargs['lat_0'], - to_meter=1, - x_0=0, - y_0=0, - a=6378137.0, # WGS84_SEMIMAJOR_AXIS (as defined in Cartopy source code) - k_0=1.0 - ) + self.projection = Proj( + proj='lcc', + ellps='WGS84', + R=6356752.3142, # WGS84_SEMIMINOR_AXIS (as defined in Cartopy source code) + lat_1=kwargs['lat_1'], + lat_2=kwargs['lat_2'], + lon_0=kwargs['lon_0'], + lat_0=kwargs['lat_0'], + to_meter=1, + x_0=0, + y_0=0, + a=6378137.0, # WGS84_SEMIMAJOR_AXIS (as defined in Cartopy source code) + k_0=1.0 + ) # Calculate centre latitudes and longitudes (UTM to LCC) - centre_lon_data, centre_lat_data = projection(x, y, inverse=True) + centre_lon_data, centre_lat_data = self.projection(x, y, inverse=True) centre_lat = {'data': centre_lat_data} centre_lon = {'data': centre_lon_data} @@ -289,39 +293,39 @@ class LCCNes(Nes): inc_y = np.abs(np.mean(np.diff(self.y['data']))) # Get bounds for rotated coordinates - y_bounds = self.create_single_spatial_bounds(self.y['data'], inc_y) - x_bounds = self.create_single_spatial_bounds(self.x['data'], inc_x) + y_bnds = self.create_single_spatial_bounds(self.y['data'], inc_y) + x_bnds = self.create_single_spatial_bounds(self.x['data'], inc_x) # Get rotated latitudes for grid edge - left_edge_y = np.append(y_bounds.flatten()[::2], y_bounds.flatten()[-1]) + left_edge_y = np.append(y_bnds.flatten()[::2], y_bnds.flatten()[-1]) right_edge_y = np.flip(left_edge_y, 0) - top_edge_y = np.repeat(y_bounds[-1][-1], len(self.x['data']) - 1) - bottom_edge_y = np.repeat(y_bounds[0][0], len(self.x['data'])) + top_edge_y = np.repeat(y_bnds[-1][-1], len(self.x['data']) - 1) + bottom_edge_y = np.repeat(y_bnds[0][0], len(self.x['data'])) y_grid_edge = np.concatenate((left_edge_y, top_edge_y, right_edge_y, bottom_edge_y)) # Get rotated longitudes for grid edge - left_edge_x = np.repeat(x_bounds[0][0], len(self.y['data']) + 1) - top_edge_x = x_bounds.flatten()[1:-1:2] - right_edge_x = np.repeat(x_bounds[-1][-1], len(self.y['data']) + 1) - bottom_edge_x = np.flip(x_bounds.flatten()[:-1:2], 0) + left_edge_x = np.repeat(x_bnds[0][0], len(self.y['data']) + 1) + top_edge_x = x_bnds.flatten()[1:-1:2] + right_edge_x = np.repeat(x_bnds[-1][-1], len(self.y['data']) + 1) + bottom_edge_x = np.flip(x_bnds.flatten()[:-1:2], 0) x_grid_edge = np.concatenate((left_edge_x, top_edge_x, right_edge_x, bottom_edge_x)) # Get edges for regular coordinates - projection = Proj( - proj='lcc', - ellps='WGS84', - R=6356752.3142, # WGS84_SEMIMINOR_AXIS (as defined in Cartopy source code) - lat_1=float(self.projection_data['standard_parallel'].split(', ')[0]), - lat_2=float(self.projection_data['standard_parallel'].split(', ')[1]), - lon_0=float(self.projection_data['longitude_of_central_meridian']), - lat_0=float(self.projection_data['latitude_of_projection_origin']), - to_meter=1, - x_0=0, - y_0=0, - a=6378137.0, # WGS84_SEMIMAJOR_AXIS (as defined in Cartopy source code) - k_0=1.0 - ) - grid_edge_lon_data, grid_edge_lat_data = projection(x_grid_edge, y_grid_edge, inverse=True) + self.projection = Proj( + proj='lcc', + ellps='WGS84', + R=6356752.3142, # WGS84_SEMIMINOR_AXIS (as defined in Cartopy source code) + lat_1=float(self.projection_data['standard_parallel'][0]), + lat_2=float(self.projection_data['standard_parallel'][1]), + lon_0=float(self.projection_data['longitude_of_central_meridian']), + lat_0=float(self.projection_data['latitude_of_projection_origin']), + to_meter=1, + x_0=0, + y_0=0, + a=6378137.0, # WGS84_SEMIMAJOR_AXIS (as defined in Cartopy source code) + k_0=1.0 + ) + grid_edge_lon_data, grid_edge_lat_data = self.projection(x_grid_edge, y_grid_edge, inverse=True) # Create grid outline by stacking the edges in both coordinates model_grid_outline = np.vstack((grid_edge_lon_data, grid_edge_lat_data)).T @@ -335,17 +339,38 @@ class LCCNes(Nes): Calculate longitude and latitude bounds and set them. """ - inc_lat = np.abs(np.mean(np.diff(self._lat['data']))) - lat_bnds = self.create_single_spatial_bounds(self._lat['data'], inc_lat, spatial_nv=4) - + # Calculate LCC coordinates bounds + inc_x = np.abs(np.mean(np.diff(self._x['data']))) + x_bnds = self.create_single_spatial_bounds(np.array([self._y['data']] * len(self._x['data'])).T, + inc_x, spatial_nv=4, inverse=True) + + inc_y = np.abs(np.mean(np.diff(self._y['data']))) + y_bnds = self.create_single_spatial_bounds(np.array([self._x['data']] * len(self._y['data'])), + inc_y, spatial_nv=4) + + # Transform LCC bounds to regular bounds + self.projection = Proj( + proj='lcc', + ellps='WGS84', + R=6356752.3142, # WGS84_SEMIMINOR_AXIS (as defined in Cartopy source code) + lat_1=float(self.projection_data['standard_parallel'][0]), + lat_2=float(self.projection_data['standard_parallel'][1]), + lon_0=float(self.projection_data['longitude_of_central_meridian']), + lat_0=float(self.projection_data['latitude_of_projection_origin']), + to_meter=1, + x_0=0, + y_0=0, + a=6378137.0, # WGS84_SEMIMAJOR_AXIS (as defined in Cartopy source code) + k_0=1.0 + ) + lon_bnds, lat_bnds = self.projection(x_bnds, y_bnds, inverse=True) + + # Obtain regular coordinates bounds self._lat_bnds = deepcopy(lat_bnds) self.lat_bnds = lat_bnds[self.write_axis_limits['y_min']:self.write_axis_limits['y_max'], self.write_axis_limits['x_min']:self.write_axis_limits['x_max'], :] - inc_lon = np.abs(np.mean(np.diff(self._lon['data']))) - lon_bnds = self.create_single_spatial_bounds(self._lon['data'], inc_lon, spatial_nv=4) - self._lon_bnds = deepcopy(lon_bnds) self.lon_bnds = lon_bnds[self.write_axis_limits['y_min']:self.write_axis_limits['y_max'], self.write_axis_limits['x_min']:self.write_axis_limits['x_max'], @@ -402,4 +427,35 @@ class LCCNes(Nes): Indicates if you want to print extra information during the process. """ - raise NotImplementedError("Grib2 format cannot be write in a Lambert Conformal Conic projection.") + raise NotImplementedError("Grib2 format cannot be written in a Lambert Conformal Conic projection.") + + def to_shapefile(self, path): + """ + Write output file with shapefile format. + + Parameters + ---------- + path : str + Path to the output file. + """ + + # Get latitude and longitude cell boundaries + if self._lat_bnds is None or self._lon_bnds is None: + self.create_spatial_bounds() + + # Reshape arrays to create geometry + aux_b_lats = self.lat_bnds.reshape((self.lat_bnds.shape[0] * self.lat_bnds.shape[1], self.lat_bnds.shape[2])) + aux_b_lons = self.lon_bnds.reshape((self.lon_bnds.shape[0] * self.lon_bnds.shape[1], self.lon_bnds.shape[2])) + + # Create dataframe cointaining all polygons + geometry = [] + for i in range(aux_b_lons.shape[0]): + geometry.append(Polygon([(aux_b_lons[i, 0], aux_b_lats[i, 0]), + (aux_b_lons[i, 1], aux_b_lats[i, 1]), + (aux_b_lons[i, 2], aux_b_lats[i, 2]), + (aux_b_lons[i, 3], aux_b_lats[i, 3]), + (aux_b_lons[i, 0], aux_b_lats[i, 0])])) + gdf = gpd.GeoDataFrame(index=range(aux_b_lons.shape[0]), geometry=geometry, crs="EPSG:4326") + gdf.to_file(path) + + return None diff --git a/nes/nc_projections/mercator_nes.py b/nes/nc_projections/mercator_nes.py index 43de48e..e5272aa 100644 --- a/nes/nc_projections/mercator_nes.py +++ b/nes/nc_projections/mercator_nes.py @@ -4,6 +4,8 @@ import numpy as np from cfunits import Units from pyproj import Proj from copy import deepcopy +import geopandas as gpd +from shapely.geometry import Polygon from nes.nc_projections.default_nes import Nes @@ -88,7 +90,6 @@ class MercatorNes(Nes): self.free_vars('crs') - @staticmethod def new(comm=None, path=None, info=False, dataset=None, xarray=False, create=False, balanced=False, parallel_method='Y', avoid_first_hours=0, avoid_last_hours=0, first_level=0, last_level=None): @@ -128,23 +129,23 @@ class MercatorNes(Nes): Returns ------- - projection : dict + projection _data: dict Dictionary with the projection data. """ if create_nes: - projection = {'data': None, - 'dimensions': (), - 'grid_mapping_name': 'mercator', - 'standard_parallel': [kwargs['lat_ts']], # TODO: Check if True - 'longitude_of_projection_origin': kwargs['lon_0'], - } + projection_data = {'data': None, + 'dimensions': (), + 'grid_mapping_name': 'mercator', + 'standard_parallel': [kwargs['lat_ts']], # TODO: Check if True + 'longitude_of_projection_origin': kwargs['lon_0'], + } else: - projection = self.variables['mercator'] + projection_data = self.variables['mercator'] self.free_vars('mercator') - return projection + return projection_data def _create_dimensions(self, netcdf): """ @@ -232,16 +233,16 @@ class MercatorNes(Nes): x = np.array([self._x['data']] * len(self._y['data'])) y = np.array([self._y['data']] * len(self._x['data'])).T - projection = Proj( - proj='merc', - a=6378137.0, # WGS84_SEMIMAJOR_AXIS (as defined in Cartopy source code) - b=6356752.3142, # WGS84_SEMIMINOR_AXIS (as defined in Cartopy source code) - lat_ts=kwargs['lat_ts'], - lon_0=kwargs['lon_0'], - ) + self.projection = Proj( + proj='merc', + a=6378137.0, # WGS84_SEMIMAJOR_AXIS (as defined in Cartopy source code) + b=6356752.3142, # WGS84_SEMIMINOR_AXIS (as defined in Cartopy source code) + lat_ts=kwargs['lat_ts'], + lon_0=kwargs['lon_0'], + ) # Calculate centre latitudes and longitudes (UTM to Mercator) - centre_lon_data, centre_lat_data = projection(x, y, inverse=True) + centre_lon_data, centre_lat_data = self.projection(x, y, inverse=True) centre_lat = {'data': centre_lat_data} centre_lon = {'data': centre_lon_data} @@ -302,14 +303,14 @@ class MercatorNes(Nes): x_grid_edge = np.concatenate((left_edge_x, top_edge_x, right_edge_x, bottom_edge_x)) # Get edges for regular coordinates - projection = Proj( - proj='merc', - a=6378137.0, # WGS84_SEMIMAJOR_AXIS (as defined in Cartopy source code) - b=6356752.3142, # WGS84_SEMIMINOR_AXIS (as defined in Cartopy source code) - lat_ts=float(self.projection_data['standard_parallel']), - lon_0=float(self.projection_data['longitude_of_central_meridian']), - ) - grid_edge_lon_data, grid_edge_lat_data = projection(x_grid_edge, y_grid_edge, inverse=True) + self.projection = Proj( + proj='merc', + a=6378137.0, # WGS84_SEMIMAJOR_AXIS (as defined in Cartopy source code) + b=6356752.3142, # WGS84_SEMIMINOR_AXIS (as defined in Cartopy source code) + lat_ts=float(self.projection_data['standard_parallel'][0]), + lon_0=float(self.projection_data['longitude_of_projection_origin']), + ) + grid_edge_lon_data, grid_edge_lat_data = self.projection(x_grid_edge, y_grid_edge, inverse=True) # Create grid outline by stacking the edges in both coordinates model_grid_outline = np.vstack((grid_edge_lon_data, grid_edge_lat_data)).T @@ -323,17 +324,31 @@ class MercatorNes(Nes): Calculate longitude and latitude bounds and set them. """ - inc_lat = np.abs(np.mean(np.diff(self._lat['data']))) - lat_bnds = self.create_single_spatial_bounds(self._lat['data'], inc_lat, spatial_nv=4) - + # Calculate Mercator coordinates bounds + inc_x = np.abs(np.mean(np.diff(self._x['data']))) + x_bnds = self.create_single_spatial_bounds(np.array([self._y['data']] * len(self._x['data'])).T, + inc_x, spatial_nv=4, inverse=True) + + inc_y = np.abs(np.mean(np.diff(self._y['data']))) + y_bnds = self.create_single_spatial_bounds(np.array([self._x['data']] * len(self._y['data'])), + inc_y, spatial_nv=4) + + # Transform Mercator bounds to regular bounds + self.projection = Proj( + proj='merc', + a=6378137.0, # WGS84_SEMIMAJOR_AXIS (as defined in Cartopy source code) + b=6356752.3142, # WGS84_SEMIMINOR_AXIS (as defined in Cartopy source code) + lat_ts=float(self.projection_data['standard_parallel'][0]), + lon_0=float(self.projection_data['longitude_of_projection_origin']), + ) + lon_bnds, lat_bnds = self.projection(x_bnds, y_bnds, inverse=True) + + # Obtain regular coordinates bounds self._lat_bnds = deepcopy(lat_bnds) self.lat_bnds = lat_bnds[self.write_axis_limits['y_min']:self.write_axis_limits['y_max'], self.write_axis_limits['x_min']:self.write_axis_limits['x_max'], :] - inc_lon = np.abs(np.mean(np.diff(self._lon['data']))) - lon_bnds = self.create_single_spatial_bounds(self._lon['data'], inc_lon, spatial_nv=4) - self._lon_bnds = deepcopy(lon_bnds) self.lon_bnds = lon_bnds[self.write_axis_limits['y_min']:self.write_axis_limits['y_max'], self.write_axis_limits['x_min']:self.write_axis_limits['x_max'], @@ -389,4 +404,35 @@ class MercatorNes(Nes): Indicates if you want to print extra information during the process. """ - raise NotImplementedError("Grib2 format cannot be write in a Mercator projection.") + raise NotImplementedError("Grib2 format cannot be written in a Mercator projection.") + + def to_shapefile(self, path): + """ + Write output file with shapefile format. + + Parameters + ---------- + path : str + Path to the output file. + """ + + # Get latitude and longitude cell boundaries + if self._lat_bnds is None or self._lon_bnds is None: + self.create_spatial_bounds() + + # Reshape arrays to create geometry + aux_b_lats = self.lat_bnds.reshape((self.lat_bnds.shape[0] * self.lat_bnds.shape[1], self.lat_bnds.shape[2])) + aux_b_lons = self.lon_bnds.reshape((self.lon_bnds.shape[0] * self.lon_bnds.shape[1], self.lon_bnds.shape[2])) + + # Create dataframe cointaining all polygons + geometry = [] + for i in range(aux_b_lons.shape[0]): + geometry.append(Polygon([(aux_b_lons[i, 0], aux_b_lats[i, 0]), + (aux_b_lons[i, 1], aux_b_lats[i, 1]), + (aux_b_lons[i, 2], aux_b_lats[i, 2]), + (aux_b_lons[i, 3], aux_b_lats[i, 3]), + (aux_b_lons[i, 0], aux_b_lats[i, 0])])) + gdf = gpd.GeoDataFrame(index=range(aux_b_lons.shape[0]), geometry=geometry, crs="EPSG:4326") + gdf.to_file(path) + + return None diff --git a/nes/nc_projections/points_nes.py b/nes/nc_projections/points_nes.py index 3c0769f..a647cc0 100644 --- a/nes/nc_projections/points_nes.py +++ b/nes/nc_projections/points_nes.py @@ -573,24 +573,6 @@ class PointsNes(Nes): raise NotImplementedError("Spatial bounds cannot be created for points datasets.") - def to_grib2(self, path, grib_keys, grib_template_path, info=False): - """ - Write output file with grib2 format. - - Parameters - ---------- - path : str - Path to the output file. - grib_keys : dict - Dictionary with the grib2 keys. - grib_template_path : str - Path to the grib2 file to use as template. - info : bool - Indicates if you want to print extra information during the process. - """ - - raise NotImplementedError("Grib2 format cannot be write with point data.") - def to_providentia(self, model_centre_lon, model_centre_lat, grid_edge_lon, grid_edge_lat): """ Transform a PointsNes into a PointsNesProvidentia object @@ -639,3 +621,33 @@ class PointsNes(Nes): points_nes_providentia.variables = variables return points_nes_providentia + + def to_grib2(self, path, grib_keys, grib_template_path, info=False): + """ + Write output file with grib2 format. + + Parameters + ---------- + path : str + Path to the output file. + grib_keys : dict + Dictionary with the grib2 keys. + grib_template_path : str + Path to the grib2 file to use as template. + info : bool + Indicates if you want to print extra information during the process. + """ + + raise NotImplementedError("Grib2 format cannot be written with point data.") + + def to_shapefile(self, path): + """ + Write output file with shapefile format. + + Parameters + ---------- + path : str + Path to the output file. + """ + + raise NotImplementedError("Shapefile format cannot be written with point data.") diff --git a/nes/nc_projections/rotated_nes.py b/nes/nc_projections/rotated_nes.py index d87f7e5..c5b2045 100644 --- a/nes/nc_projections/rotated_nes.py +++ b/nes/nc_projections/rotated_nes.py @@ -4,6 +4,8 @@ import numpy as np import math from cfunits import Units from copy import deepcopy +import geopandas as gpd +from shapely.geometry import Polygon from .default_nes import Nes @@ -59,7 +61,7 @@ class RotatedNes(Nes): avoid_last_hours : int Number of hours to remove from last time steps. """ - + super(RotatedNes, self).__init__(comm=comm, path=path, info=info, dataset=dataset, balanced=balanced, xarray=xarray, parallel_method=parallel_method, @@ -128,23 +130,23 @@ class RotatedNes(Nes): Returns ------- - projection : dict + projection_data : dict Dictionary with the projection data. """ if create_nes: - projection = {'data': None, - 'dimensions': (), - 'grid_mapping_name': 'rotated_latitude_longitude', - 'grid_north_pole_latitude': 90 - kwargs['centre_lat'], - 'grid_north_pole_longitude': -180 + kwargs['centre_lon'], - } + projection_data = {'data': None, + 'dimensions': (), + 'grid_mapping_name': 'rotated_latitude_longitude', + 'grid_north_pole_latitude': 90 - kwargs['centre_lat'], + 'grid_north_pole_longitude': -180 + kwargs['centre_lon'], + } else: - projection = self.variables['rotated_pole'] + projection_data = self.variables['rotated_pole'] self.free_vars('rotated_pole') - return projection + return projection_data def _create_dimensions(self, netcdf): """ @@ -379,17 +381,24 @@ class RotatedNes(Nes): Calculate longitude and latitude bounds and set them. """ - inc_lat = np.abs(np.mean(np.diff(self._lat['data']))) - lat_bnds = self.create_single_spatial_bounds(self._lat['data'], inc_lat, spatial_nv=4) + # Calculate rotated coordinates bounds + inc_rlat = np.abs(np.mean(np.diff(self._rlat['data']))) + rlat_bnds = self.create_single_spatial_bounds(np.array([self._rlat['data']] * len(self._rlon['data'])).T, + inc_rlat, spatial_nv=4, inverse=True) + + inc_rlon = np.abs(np.mean(np.diff(self._rlon['data']))) + rlon_bnds = self.create_single_spatial_bounds(np.array([self._rlon['data']] * len(self._rlat['data'])), + inc_rlon, spatial_nv=4) + # Transform rotated bounds to regular bounds + lon_bnds, lat_bnds = self.rotated2latlon(rlon_bnds, rlat_bnds) + + # Obtain regular coordinates bounds self._lat_bnds = deepcopy(lat_bnds) self.lat_bnds = lat_bnds[self.write_axis_limits['y_min']:self.write_axis_limits['y_max'], self.write_axis_limits['x_min']:self.write_axis_limits['x_max'], :] - inc_lon = np.abs(np.mean(np.diff(self._lon['data']))) - lon_bnds = self.create_single_spatial_bounds(self._lon['data'], inc_lon, spatial_nv=4) - self._lon_bnds = deepcopy(lon_bnds) self.lon_bnds = lon_bnds[self.write_axis_limits['y_min']:self.write_axis_limits['y_max'], self.write_axis_limits['x_min']:self.write_axis_limits['x_max'], @@ -445,4 +454,34 @@ class RotatedNes(Nes): Indicates if you want to print extra information during the process. """ - raise NotImplementedError("Grib2 format cannot be write in a Rotated pole projection.") + raise NotImplementedError("Grib2 format cannot be written in a Rotated pole projection.") + + def to_shapefile(self, path): + """ + Write output file with shapefile format. + + Parameters + ---------- + path : str + Path to the output file. + """ + + if self._lat_bnds is None or self._lon_bnds is None: + self.create_spatial_bounds() + + # Reshape arrays to create geometry + aux_b_lats = self.lat_bnds.reshape((self.lat_bnds.shape[0] * self.lat_bnds.shape[1], self.lat_bnds.shape[2])) + aux_b_lons = self.lon_bnds.reshape((self.lon_bnds.shape[0] * self.lon_bnds.shape[1], self.lon_bnds.shape[2])) + + # Create dataframe cointaining all polygons + geometry = [] + for i in range(aux_b_lons.shape[0]): + geometry.append(Polygon([(aux_b_lons[i, 0], aux_b_lats[i, 0]), + (aux_b_lons[i, 1], aux_b_lats[i, 1]), + (aux_b_lons[i, 2], aux_b_lats[i, 2]), + (aux_b_lons[i, 3], aux_b_lats[i, 3]), + (aux_b_lons[i, 0], aux_b_lats[i, 0])])) + gdf = gpd.GeoDataFrame(index=range(aux_b_lons.shape[0]), geometry=geometry, crs="EPSG:4326") + gdf.to_file(path) + + return None diff --git a/tests/scalability_test_nord3v2.bash b/tests/scalability_test_nord3v2.bash index 7e39810..c725675 100644 --- a/tests/scalability_test_nord3v2.bash +++ b/tests/scalability_test_nord3v2.bash @@ -15,7 +15,8 @@ module load mpi4py/3.0.3-foss-2019b-Python-3.7.4 module load filelock/3.7.1-foss-2019b-Python-3.7.4 module load pyproj/2.5.0-foss-2019b-Python-3.7.4 module load eccodes-python/0.9.5-foss-2019b-Python-3.7.4 - +module load geopandas/0.8.1-foss-2019b-Python-3.7.4 +module load Shapely/1.7.1-foss-2019b-Python-3.7.4 for nprocs in 1 2 4 8 do diff --git a/tests/test_bash_mn4.cmd b/tests/test_bash_mn4.cmd index 1aa487e..db2758f 100644 --- a/tests/test_bash_mn4.cmd +++ b/tests/test_bash_mn4.cmd @@ -23,7 +23,8 @@ module load OpenMPI/4.0.5-GCC-8.3.0-mn4 module load filelock/3.7.1-foss-2019b-Python-3.7.4 module load eccodes-python/0.9.5-foss-2019b-Python-3.7.4 module load pyproj/2.5.0-foss-2019b-Python-3.7.4 - +module load geopandas/0.8.1-foss-2019b-Python-3.7.4 +module load Shapely/1.7.1-foss-2019b-Python-3.7.4 export PYTHONPATH=/gpfs/projects/bsc32/models/NES:${PYTHONPATH} cd /gpfs/projects/bsc32/models/NES/tests diff --git a/tests/test_bash_nord3v2.cmd b/tests/test_bash_nord3v2.cmd index c177698..da0f4d7 100644 --- a/tests/test_bash_nord3v2.cmd +++ b/tests/test_bash_nord3v2.cmd @@ -23,6 +23,8 @@ module load mpi4py/3.0.3-foss-2019b-Python-3.7.4 module load filelock/3.7.1-foss-2019b-Python-3.7.4 module load eccodes-python/0.9.5-foss-2019b-Python-3.7.4 module load pyproj/2.5.0-foss-2019b-Python-3.7.4 +module load geopandas/0.8.1-foss-2019b-Python-3.7.4 +module load Shapely/1.7.1-foss-2019b-Python-3.7.4 export PYTHONPATH=/gpfs/projects/bsc32/models/NES:${PYTHONPATH} cd /gpfs/projects/bsc32/models/NES/tests diff --git a/tutorials/1.Introduction/1.5.Read_Write_Mercator.ipynb b/tutorials/1.Introduction/1.5.Read_Write_Mercator.ipynb index 67f9bb4..d4cb3f3 100644 --- a/tutorials/1.Introduction/1.5.Read_Write_Mercator.ipynb +++ b/tutorials/1.Introduction/1.5.Read_Write_Mercator.ipynb @@ -413,7 +413,7 @@ " lon_bnds (y, x, nv) float32 ...\n", " var_aux (time, y, x) float32 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0\n", " mercator int32 -2147483647\n", - " cell_area (y, x) float32 1.316e+09 1.316e+09 ... 1.051e+09 1.051e+09
    • " ], "text/plain": [ "\n", @@ -498,7 +498,7 @@ { "data": { "text/plain": [ - "" + "" ] }, "execution_count": 4, @@ -1202,91 +1202,6 @@ "text": [ "Rank 000: Creating mercator_file_1.nc\n", "Rank 000: NetCDF ready to write\n", - "[[[-43.68109893798828 -43.68109893798828 -43.354862213134766\n", - " -43.354862213134766]\n", - " [-43.68109893798828 -43.68109893798828 -43.354862213134766\n", - " -43.354862213134766]\n", - " [-43.68109893798828 -43.68109893798828 -43.354862213134766\n", - " -43.354862213134766]\n", - " ...\n", - " [-43.68109893798828 -43.68109893798828 -43.354862213134766\n", - " -43.354862213134766]\n", - " [-43.68109893798828 -43.68109893798828 -43.354862213134766\n", - " -43.354862213134766]\n", - " [-43.68109893798828 -43.68109893798828 -43.354862213134766\n", - " -43.354862213134766]]\n", - "\n", - " [[-43.354862213134766 -43.354862213134766 -43.02686309814453\n", - " -43.02686309814453]\n", - " [-43.354862213134766 -43.354862213134766 -43.02686309814453\n", - " -43.02686309814453]\n", - " [-43.354862213134766 -43.354862213134766 -43.02686309814453\n", - " -43.02686309814453]\n", - " ...\n", - " [-43.354862213134766 -43.354862213134766 -43.02686309814453\n", - " -43.02686309814453]\n", - " [-43.354862213134766 -43.354862213134766 -43.02686309814453\n", - " -43.02686309814453]\n", - " [-43.354862213134766 -43.354862213134766 -43.02686309814453\n", - " -43.02686309814453]]\n", - "\n", - " [[-43.02686309814453 -43.02686309814453 -42.69709777832031\n", - " -42.69709777832031]\n", - " [-43.02686309814453 -43.02686309814453 -42.69709777832031\n", - " -42.69709777832031]\n", - " [-43.02686309814453 -43.02686309814453 -42.69709777832031\n", - " -42.69709777832031]\n", - " ...\n", - " [-43.02686309814453 -43.02686309814453 -42.69709777832031\n", - " -42.69709777832031]\n", - " [-43.02686309814453 -43.02686309814453 -42.69709777832031\n", - " -42.69709777832031]\n", - " [-43.02686309814453 -43.02686309814453 -42.69709777832031\n", - " -42.69709777832031]]\n", - "\n", - " ...\n", - "\n", - " [[48.86896514892578 48.86896514892578 49.16401672363281\n", - " 49.16401672363281]\n", - " [48.86896514892578 48.86896514892578 49.16401672363281\n", - " 49.16401672363281]\n", - " [48.86896514892578 48.86896514892578 49.16401672363281\n", - " 49.16401672363281]\n", - " ...\n", - " [48.86896514892578 48.86896514892578 49.16401672363281\n", - " 49.16401672363281]\n", - " [48.86896514892578 48.86896514892578 49.16401672363281\n", - " 49.16401672363281]\n", - " [48.86896514892578 48.86896514892578 49.16401672363281\n", - " 49.16401672363281]]\n", - "\n", - " [[49.16401672363281 49.16401672363281 49.45732498168945\n", - " 49.45732498168945]\n", - " [49.16401672363281 49.16401672363281 49.45732498168945\n", - " 49.45732498168945]\n", - " [49.16401672363281 49.16401672363281 49.45732498168945\n", - " 49.45732498168945]\n", - " ...\n", - " [49.16401672363281 49.16401672363281 49.45732498168945\n", - " 49.45732498168945]\n", - " [49.16401672363281 49.16401672363281 49.45732498168945\n", - " 49.45732498168945]\n", - " [49.16401672363281 49.16401672363281 49.45732498168945\n", - " 49.45732498168945]]\n", - "\n", - " [[49.45732498168945 49.45732498168945 49.74888229370117\n", - " 49.74888229370117]\n", - " [49.45732498168945 49.45732498168945 49.74888229370117\n", - " 49.74888229370117]\n", - " [49.45732498168945 49.45732498168945 49.74888229370117\n", - " 49.74888229370117]\n", - " ...\n", - " [49.45732498168945 49.45732498168945 49.74888229370117\n", - " 49.74888229370117]\n", - " [49.45732498168945 49.45732498168945 49.74888229370117\n", - " 49.74888229370117]\n", - " [49.45732498168945 49.45732498168945 49.74888229370117\n", - " 49.74888229370117]]]\n", "Rank 000: Dimensions done\n", "Rank 000: Writing var_aux var (1/2)\n", "Rank 000: Var var_aux created (1/2)\n", @@ -1320,7 +1235,7 @@ { "data": { "text/plain": [ - "" + "" ] }, "execution_count": 16, @@ -1716,7 +1631,7 @@ " cell_area (time, lev, y, x) float32 1.316e+09 1.316e+09 ... 1.051e+09\n", " mercator |S1 b''\n", "Attributes:\n", - " Conventions: CF-1.7
  • Conventions :
    CF-1.7
    • " ], "text/plain": [ "\n", diff --git a/tutorials/2.Creation/2.2.Create_Rotated.ipynb b/tutorials/2.Creation/2.2.Create_Rotated.ipynb index a26bdea..5493bd5 100644 --- a/tutorials/2.Creation/2.2.Create_Rotated.ipynb +++ b/tutorials/2.Creation/2.2.Create_Rotated.ipynb @@ -9,7 +9,7 @@ }, { "cell_type": "code", - "execution_count": 5, + "execution_count": 1, "metadata": {}, "outputs": [], "source": [ @@ -19,7 +19,7 @@ }, { "cell_type": "code", - "execution_count": 6, + "execution_count": 2, "metadata": {}, "outputs": [], "source": [ @@ -37,7 +37,7 @@ }, { "cell_type": "code", - "execution_count": 7, + "execution_count": 3, "metadata": {}, "outputs": [ { @@ -56,7 +56,7 @@ }, { "cell_type": "code", - "execution_count": 9, + "execution_count": 4, "metadata": {}, "outputs": [ { @@ -425,13 +425,13 @@ " lon (rlat, rlon) float64 -22.18 -22.02 -21.85 ... 88.05 88.23\n", " rotated_pole |S1 b''\n", "Attributes:\n", - " Conventions: CF-1.7
  • Conventions :
    CF-1.7
    • " ], "text/plain": [ "\n", @@ -461,7 +461,7 @@ " Conventions: CF-1.7" ] }, - "execution_count": 9, + "execution_count": 4, "metadata": {}, "output_type": "execute_result" } diff --git a/tutorials/2.Creation/2.5.Create-LCC.ipynb b/tutorials/2.Creation/2.5.Create-LCC.ipynb index 3f9fd7d..e64d520 100644 --- a/tutorials/2.Creation/2.5.Create-LCC.ipynb +++ b/tutorials/2.Creation/2.5.Create-LCC.ipynb @@ -9,7 +9,7 @@ }, { "cell_type": "code", - "execution_count": 2, + "execution_count": 1, "metadata": {}, "outputs": [], "source": [ @@ -19,7 +19,7 @@ }, { "cell_type": "code", - "execution_count": 3, + "execution_count": 2, "metadata": {}, "outputs": [], "source": [ @@ -40,7 +40,7 @@ }, { "cell_type": "code", - "execution_count": 4, + "execution_count": 3, "metadata": {}, "outputs": [ { @@ -59,7 +59,7 @@ }, { "cell_type": "code", - "execution_count": 5, + "execution_count": 4, "metadata": {}, "outputs": [ { @@ -428,9 +428,9 @@ " lon (y, x) float64 ...\n", " Lambert_conformal |S1 b''\n", "Attributes:\n", - " Conventions: CF-1.7" + " Conventions: CF-1.7" ], "text/plain": [ "\n", @@ -448,7 +448,7 @@ " Conventions: CF-1.7" ] }, - "execution_count": 5, + "execution_count": 4, "metadata": {}, "output_type": "execute_result" } diff --git a/tutorials/2.Creation/2.6.Create_Mercator.ipynb b/tutorials/2.Creation/2.6.Create_Mercator.ipynb index d53476d..e2c7399 100644 --- a/tutorials/2.Creation/2.6.Create_Mercator.ipynb +++ b/tutorials/2.Creation/2.6.Create_Mercator.ipynb @@ -9,7 +9,7 @@ }, { "cell_type": "code", - "execution_count": 9, + "execution_count": 1, "metadata": {}, "outputs": [], "source": [ @@ -19,7 +19,7 @@ }, { "cell_type": "code", - "execution_count": 10, + "execution_count": 2, "metadata": {}, "outputs": [], "source": [ @@ -38,7 +38,7 @@ }, { "cell_type": "code", - "execution_count": 11, + "execution_count": 3, "metadata": {}, "outputs": [ { @@ -57,7 +57,7 @@ }, { "cell_type": "code", - "execution_count": 12, + "execution_count": 4, "metadata": {}, "outputs": [ { @@ -422,36 +422,36 @@ " * y (y) float64 -5.382e+06 -5.332e+06 ... 6.318e+06 6.368e+06\n", " * x (x) float64 -1.01e+05 -5.102e+04 -1.018e+03 ... 1.03e+07 1.035e+07\n", "Data variables:\n", - " lat (y, x) float64 -43.52 -43.52 -43.52 -43.52 ... 49.6 49.6 49.6 49.6\n", - " lon (y, x) float64 -18.91 -18.46 -18.01 -17.56 ... 74.22 74.67 75.12\n", + " lat (y, x) float64 -43.67 -43.67 -43.67 -43.67 ... 49.75 49.75 49.75\n", + " lon (y, x) float64 -18.91 -18.46 -18.01 -17.56 ... 74.1 74.55 75.0\n", " mercator |S1 b''\n", "Attributes:\n", - " Conventions: CF-1.7
  • Conventions :
    CF-1.7
    • " ], "text/plain": [ "\n", @@ -469,7 +469,7 @@ " Conventions: CF-1.7" ] }, - "execution_count": 12, + "execution_count": 4, "metadata": {}, "output_type": "execute_result" } diff --git a/tutorials/4.Interpolation/4.3.Providentia_Interpolation.ipynb b/tutorials/4.Interpolation/4.3.Providentia_Interpolation.ipynb index 067c5f3..9c9e928 100644 --- a/tutorials/4.Interpolation/4.3.Providentia_Interpolation.ipynb +++ b/tutorials/4.Interpolation/4.3.Providentia_Interpolation.ipynb @@ -26,7 +26,7 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 2, "metadata": {}, "outputs": [], "source": [ @@ -41,7 +41,7 @@ { "data": { "text/plain": [ - "" + "" ] }, "execution_count": 3, @@ -239,7 +239,7 @@ { "data": { "text/plain": [ - "" + "" ] }, "execution_count": 8, @@ -607,7 +607,7 @@ { "data": { "text/plain": [ - "" + "" ] }, "execution_count": 15, @@ -700,7 +700,7 @@ "name": "stderr", "output_type": "stream", "text": [ - "/esarchive/scratch/avilanova/software/NES/nes/nc_projections/points_nes_providentia.py:581: UserWarning: WARNING!!! Providentia datasets cannot be written in parallel yet. Changing to serial mode.\n", + "/esarchive/scratch/avilanova/software/NES/nes/nc_projections/points_nes_providentia.py:585: UserWarning: WARNING!!! Providentia datasets cannot be written in parallel yet. Changing to serial mode.\n", " warnings.warn(msg)\n" ] } @@ -1083,10 +1083,10 @@ " grid_edge_latitude (grid_edge) float64 29.9 30.1 30.3 ... 29.9 29.9\n", " sconco3 (station, time) float64 ...\n", "Attributes:\n", - " Conventions: CF-1.7
  • Conventions :
    CF-1.7
    • " ], "text/plain": [ "\n", @@ -1246,7 +1246,7 @@ { "data": { "text/plain": [ - "" + "" ] }, "execution_count": 22, @@ -1731,11 +1731,11 @@ " conventions: CF-1.7\n", " data_version: 1.0\n", " history: Thu Feb 11 10:19:01 2021: ncks -O --dfl_lvl 1 /gpfs/proje...\n", - " NCO: 4.7.2
  • title :
    Inverse distance weighting (4 neighbours) interpolated cams61_chimere_ph2 experiment data for the component sconco3 with reference to the measurement stations in the EBAS network in 2018-04.
    institution :
    Barcelona Supercomputing Center
    source :
    Experiment cams61_chimere_ph2
    creator_name :
    Dene R. Bowdalo
    creator_email :
    dene.bowdalo@bsc.es
    conventions :
    CF-1.7
    data_version :
    1.0
    history :
    Thu Feb 11 10:19:01 2021: ncks -O --dfl_lvl 1 /gpfs/projects/bsc32/AC_cache/recon/exp_interp/1.3.3/cams61_chimere_ph2-eu-000/hourly/sconco3/EBAS/sconco3_201804.nc /gpfs/projects/bsc32/AC_cache/recon/exp_interp/1.3.3/cams61_chimere_ph2-eu-000/hourly/sconco3/EBAS/sconco3_201804.nc
    NCO :
    4.7.2
    • " ], "text/plain": [ "\n", diff --git a/tutorials/5.Others/5.2.Add_Coordinates_Bounds.ipynb b/tutorials/5.Others/5.2.Add_Coordinates_Bounds.ipynb index 8ad4167..29ca013 100644 --- a/tutorials/5.Others/5.2.Add_Coordinates_Bounds.ipynb +++ b/tutorials/5.Others/5.2.Add_Coordinates_Bounds.ipynb @@ -78,58 +78,55 @@ { "data": { "text/plain": [ - "masked_array(\n", - " data=[[[16.35033798, 16.35033798, 16.35033798, 16.35033798],\n", - " [16.43292999, 16.43292999, 16.43292999, 16.43292999],\n", - " [16.51514626, 16.51514626, 16.51514626, 16.51514626],\n", - " ...,\n", - " [16.51514626, 16.51514626, 16.51514626, 16.51514626],\n", - " [16.43292999, 16.43292999, 16.43292999, 16.43292999],\n", - " [16.35033798, 16.35033798, 16.35033798, 16.35033798]],\n", - "\n", - " [[16.52742577, 16.52742577, 16.52742577, 16.52742577],\n", - " [16.61023903, 16.61023903, 16.61023903, 16.61023903],\n", - " [16.69267654, 16.69267654, 16.69267654, 16.69267654],\n", - " ...,\n", - " [16.69267654, 16.69267654, 16.69267654, 16.69267654],\n", - " [16.61024284, 16.61024284, 16.61024284, 16.61024284],\n", - " [16.52742577, 16.52742577, 16.52742577, 16.52742577]],\n", + "array([[[16.2203979 , 16.30306824, 16.48028979, 16.39739715],\n", + " [16.30306855, 16.3853609 , 16.56280424, 16.48029011],\n", + " [16.38536121, 16.46727425, 16.64493885, 16.56280455],\n", + " ...,\n", + " [16.46727269, 16.38535964, 16.56280298, 16.64493728],\n", + " [16.3853609 , 16.30306855, 16.48029011, 16.56280424],\n", + " [16.30306824, 16.2203979 , 16.39739715, 16.48028979]],\n", "\n", - " [[16.70447159, 16.70447159, 16.70447159, 16.70447159],\n", - " [16.78750801, 16.78750801, 16.78750801, 16.78750801],\n", - " [16.87016678, 16.87016678, 16.87016678, 16.87016678],\n", - " ...,\n", - " [16.87016678, 16.87016678, 16.87016678, 16.87016678],\n", - " [16.78750992, 16.78750992, 16.78750992, 16.78750992],\n", - " [16.70447159, 16.70447159, 16.70447159, 16.70447159]],\n", + " [[16.39739783, 16.48029047, 16.65746762, 16.57435251],\n", + " [16.48029079, 16.56280491, 16.74020402, 16.65746794],\n", + " [16.56280523, 16.64493952, 16.82256006, 16.74020434],\n", + " ...,\n", + " [16.64493796, 16.56280366, 16.74020276, 16.82255849],\n", + " [16.56280491, 16.48029079, 16.65746794, 16.74020402],\n", + " [16.48029047, 16.39739783, 16.57435251, 16.65746762]],\n", "\n", + " [[16.57435149, 16.65746661, 16.83459876, 16.751261 ],\n", + " [16.65746692, 16.74020301, 16.91755729, 16.83459908],\n", + " [16.74020332, 16.82255904, 17.00013494, 16.91755761],\n", " ...,\n", + " [16.82255748, 16.74020175, 16.91755603, 17.00013337],\n", + " [16.74020301, 16.65746692, 16.83459908, 16.91755729],\n", + " [16.65746661, 16.57435149, 16.751261 , 16.83459876]],\n", "\n", - " [[58.32094955, 58.32094955, 58.32094955, 58.32094955],\n", - " [58.47268295, 58.47268295, 58.47268295, 58.47268295],\n", - " [58.62430954, 58.62430954, 58.62430954, 58.62430954],\n", - " ...,\n", - " [58.62430954, 58.62430954, 58.62430954, 58.62430954],\n", - " [58.47268295, 58.47268295, 58.47268295, 58.47268295],\n", - " [58.32094955, 58.32094955, 58.32094955, 58.32094955]],\n", + " ...,\n", "\n", - " [[58.42628479, 58.42628479, 58.42628479, 58.42628479],\n", - " [58.57820129, 58.57820129, 58.57820129, 58.57820129],\n", - " [58.73002625, 58.73002625, 58.73002625, 58.73002625],\n", - " ...,\n", - " [58.73002625, 58.73002625, 58.73002625, 58.73002625],\n", - " [58.57820129, 58.57820129, 58.57820129, 58.57820129],\n", - " [58.42628479, 58.42628479, 58.42628479, 58.42628479]],\n", + " [[58.19210948, 58.34380497, 58.44964444, 58.29776032],\n", + " [58.34380555, 58.49539321, 58.6014247 , 58.44964502],\n", + " [58.49539378, 58.64687141, 58.75309835, 58.60142528],\n", + " ...,\n", + " [58.64686852, 58.49539089, 58.60142239, 58.75309546],\n", + " [58.49539321, 58.34380555, 58.44964502, 58.6014247 ],\n", + " [58.34380497, 58.19210948, 58.29776032, 58.44964444]],\n", "\n", - " [[58.53079224, 58.53079224, 58.53079224, 58.53079224],\n", - " [58.68289948, 58.68289948, 58.68289948, 58.68289948],\n", - " [58.83491898, 58.83491898, 58.83491898, 58.83491898],\n", - " ...,\n", - " [58.83491898, 58.83491898, 58.83491898, 58.83491898],\n", - " [58.68290329, 58.68290329, 58.68290329, 58.68290329],\n", - " [58.53079224, 58.53079224, 58.53079224, 58.53079224]]],\n", - " mask=False,\n", - " fill_value=1e+20)" + " [[58.29776072, 58.44964485, 58.55466327, 58.40259426],\n", + " [58.44964543, 58.6014251 , 58.7066318 , 58.55466385],\n", + " [58.60142568, 58.75309876, 58.85849715, 58.70663238],\n", + " ...,\n", + " [58.75309587, 58.60142279, 58.70662948, 58.85849425],\n", + " [58.6014251 , 58.44964543, 58.55466385, 58.7066318 ],\n", + " [58.44964485, 58.29776072, 58.40259426, 58.55466327]],\n", + "\n", + " [[58.40259366, 58.55466267, 58.65885172, 58.50660166],\n", + " [58.55466325, 58.7066312 , 58.81100467, 58.6588523 ],\n", + " [58.70663178, 58.85849655, 58.96305787, 58.81100525],\n", + " ...,\n", + " [58.85849365, 58.70662888, 58.81100235, 58.96305497],\n", + " [58.7066312 , 58.55466325, 58.6588523 , 58.81100467],\n", + " [58.55466267, 58.40259366, 58.50660166, 58.65885172]]])" ] }, "execution_count": 5, @@ -149,58 +146,55 @@ { "data": { "text/plain": [ - "masked_array(\n", - " data=[[[-22.32515679, -22.03737297, -22.03737297, -22.32515679],\n", - " [-22.16056404, -21.87278023, -21.87278023, -22.16056404],\n", - " [-21.99569092, -21.7079071 , -21.7079071 , -21.99569092],\n", - " ...,\n", - " [ 41.70790329, 41.99568711, 41.99568711, 41.70790329],\n", - " [ 41.8727745 , 42.16055832, 42.16055832, 41.8727745 ],\n", - " [ 42.03736725, 42.32515106, 42.32515106, 42.03736725]],\n", - "\n", - " [[-22.42207108, -22.13428726, -22.13428726, -22.42207108],\n", - " [-22.25707779, -21.96929397, -21.96929397, -22.25707779],\n", - " [-22.0917965 , -21.80401268, -21.80401268, -22.0917965 ],\n", - " ...,\n", - " [ 41.80400696, 42.09179077, 42.09179077, 41.80400696],\n", - " [ 41.96928253, 42.25706635, 42.25706635, 41.96928253],\n", - " [ 42.13427963, 42.42206345, 42.42206345, 42.13427963]],\n", + "array([[[-22.21497021, -22.05071303, -22.14733617, -22.31199395],\n", + " [-22.0507124 , -21.88618013, -21.9824008 , -22.14733554],\n", + " [-21.8861795 , -21.72137239, -21.81718872, -21.98240017],\n", + " ...,\n", + " [ 41.72137553, 41.88618264, 41.98240332, 41.81719187],\n", + " [ 41.88618013, 42.0507124 , 42.14733554, 41.9824008 ],\n", + " [ 42.05071303, 42.21497021, 42.31199395, 42.14733617]],\n", "\n", - " [[-22.51915894, -22.23137512, -22.23137512, -22.51915894],\n", - " [-22.35376511, -22.06598129, -22.06598129, -22.35376511],\n", - " [-22.18808136, -21.90029754, -21.90029754, -22.18808136],\n", - " ...,\n", - " [ 41.90029373, 42.18807755, 42.18807755, 41.90029373],\n", - " [ 42.06598129, 42.35376511, 42.35376511, 42.06598129],\n", - " [ 42.23137131, 42.51915512, 42.51915512, 42.23137131]],\n", + " [[-22.31199432, -22.14733654, -22.24413665, -22.4091946 ],\n", + " [-22.14733591, -21.98240117, -22.07879923, -22.24413602],\n", + " [-21.98240054, -21.81718908, -21.91318321, -22.0787986 ],\n", + " ...,\n", + " [ 41.81719223, 41.98240369, 42.07880176, 41.91318637],\n", + " [ 41.98240117, 42.14733591, 42.24413602, 42.07879923],\n", + " [ 42.14733654, 42.31199432, 42.4091946 , 42.24413665]],\n", "\n", + " [[-22.40919405, -22.2441361 , -22.34111548, -22.50657316],\n", + " [-22.24413547, -22.07879868, -22.17537644, -22.34111485],\n", + " [-22.07879805, -21.91318266, -22.00935688, -22.1753758 ],\n", " ...,\n", + " [ 41.91318582, 42.07880121, 42.17537897, 42.00936005],\n", + " [ 42.07879868, 42.24413547, 42.34111485, 42.17537644],\n", + " [ 42.2441361 , 42.40919405, 42.50657316, 42.34111548]],\n", "\n", - " [[-67.72155914, -67.43377533, -67.43377533, -67.72155914],\n", - " [-67.54095612, -67.2531723 , -67.2531723 , -67.54095612],\n", - " [-67.3592392 , -67.07145538, -67.07145538, -67.3592392 ],\n", - " ...,\n", - " [ 87.07144775, 87.35923157, 87.35923157, 87.07144775],\n", - " [ 87.25316467, 87.54094849, 87.54094849, 87.25316467],\n", - " [ 87.4337677 , 87.72155152, 87.72155152, 87.4337677 ]],\n", + " ...,\n", "\n", - " [[-68.04577027, -67.75798645, -67.75798645, -68.04577027],\n", - " [-67.86636505, -67.57858124, -67.57858124, -67.86636505],\n", - " [-67.68583069, -67.39804687, -67.39804687, -67.68583069],\n", - " ...,\n", - " [ 87.39804687, 87.68583069, 87.68583069, 87.39804687],\n", - " [ 87.57856598, 87.86634979, 87.86634979, 87.57856598],\n", - " [ 87.75797882, 88.04576264, 88.04576264, 87.75797882]],\n", + " [[-67.50645709, -67.32583243, -67.64966627, -67.82912696],\n", + " [-67.32583174, -67.14410165, -67.46910621, -67.64966558],\n", + " [-67.14410095, -66.96124932, -67.28743133, -67.46910552],\n", + " ...,\n", + " [ 86.96125282, 87.14410443, 87.46910897, 87.28743481],\n", + " [ 87.14410165, 87.32583174, 87.64966558, 87.46910621],\n", + " [ 87.32583243, 87.50645709, 87.82912696, 87.64966627]],\n", "\n", - " [[-68.37192688, -68.08414306, -68.08414306, -68.37192688],\n", - " [-68.19371185, -67.90592804, -67.90592804, -68.19371185],\n", - " [-68.01440582, -67.72662201, -67.72662201, -68.01440582],\n", - " ...,\n", - " [ 87.72661438, 88.0143982 , 88.0143982 , 87.72661438],\n", - " [ 87.90592804, 88.19371185, 88.19371185, 87.90592804],\n", - " [ 88.08414306, 88.37192688, 88.37192688, 88.08414306]]],\n", - " mask=False,\n", - " fill_value=1e+20)" + " [[-67.82912819, -67.64966751, -67.97544812, -68.1537229 ],\n", + " [-67.64966682, -67.46910745, -67.79608108, -67.97544744],\n", + " [-67.46910676, -67.28743258, -67.6156063 , -67.79608039],\n", + " ...,\n", + " [ 87.28743606, 87.46911022, 87.79608382, 87.61560976],\n", + " [ 87.46910745, 87.64966682, 87.97544744, 87.79608108],\n", + " [ 87.64966751, 87.82912819, 88.1537229 , 87.97544812]],\n", + "\n", + " [[-68.15372103, -67.97544625, -68.30317799, -68.48024479],\n", + " [-67.97544557, -67.7960792 , -68.12502637, -68.30317732],\n", + " [-67.79607851, -67.61560442, -67.94577447, -68.12502569],\n", + " ...,\n", + " [ 87.61560787, 87.79608195, 88.1250291 , 87.9457779 ],\n", + " [ 87.7960792 , 87.97544557, 88.30317732, 88.12502637],\n", + " [ 87.97544625, 88.15372103, 88.48024479, 88.30317799]]])" ] }, "execution_count": 6, diff --git a/tutorials/5.Others/5.3.To_Shapefile.ipynb b/tutorials/5.Others/5.3.To_Shapefile.ipynb new file mode 100644 index 0000000..07df9d2 --- /dev/null +++ b/tutorials/5.Others/5.3.To_Shapefile.ipynb @@ -0,0 +1,210 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# How to save objects as shapefiles" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [], + "source": [ + "from nes import *" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## 1. From grids that have already been created" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [], + "source": [ + "regular_grid_path = '/gpfs/scratch/bsc32/bsc32538/mr_multiplyby/original_file/MONARCH_d01_2008123100.nc'\n", + "regular_grid = open_netcdf(path=regular_grid_path, info=True)\n", + "regular_grid.to_shapefile('regular.shp')" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [], + "source": [ + "rotated_grid_path = '/gpfs/scratch/bsc32/bsc32538/original_files/CAMS_MONARCH_d01_2022070412.nc'\n", + "rotated_grid = open_netcdf(path=rotated_grid_path, info=True)\n", + "rotated_grid.to_shapefile('rotated.shp')" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [], + "source": [ + "lcc_grid_path = '/esarchive/exp/wrf-hermes-cmaq/b075/eu/hourly/sconco3/sconco3_2021010100.nc'\n", + "lcc_grid = open_netcdf(path=lcc_grid_path, info=True)\n", + "lcc_grid.to_shapefile('lcc.shp')" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [], + "source": [ + "mercator_grid_path = '/gpfs/scratch/bsc32/bsc32538/a4mg/HERMESv3/auxiliar_files/d01/temporal_coords.nc'\n", + "mercator_grid = open_netcdf(path=mercator_grid_path, info=True)\n", + "mercator_grid.to_shapefile('mercator.shp')" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## 2. From grids created from scratch" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": {}, + "outputs": [], + "source": [ + "lat_orig = 40\n", + "lon_orig = 0\n", + "inc_lat = 0.5\n", + "inc_lon = 0.5\n", + "n_lat = 100\n", + "n_lon = 100\n", + "regular_grid = create_nes(comm=None, info=False, projection='regular',\n", + " lat_orig=lat_orig, lon_orig=lon_orig, inc_lat=inc_lat, inc_lon=inc_lon, \n", + " n_lat=n_lat, n_lon=n_lon)" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": {}, + "outputs": [], + "source": [ + "regular_grid.to_shapefile('create_regular.shp')" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": {}, + "outputs": [], + "source": [ + "centre_lat = 51\n", + "centre_lon = 10\n", + "west_boundary = -35\n", + "south_boundary = -27\n", + "inc_rlat = 0.2\n", + "inc_rlon = 0.2\n", + "rotated_grid = create_nes(comm=None, info=False, projection='rotated',\n", + " centre_lat=centre_lat, centre_lon=centre_lon,\n", + " west_boundary=west_boundary, south_boundary=south_boundary,\n", + " inc_rlat=inc_rlat, inc_rlon=inc_rlon)" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": {}, + "outputs": [], + "source": [ + "rotated_grid.to_shapefile('create_rotated.shp')" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "metadata": {}, + "outputs": [], + "source": [ + "lat_1 = 37\n", + "lat_2 = 43\n", + "lon_0 = -3\n", + "lat_0 = 40\n", + "nx = 397\n", + "ny = 397\n", + "inc_x = 4000\n", + "inc_y = 4000\n", + "x_0 = -807847.688\n", + "y_0 = -797137.125\n", + "lcc_grid = create_nes(comm=None, info=False, projection='lcc',\n", + " lat_1=lat_1, lat_2=lat_2, lon_0=lon_0, lat_0=lat_0, \n", + " nx=nx, ny=ny, inc_x=inc_x, inc_y=inc_y, x_0=x_0, y_0=y_0)" + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "metadata": {}, + "outputs": [], + "source": [ + "lcc_grid.to_shapefile('create_lcc.shp')" + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "metadata": {}, + "outputs": [], + "source": [ + "lat_ts = -1.5\n", + "lon_0 = -18.0\n", + "nx = 210\n", + "ny = 236\n", + "inc_x = 50000\n", + "inc_y = 50000\n", + "x_0 = -126017.5\n", + "y_0 = -5407460.0\n", + "mercator_grid = create_nes(comm=None, info=False, projection='mercator',\n", + " lat_ts=lat_ts, lon_0=lon_0, nx=nx, ny=ny, \n", + " inc_x=inc_x, inc_y=inc_y, x_0=x_0, y_0=y_0)" + ] + }, + { + "cell_type": "code", + "execution_count": 13, + "metadata": {}, + "outputs": [], + "source": [ + "mercator_grid.to_shapefile('create_mercator.shp')" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.7.4" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +} diff --git a/tutorials/Jupyter_bash_nord3v2.cmd b/tutorials/Jupyter_bash_nord3v2.cmd index 563043a..8138ab9 100644 --- a/tutorials/Jupyter_bash_nord3v2.cmd +++ b/tutorials/Jupyter_bash_nord3v2.cmd @@ -32,6 +32,8 @@ module load cfunits/1.8-foss-2019b-Python-3.7.4 module load filelock/3.7.1-foss-2019b-Python-3.7.4 module load pyproj/2.5.0-foss-2019b-Python-3.7.4 module load eccodes-python/0.9.5-foss-2019b-Python-3.7.4 +module load geopandas/0.8.1-foss-2019b-Python-3.7.4 +module load Shapely/1.7.1-foss-2019b-Python-3.7.4 #module load NES/0.9.0-foss-2019b-Python-3.7.4 export PYTHONPATH=/esarchive/scratch/avilanova/software/NES:${PYTHONPATH} -- GitLab