diff --git a/nes/nc_projections/default_nes.py b/nes/nc_projections/default_nes.py index aca209c9ff1625af3254d1290a410ef261feb7fc..493a1174c2d116fab11ed8b24e4ef03dacd5d70d 100644 --- a/nes/nc_projections/default_nes.py +++ b/nes/nc_projections/default_nes.py @@ -4,6 +4,7 @@ import sys import os import warnings import numpy as np +import pandas as pd import datetime from xarray import open_dataset from netCDF4 import Dataset, num2date, date2num @@ -136,6 +137,7 @@ class Nes(object): self.info = info self.is_xarray = xarray self.__ini_path = path + self.shapefile = None # Selecting info self.hours_start = avoid_first_hours @@ -2063,14 +2065,9 @@ class Nes(object): return None - def to_shapefile(self, path): + def create_shapefile(self): """ - Write output file with shapefile format. - - Parameters - ---------- - path : str - Path to the output file. + Create spatial geodataframe (shapefile). """ if self._lat_bnds is None or self._lon_bnds is None: @@ -2107,9 +2104,37 @@ class Nes(object): (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) + fids = np.arange(len(self._lat['data']) * len(self._lon['data'])) + fids = fids.reshape((len(self._lat['data']), len(self._lon['data']))) + fids = fids[self.read_axis_limits['y_min']:self.read_axis_limits['y_max'], + self.read_axis_limits['x_min']:self.read_axis_limits['x_max']] + gdf = gpd.GeoDataFrame(index=pd.Index(name='FID', data=fids.ravel()), geometry=geometry, crs="EPSG:4326") + self.shapefile = gdf + return gdf + + def write_shapefile(self, path): + """Save spatial geodataframe (shapefile). + + Parameters + ---------- + path : str + Path to the output file. + """ + + if self.shapefile is None: + raise ValueError('Shapefile was not created.') + + if self.size == 1: + # In serial, avoid gather + self.shapefile.to_file(path) + else: + # In parallel + data = self.comm.gather(self.shapefile, root=0) + if self.master: + data = pd.concat(data) + data.to_file(path) + return None def __gather_data_py_object(self): diff --git a/nes/nc_projections/lcc_nes.py b/nes/nc_projections/lcc_nes.py index f6605a6f08c60acb3a0508a25a2945452909d8ca..e8048044774c87488130d87ade90335717c5663d 100644 --- a/nes/nc_projections/lcc_nes.py +++ b/nes/nc_projections/lcc_nes.py @@ -1,6 +1,7 @@ #!/usr/bin/env python import numpy as np +import pandas as pd from cfunits import Units from pyproj import Proj from copy import deepcopy @@ -429,14 +430,9 @@ class LCCNes(Nes): raise NotImplementedError("Grib2 format cannot be written in a Lambert Conformal Conic projection.") - def to_shapefile(self, path): + def create_shapefile(self): """ - Write output file with shapefile format. - - Parameters - ---------- - path : str - Path to the output file. + Create spatial geodataframe (shapefile). """ # Get latitude and longitude cell boundaries @@ -455,7 +451,11 @@ class LCCNes(Nes): (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 + fids = np.arange(self._lat['data'].shape[0] * self._lat['data'].shape[1]) + fids = fids.reshape((self._lat['data'].shape[0],self._lat['data'].shape[1])) + fids = fids[self.read_axis_limits['y_min']:self.read_axis_limits['y_max'], + self.read_axis_limits['x_min']:self.read_axis_limits['x_max']] + gdf = gpd.GeoDataFrame(index=pd.Index(name='FID', data=fids.ravel()), geometry=geometry, crs="EPSG:4326") + self.shapefile = gdf + + return gdf diff --git a/nes/nc_projections/mercator_nes.py b/nes/nc_projections/mercator_nes.py index e5272aa08d041cdcc5b18fd58dcb30c4942dbd77..4ca7d4e3cb91e18d0acc67a7ab784700377c06bc 100644 --- a/nes/nc_projections/mercator_nes.py +++ b/nes/nc_projections/mercator_nes.py @@ -1,6 +1,7 @@ #!/usr/bin/env python import numpy as np +import pandas as pd from cfunits import Units from pyproj import Proj from copy import deepcopy @@ -406,14 +407,9 @@ class MercatorNes(Nes): raise NotImplementedError("Grib2 format cannot be written in a Mercator projection.") - def to_shapefile(self, path): + def create_shapefile(self): """ - Write output file with shapefile format. - - Parameters - ---------- - path : str - Path to the output file. + Create spatial geodataframe (shapefile). """ # Get latitude and longitude cell boundaries @@ -432,7 +428,11 @@ class MercatorNes(Nes): (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 + fids = np.arange(self._lat['data'].shape[0] * self._lat['data'].shape[1]) + fids = fids.reshape((self._lat['data'].shape[0],self._lat['data'].shape[1])) + fids = fids[self.read_axis_limits['y_min']:self.read_axis_limits['y_max'], + self.read_axis_limits['x_min']:self.read_axis_limits['x_max']] + gdf = gpd.GeoDataFrame(index=pd.Index(name='FID', data=fids.ravel()), geometry=geometry, crs="EPSG:4326") + self.shapefile = gdf + + return gdf diff --git a/nes/nc_projections/points_nes.py b/nes/nc_projections/points_nes.py index a647cc035149cdf7a4fbaf736ea215da1ece81ed..8d5ccff15b5b5ff5ce645ea58e8dbb90eac5cadc 100644 --- a/nes/nc_projections/points_nes.py +++ b/nes/nc_projections/points_nes.py @@ -640,14 +640,9 @@ class PointsNes(Nes): raise NotImplementedError("Grib2 format cannot be written with point data.") - def to_shapefile(self, path): + def create_shapefile(self): """ - Write output file with shapefile format. - - Parameters - ---------- - path : str - Path to the output file. + Create spatial geodataframe (shapefile). """ 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 c5b2045b924c47c4d1029f309e72efd0a1f6a819..54f9e44d2a4ef4068e0303878a4bf8a27a07d0f9 100644 --- a/nes/nc_projections/rotated_nes.py +++ b/nes/nc_projections/rotated_nes.py @@ -1,6 +1,7 @@ #!/usr/bin/env python import numpy as np +import pandas as pd import math from cfunits import Units from copy import deepcopy @@ -456,14 +457,9 @@ class RotatedNes(Nes): raise NotImplementedError("Grib2 format cannot be written in a Rotated pole projection.") - def to_shapefile(self, path): + def create_shapefile(self): """ - Write output file with shapefile format. - - Parameters - ---------- - path : str - Path to the output file. + Create spatial geodataframe (shapefile). """ if self._lat_bnds is None or self._lon_bnds is None: @@ -481,7 +477,11 @@ class RotatedNes(Nes): (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 + fids = np.arange(self._lat['data'].shape[0] * self._lat['data'].shape[1]) + fids = fids.reshape((self._lat['data'].shape[0],self._lat['data'].shape[1])) + fids = fids[self.read_axis_limits['y_min']:self.read_axis_limits['y_max'], + self.read_axis_limits['x_min']:self.read_axis_limits['x_max']] + gdf = gpd.GeoDataFrame(index=pd.Index(name='FID', data=fids.ravel()), geometry=geometry, crs="EPSG:4326") + self.shapefile = gdf + + return gdf diff --git a/tutorials/5.Others/5.3.Shapefiles.ipynb b/tutorials/5.Others/5.3.Shapefiles.ipynb new file mode 100644 index 0000000000000000000000000000000000000000..6fe590da742d02b786b853f037cf007edfe3e2fb --- /dev/null +++ b/tutorials/5.Others/5.3.Shapefiles.ipynb @@ -0,0 +1,504 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# How to create and save shapefiles" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [], + "source": [ + "from nes import *\n", + "import xarray as xr" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## 1. From grids that have already been created" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [], + "source": [ + "grid_path = '/gpfs/scratch/bsc32/bsc32538/original_files/CAMS_MONARCH_d01_2022070412.nc'\n", + "grid = open_netcdf(path=grid_path, info=True)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Create shapefile" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "
\n", + " | geometry | \n", + "
---|---|
FID | \n", + "\n", + " |
0 | \n", + "POLYGON ((-22.16553 16.27358, -22.04224 16.335... | \n", + "
1 | \n", + "POLYGON ((-22.04224 16.33554, -21.91879 16.397... | \n", + "
2 | \n", + "POLYGON ((-21.91879 16.39727, -21.79523 16.458... | \n", + "
3 | \n", + "POLYGON ((-21.79523 16.45880, -21.67145 16.520... | \n", + "
4 | \n", + "POLYGON ((-21.67145 16.52011, -21.54755 16.581... | \n", + "
... | \n", + "... | \n", + "
168582 | \n", + "POLYGON ((87.45258 59.04235, 87.58887 58.92854... | \n", + "
168583 | \n", + "POLYGON ((87.58887 58.92854, 87.72452 58.81466... | \n", + "
168584 | \n", + "POLYGON ((87.72452 58.81466, 87.85956 58.70073... | \n", + "
168585 | \n", + "POLYGON ((87.85956 58.70073, 87.99396 58.58673... | \n", + "
168586 | \n", + "POLYGON ((87.99396 58.58673, 88.12778 58.47268... | \n", + "
168587 rows × 1 columns
\n", + "\n", + " | geometry | \n", + "O3 | \n", + "
---|---|---|
FID | \n", + "\n", + " | \n", + " |
0 | \n", + "POLYGON ((-22.16553 16.27358, -22.04224 16.335... | \n", + "0.037919 | \n", + "
1 | \n", + "POLYGON ((-22.04224 16.33554, -21.91879 16.397... | \n", + "0.038064 | \n", + "
2 | \n", + "POLYGON ((-21.91879 16.39727, -21.79523 16.458... | \n", + "0.038086 | \n", + "
3 | \n", + "POLYGON ((-21.79523 16.45880, -21.67145 16.520... | \n", + "0.038092 | \n", + "
4 | \n", + "POLYGON ((-21.67145 16.52011, -21.54755 16.581... | \n", + "0.038098 | \n", + "
... | \n", + "... | \n", + "... | \n", + "
168582 | \n", + "POLYGON ((87.45258 59.04235, 87.58887 58.92854... | \n", + "0.026542 | \n", + "
168583 | \n", + "POLYGON ((87.58887 58.92854, 87.72452 58.81466... | \n", + "0.026310 | \n", + "
168584 | \n", + "POLYGON ((87.72452 58.81466, 87.85956 58.70073... | \n", + "0.027037 | \n", + "
168585 | \n", + "POLYGON ((87.85956 58.70073, 87.99396 58.58673... | \n", + "0.027271 | \n", + "
168586 | \n", + "POLYGON ((87.99396 58.58673, 88.12778 58.47268... | \n", + "0.026327 | \n", + "
168587 rows × 2 columns
\n", + "\n", + " | geometry | \n", + "
---|---|
FID | \n", + "\n", + " |
0 | \n", + "POLYGON ((-22.21497 16.22040, -22.05071 16.303... | \n", + "
1 | \n", + "POLYGON ((-22.05071 16.30307, -21.88618 16.385... | \n", + "
2 | \n", + "POLYGON ((-21.88618 16.38536, -21.72137 16.467... | \n", + "
3 | \n", + "POLYGON ((-21.72137 16.46727, -21.55629 16.548... | \n", + "
4 | \n", + "POLYGON ((-21.55629 16.54881, -21.39094 16.629... | \n", + "
... | \n", + "... | \n", + "
95116 | \n", + "POLYGON ((87.25127 59.16191, 87.43401 59.01025... | \n", + "
95117 | \n", + "POLYGON ((87.43401 59.01025, 87.61561 58.85850... | \n", + "
95118 | \n", + "POLYGON ((87.61561 58.85850, 87.79608 58.70663... | \n", + "
95119 | \n", + "POLYGON ((87.79608 58.70663, 87.97545 58.55466... | \n", + "
95120 | \n", + "POLYGON ((87.97545 58.55466, 88.15372 58.40259... | \n", + "
95121 rows × 1 columns
\n", + "