From 8dc224b8eed9eb36925727056b981240b1f6a365 Mon Sep 17 00:00:00 2001 From: Alba Vilanova Date: Tue, 14 Feb 2023 07:18:35 +0100 Subject: [PATCH 1/2] Add global projection --- nes/create_nes.py | 10 +++++++++- nes/nc_projections/latlon_nes.py | 31 +++++++++++++++++++++--------- nes/nc_projections/lcc_nes.py | 6 ++---- nes/nc_projections/mercator_nes.py | 6 ++---- nes/nc_projections/points_nes.py | 8 +++----- nes/nc_projections/rotated_nes.py | 25 ++++++++++++------------ 6 files changed, 50 insertions(+), 36 deletions(-) diff --git a/nes/create_nes.py b/nes/create_nes.py index 8c5980b..33f0bf6 100644 --- a/nes/create_nes.py +++ b/nes/create_nes.py @@ -61,6 +61,8 @@ def create_nes(comm=None, info=False, projection=None, parallel_method='Y', bala required_vars = ['lat', 'lon'] elif projection == 'regular': required_vars = ['lat_orig', 'lon_orig', 'inc_lat', 'inc_lon', 'n_lat', 'n_lon'] + elif projection == 'global': + required_vars = ['inc_lat', 'inc_lon'] elif projection == 'rotated': required_vars = ['centre_lat', 'centre_lon', 'west_boundary', 'south_boundary', 'inc_rlat', 'inc_rlon'] elif projection == 'lcc': @@ -76,6 +78,12 @@ def create_nes(comm=None, info=False, projection=None, parallel_method='Y', bala msg += 'For a {} projection, it is necessary to define {}'.format(projection, required_vars) raise ValueError(msg) + for var in kwargs_list: + if var not in required_vars: + msg = 'Variable {0} has been defined. '.format(var) + msg += 'For a {} projection, you can only define {}'.format(projection, required_vars) + raise ValueError(msg) + if projection is None: if parallel_method == 'Y': warnings.warn("Parallel method cannot be 'Y' to create points NES. Setting it to 'X'") @@ -86,7 +94,7 @@ def create_nes(comm=None, info=False, projection=None, parallel_method='Y', bala avoid_first_hours=avoid_first_hours, avoid_last_hours=avoid_last_hours, first_level=first_level, last_level=last_level, balanced=balanced, create_nes=True, strlen=strlen, times=times, **kwargs) - elif projection == 'regular': + elif projection in ['regular', 'global']: nessy = LatLonNes(comm=comm, dataset=None, xarray=False, info=info, parallel_method=parallel_method, avoid_first_hours=avoid_first_hours, avoid_last_hours=avoid_last_hours, first_level=first_level, last_level=last_level, balanced=balanced, diff --git a/nes/nc_projections/latlon_nes.py b/nes/nc_projections/latlon_nes.py index 31acd58..32756f1 100644 --- a/nes/nc_projections/latlon_nes.py +++ b/nes/nc_projections/latlon_nes.py @@ -129,19 +129,32 @@ class LatLonNes(Nes): Dictionary with data of centre longitudes in 1D """ + inc_lat = kwargs['inc_lat'] + inc_lon = kwargs['inc_lon'] + + # Global domain + if len(kwargs) == 2: + lat_orig = -90 + lon_orig = -180 + n_lat = int(180 // inc_lat) + n_lon = int(360 // inc_lon) + + # Other domains + else: + lat_orig = kwargs['lat_orig'] + lon_orig = kwargs['lon_orig'] + n_lat = kwargs['n_lon'] + n_lon = kwargs['n_lon'] + # Calculate centre latitudes - lat_c_orig = kwargs['lat_orig'] + (kwargs['inc_lat'] / 2) - centre_lat_data = np.linspace( - lat_c_orig, lat_c_orig + (kwargs['inc_lat'] * (kwargs['n_lat'] - 1)), kwargs['n_lat']) - centre_lat = {'data': centre_lat_data} + lat_c_orig = lat_orig + (inc_lat / 2) + centre_lat = np.linspace(lat_c_orig, lat_c_orig + (inc_lat * (n_lat - 1)), n_lat) # Calculate centre longitudes - lon_c_orig = kwargs['lon_orig'] + (kwargs['inc_lon'] / 2) - centre_lon_data = np.linspace( - lon_c_orig, lon_c_orig + (kwargs['inc_lon'] * (kwargs['n_lon'] - 1)), kwargs['n_lon']) - centre_lon = {'data': centre_lon_data} + lon_c_orig = lon_orig + (inc_lon / 2) + centre_lon = np.linspace(lon_c_orig, lon_c_orig + (inc_lon * (n_lon - 1)), n_lon) - return centre_lat, centre_lon + return {'data': centre_lat}, {'data': centre_lon} def create_providentia_exp_centre_coordinates(self): """ diff --git a/nes/nc_projections/lcc_nes.py b/nes/nc_projections/lcc_nes.py index 6f107ea..ad4f1ab 100644 --- a/nes/nc_projections/lcc_nes.py +++ b/nes/nc_projections/lcc_nes.py @@ -272,11 +272,9 @@ class LCCNes(Nes): ) # Calculate centre latitudes and longitudes (UTM to LCC) - centre_lon_data, centre_lat_data = self.projection(x, y, inverse=True) - centre_lat = {'data': centre_lat_data} - centre_lon = {'data': centre_lon_data} + centre_lon, centre_lat = self.projection(x, y, inverse=True) - return centre_lat, centre_lon + return {'data': centre_lat}, {'data': centre_lon} def create_providentia_exp_centre_coordinates(self): """ diff --git a/nes/nc_projections/mercator_nes.py b/nes/nc_projections/mercator_nes.py index c539bcd..433a528 100644 --- a/nes/nc_projections/mercator_nes.py +++ b/nes/nc_projections/mercator_nes.py @@ -263,11 +263,9 @@ class MercatorNes(Nes): ) # Calculate centre latitudes and longitudes (UTM to Mercator) - centre_lon_data, centre_lat_data = self.projection(x, y, inverse=True) - centre_lat = {'data': centre_lat_data} - centre_lon = {'data': centre_lon_data} + centre_lon, centre_lat = self.projection(x, y, inverse=True) - return centre_lat, centre_lon + return {'data': centre_lat}, {'data': centre_lon} def create_providentia_exp_centre_coordinates(self): """ diff --git a/nes/nc_projections/points_nes.py b/nes/nc_projections/points_nes.py index 1ca8c86..55d9a17 100644 --- a/nes/nc_projections/points_nes.py +++ b/nes/nc_projections/points_nes.py @@ -552,14 +552,12 @@ class PointsNes(Nes): """ # Calculate centre latitudes - centre_lat_data = kwargs['lat'] - centre_lat = {'data': centre_lat_data} + centre_lat = kwargs['lat'] # Calculate centre longitudes - centre_lon_data = kwargs['lon'] - centre_lon = {'data': centre_lon_data} + centre_lon = kwargs['lon'] - return centre_lat, centre_lon + return {'data': centre_lat}, {'data': centre_lon} def _create_metadata(self, netcdf): """ diff --git a/nes/nc_projections/rotated_nes.py b/nes/nc_projections/rotated_nes.py index 42b6ebe..36b1ddd 100644 --- a/nes/nc_projections/rotated_nes.py +++ b/nes/nc_projections/rotated_nes.py @@ -240,15 +240,16 @@ class RotatedNes(Nes): """ # Calculate rotated latitudes - self.n_lat = int((abs(kwargs['south_boundary']) / kwargs['inc_rlat']) * 2 + 1) - self.rotated_lat = np.linspace(kwargs['south_boundary'], kwargs['south_boundary'] + - (kwargs['inc_rlat'] * (self.n_lat - 1)), self.n_lat) + n_lat = int((abs(kwargs['south_boundary']) / kwargs['inc_rlat']) * 2 + 1) + rlat = np.linspace(kwargs['south_boundary'], kwargs['south_boundary'] + + (kwargs['inc_rlat'] * (n_lat - 1)), n_lat) + # Calculate rotated longitudes - self.n_lon = int((abs(kwargs['west_boundary']) / kwargs['inc_rlon']) * 2 + 1) - self.rotated_lon = np.linspace(kwargs['west_boundary'], kwargs['west_boundary'] + - (kwargs['inc_rlon'] * (self.n_lon - 1)), self.n_lon) + n_lon = int((abs(kwargs['west_boundary']) / kwargs['inc_rlon']) * 2 + 1) + rlon = np.linspace(kwargs['west_boundary'], kwargs['west_boundary'] + + (kwargs['inc_rlon'] * (n_lon - 1)), n_lon) - return {'data': self.rotated_lat}, {'data': self.rotated_lon} + return {'data': rlat}, {'data': rlon} def rotated2latlon(self, lon_deg, lat_deg, lon_min=-180, **kwargs): """ @@ -324,13 +325,11 @@ class RotatedNes(Nes): self._rlat, self._rlon = self._create_rotated_coordinates(**kwargs) # Calculate centre latitudes and longitudes (1D to 2D) - centre_lon_data, centre_lat_data = self.rotated2latlon(np.array([self.rotated_lon] * len(self.rotated_lat)), - np.array([self.rotated_lat] * len(self.rotated_lon)).T, - **kwargs) - centre_lon = {'data': centre_lon_data} - centre_lat = {'data': centre_lat_data} + centre_lon, centre_lat = self.rotated2latlon(np.array([self._rlon['data']] * len(self._rlat['data'])), + np.array([self._rlat['data']] * len(self._rlon['data'])).T, + **kwargs) - return centre_lat, centre_lon + return {'data': centre_lat}, {'data': centre_lon} def create_providentia_exp_centre_coordinates(self): """ -- GitLab From c08a4e5676f83b37152ef34f8ad4d3d5b6ff06d2 Mon Sep 17 00:00:00 2001 From: ctena Date: Tue, 14 Feb 2023 13:20:04 +0100 Subject: [PATCH 2/2] Create global grid: Added test & tutorial --- tests/1.2-test_create_projection.py | 196 ++++++++ tutorials/2.Creation/2.8.Create_Global.ipynb | 475 +++++++++++++++++++ 2 files changed, 671 insertions(+) create mode 100644 tests/1.2-test_create_projection.py create mode 100644 tutorials/2.Creation/2.8.Create_Global.ipynb diff --git a/tests/1.2-test_create_projection.py b/tests/1.2-test_create_projection.py new file mode 100644 index 0000000..5b26ce0 --- /dev/null +++ b/tests/1.2-test_create_projection.py @@ -0,0 +1,196 @@ +#!/usr/bin/env python + +import sys +from mpi4py import MPI +import pandas as pd +import timeit +from nes import * + +comm = MPI.COMM_WORLD +rank = comm.Get_rank() +size = comm.Get_size() + +# Original path: /esarchive/exp/monarch/a4dd/original_files/000/2022111512/MONARCH_d01_2022111512.nc +# Rotated grid from MONARCH +original_path = '/gpfs/projects/bsc32/models/NES_tutorial_data/MONARCH_d01_2022111512.nc' + +parallel_method = 'Y' + +result_path = "Times_test_1.2_create_projection_{0}_{1:03d}.csv".format(parallel_method, size) +result = pd.DataFrame(index=['read', 'calculate', 'write'], + columns=['1.2.1.Regular', '1.2.2.Rotated', '1.2.3.LCC', '1.2.4.Mercator', '1.2.5.Global']) + +# ===== PATH TO MASK ===== # +shapefile_path = '/gpfs/projects/bsc32/models/NES_tutorial_data/timezones_2021c/timezones_2021c.shp' + +# ====================================================================================================================== +# ============================================= REGULAR ======================================================== +# ====================================================================================================================== + +test_name = '1.2.1.Regular' +if rank == 0: + print(test_name) +# READ +st_time = timeit.default_timer() + +lat_orig = 41.1 +lon_orig = 1.8 +inc_lat = 0.01 +inc_lon = 0.01 +n_lat = 100 +n_lon = 100 +nessy = create_nes(projection='regular', parallel_method=parallel_method, + lat_orig=lat_orig, lon_orig=lon_orig, inc_lat=inc_lat, inc_lon=inc_lon, n_lat=n_lat, n_lon=n_lon) + +comm.Barrier() +result.loc['read', test_name] = timeit.default_timer() - st_time + +st_time = timeit.default_timer() +# Method can be centroid, nearest and intersection +nessy.to_netcdf(test_name.replace(' ', '_') + "_{0:03d}.nc".format(size)) +comm.Barrier() +result.loc['write', test_name] = timeit.default_timer() - st_time + +comm.Barrier() +if rank == 0: + print(result.loc[:, test_name]) +sys.stdout.flush() + +# ====================================================================================================================== +# ============================================= ROTATED ======================================================== +# ====================================================================================================================== + +test_name = '1.2.2.Rotated' +if rank == 0: + print(test_name) +# READ +st_time = timeit.default_timer() + +centre_lat = 51 +centre_lon = 10 +west_boundary = -35 +south_boundary = -27 +inc_rlat = 0.2 +inc_rlon = 0.2 +nessy = create_nes(projection='rotated', parallel_method=parallel_method, + centre_lat=centre_lat, centre_lon=centre_lon, + west_boundary=west_boundary, south_boundary=south_boundary, + inc_rlat=inc_rlat, inc_rlon=inc_rlon) + +comm.Barrier() +result.loc['read', test_name] = timeit.default_timer() - st_time + +st_time = timeit.default_timer() +# Method can be centroid, nearest and intersection +nessy.to_netcdf(test_name.replace(' ', '_') + "_{0:03d}.nc".format(size)) +comm.Barrier() +result.loc['write', test_name] = timeit.default_timer() - st_time + +comm.Barrier() +if rank == 0: + print(result.loc[:, test_name]) +sys.stdout.flush() + +# ====================================================================================================================== +# ============================================= LCC ======================================================== +# ====================================================================================================================== + +test_name = '1.2.3.LCC' +if rank == 0: + print(test_name) +# READ +st_time = timeit.default_timer() + +lat_1 = 37 +lat_2 = 43 +lon_0 = -3 +lat_0 = 40 +nx = 397 +ny = 397 +inc_x = 4000 +inc_y = 4000 +x_0 = -807847.688 +y_0 = -797137.125 +nessy = create_nes(projection='lcc', parallel_method=parallel_method, + lat_1=lat_1, lat_2=lat_2, lon_0=lon_0, lat_0=lat_0, + nx=nx, ny=ny, inc_x=inc_x, inc_y=inc_y, x_0=x_0, y_0=y_0) + +comm.Barrier() +result.loc['read', test_name] = timeit.default_timer() - st_time + +st_time = timeit.default_timer() +# Method can be centroid, nearest and intersection +nessy.to_netcdf(test_name.replace(' ', '_') + "_{0:03d}.nc".format(size)) +comm.Barrier() +result.loc['write', test_name] = timeit.default_timer() - st_time + +comm.Barrier() +if rank == 0: + print(result.loc[:, test_name]) +sys.stdout.flush() + +# ====================================================================================================================== +# ============================================= MERCATOR ======================================================== +# ====================================================================================================================== + +test_name = '1.2.4.Mercator' +if rank == 0: + print(test_name) +# READ +st_time = timeit.default_timer() + +lat_ts = -1.5 +lon_0 = -18.0 +nx = 210 +ny = 236 +inc_x = 50000 +inc_y = 50000 +x_0 = -126017.5 +y_0 = -5407460.0 +nessy = create_nes(projection='mercator', parallel_method=parallel_method, + lat_ts=lat_ts, lon_0=lon_0, nx=nx, ny=ny, inc_x=inc_x, inc_y=inc_y, x_0=x_0, y_0=y_0) + +comm.Barrier() +result.loc['read', test_name] = timeit.default_timer() - st_time + +st_time = timeit.default_timer() +# Method can be centroid, nearest and intersection +nessy.to_netcdf(test_name.replace(' ', '_') + "_{0:03d}.nc".format(size)) +comm.Barrier() +result.loc['write', test_name] = timeit.default_timer() - st_time + +comm.Barrier() +if rank == 0: + print(result.loc[:, test_name]) +sys.stdout.flush() + +# ====================================================================================================================== +# ============================================== GLOBAL ======================================================== +# ====================================================================================================================== + +test_name = '1.2.5.Global' +if rank == 0: + print(test_name) +# READ +st_time = timeit.default_timer() + +inc_lat = 0.1 +inc_lon = 0.1 +nessy = create_nes(projection='global', parallel_method=parallel_method, inc_lat=inc_lat, inc_lon=inc_lon) + +comm.Barrier() +result.loc['read', test_name] = timeit.default_timer() - st_time + +st_time = timeit.default_timer() +# Method can be centroid, nearest and intersection +nessy.to_netcdf(test_name.replace(' ', '_') + "_{0:03d}.nc".format(size)) +comm.Barrier() +result.loc['write', test_name] = timeit.default_timer() - st_time + +comm.Barrier() +if rank == 0: + print(result.loc[:, test_name]) +sys.stdout.flush() + +if rank == 0: + result.to_csv(result_path) diff --git a/tutorials/2.Creation/2.8.Create_Global.ipynb b/tutorials/2.Creation/2.8.Create_Global.ipynb new file mode 100644 index 0000000..5d3bf9a --- /dev/null +++ b/tutorials/2.Creation/2.8.Create_Global.ipynb @@ -0,0 +1,475 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# How to create global grids" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [], + "source": [ + "import xarray as xr\n", + "from nes import *" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [], + "source": [ + "inc_lat = 0.1\n", + "inc_lon = 0.1\n", + "grid = create_nes(projection='global', \n", + " inc_lat=inc_lat, inc_lon=inc_lon)" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Rank 000: Creating global_grid.nc\n", + "Rank 000: NetCDF ready to write\n", + "Rank 000: Dimensions done\n", + "Rank 000: Cell measures done\n" + ] + } + ], + "source": [ + "grid.to_netcdf('global_grid.nc', info=True)" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "
\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "
<xarray.Dataset>\n",
+       "Dimensions:  (time: 1, lev: 1, lat: 1799, lon: 3599)\n",
+       "Coordinates:\n",
+       "  * time     (time) datetime64[ns] 1996-12-31\n",
+       "  * lev      (lev) float64 0.0\n",
+       "  * lat      (lat) float64 -89.95 -89.85 -89.75 -89.65 ... 89.65 89.75 89.85\n",
+       "  * lon      (lon) float64 -179.9 -179.8 -179.8 -179.6 ... 179.7 179.8 179.9\n",
+       "Data variables:\n",
+       "    crs      |S1 b''\n",
+       "Attributes:\n",
+       "    Conventions:  CF-1.7
" + ], + "text/plain": [ + "\n", + "Dimensions: (time: 1, lev: 1, lat: 1799, lon: 3599)\n", + "Coordinates:\n", + " * time (time) datetime64[ns] 1996-12-31\n", + " * lev (lev) float64 0.0\n", + " * lat (lat) float64 -89.95 -89.85 -89.75 -89.65 ... 89.65 89.75 89.85\n", + " * lon (lon) float64 -179.9 -179.8 -179.8 -179.6 ... 179.7 179.8 179.9\n", + "Data variables:\n", + " crs |S1 ...\n", + "Attributes:\n", + " Conventions: CF-1.7" + ] + }, + "execution_count": 4, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "xr.open_dataset('global_grid.nc')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "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 +} -- GitLab