{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Creating a custom grid \n",
"\n",
"Emiproc provides many useful tools to work with gridded data.\n",
"The only thing you need is to tell emiproc how your grid is defined.\n",
"For this purpose emiproc implements an abstract class `Grid`. \n",
"\n",
"This tutorial aims to show how to create a custom grid for your own needs.\n",
"\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Understanding geopandas and shapely \n",
"\n",
"Geopandas will help us generate the geometry of the grid. \n",
"\n",
"The principle is that each grid cell is a polygon. \n",
"\n",
"Polygons are simply a list of 2d point coordinates that are connected by lines.\n",
"So only by defining the points of the polygon, we can create a grid cell.\n",
"\n",
"In geopandas you can create Series of polygons."
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"
Make this Notebook Trusted to load map: File -> Trust Notebook
"
],
"text/plain": [
""
]
},
"execution_count": 6,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"import geopandas as gpd\n",
"import shapely\n",
"\n",
"\n",
"polygon = shapely.geometry.Polygon([(0, 0), (1, 0), (1.5, 0.5), (1, 1), (0, 1)])\n",
"polygon2 = shapely.geometry.Polygon([(1.5, 0.5), (1.5, 1.5), (0.5, 1.5)])\n",
"\n",
"goeserie = gpd.GeoSeries([polygon, polygon2])\n",
"goeserie.explore()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Georeferencing the grid\n",
"\n",
"The grid can be georeferenced by defining the coordinate reference system (CRS) of the grid.\n",
"Most of the time we use the World Geodetic System 1984 (WGS84). https://en.wikipedia.org/wiki/World_Geodetic_System_1984 \n",
"\n",
"Thanks to geopandas, we can easily define the CRS of the grid.\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": 75,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"
Make this Notebook Trusted to load map: File -> Trust Notebook
"
],
"text/plain": [
""
]
},
"execution_count": 75,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"goeserie = gpd.GeoSeries([polygon, polygon2], crs=\"WGS84\")\n",
"goeserie.explore()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Getting the bounds of the grid\n",
"\n",
"Usually grids are defined in the input files of the model or inventory. \n",
"This is your job to get the data and to read it.\n",
"\n",
"In many cases the inventories and models are defined on regular lat/lon grids.\n",
"However, in some cases, we can face more exotic kinds of grids. \n",
"\n",
"In this example we will create a small regular grid."
]
},
{
"cell_type": "code",
"execution_count": 69,
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"\n",
"right_x = np.arange(1.0, 4.0)\n",
"left_x = right_x - 1.0\n",
"bottom_y = np.arange(1.0, 5.0)\n",
"top_y = bottom_y + 1.0\n",
"\n",
"\n",
"left_bot_x, left_bot_y = np.meshgrid(left_x, bottom_y)\n",
"right_top_x, right_top_y = np.meshgrid(right_x, top_y)\n",
"left_top_x, left_top_y = np.meshgrid(left_x, top_y)\n",
"right_bot_x, right_bot_y = np.meshgrid(right_x, bottom_y)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now that we have the coordinates we can create the grid. \n",
"\n",
"In general we don't want to create each polygon one by one for performance reasons.\n",
"We can use some functions to create them all at once."
]
},
{
"cell_type": "code",
"execution_count": 79,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"0 POLYGON ((0.00000 1.00000, 1.00000 1.00000, 1....\n",
"1 POLYGON ((1.00000 1.00000, 2.00000 1.00000, 2....\n",
"2 POLYGON ((2.00000 1.00000, 3.00000 1.00000, 3....\n",
"3 POLYGON ((0.00000 2.00000, 1.00000 2.00000, 1....\n",
"4 POLYGON ((1.00000 2.00000, 2.00000 2.00000, 2....\n",
"5 POLYGON ((2.00000 2.00000, 3.00000 2.00000, 3....\n",
"6 POLYGON ((0.00000 3.00000, 1.00000 3.00000, 1....\n",
"7 POLYGON ((1.00000 3.00000, 2.00000 3.00000, 2....\n",
"8 POLYGON ((2.00000 3.00000, 3.00000 3.00000, 3....\n",
"9 POLYGON ((0.00000 4.00000, 1.00000 4.00000, 1....\n",
"10 POLYGON ((1.00000 4.00000, 2.00000 4.00000, 2....\n",
"11 POLYGON ((2.00000 4.00000, 3.00000 4.00000, 3....\n",
"dtype: geometry"
]
},
"execution_count": 79,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"from shapely.creation import polygons\n",
"\n",
"geometries = np.asarray(\n",
" [\n",
" np.stack(\n",
" (\n",
" left_bot_x.ravel(),\n",
" right_bot_x.ravel(),\n",
" right_top_x.ravel(),\n",
" left_top_x.ravel(),\n",
" )\n",
" ),\n",
" np.stack(\n",
" (\n",
" left_bot_y.ravel(),\n",
" right_bot_y.ravel(),\n",
" right_top_y.ravel(),\n",
" left_top_y.ravel(),\n",
" )\n",
" ),\n",
" ]\n",
")\n",
"# Here when creatin the polygons, we need to transpose the geometries to have\n",
"# the shape = (polygons, points, coordinates)\n",
"geoserie = gpd.GeoSeries(polygons(geometries.T))\n",
"geoserie"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Creating the grid object\n",
"\n",
"emiproc provides the Grid class to create a grid object.\n",
"\n",
"If you create the geoseries of polygons, you can use the class \n",
":py:class:`emiproc.grids.GeoPandasGrid` to create the grid object.\n",
"\n",
"What is usually recommended if you need the data from a file is to create \n",
"sub class from the :py:class:`emiproc.grids.Grid` class and to implement\n",
"reading the file and generating the grid in the `__init__` method.\n"
]
},
{
"cell_type": "code",
"execution_count": 76,
"metadata": {},
"outputs": [],
"source": [
"from emiproc.grids import GeoPandasGrid\n",
"\n",
"\n",
"class MyModelGrid(GeoPandasGrid):\n",
" \"\"\"Mymodel grid.\n",
"\n",
" A grid for my model.\n",
"\n",
" Put any documentation here.\n",
" \"\"\"\n",
"\n",
" def __init__(self, grid_file_path):\n",
"\n",
" # Read the grid file\n",
" ...\n",
"\n",
" # Extract the coordinates\n",
"\n",
" coordinates_arrays = ...\n",
"\n",
" # Create the polygons\n",
"\n",
" geometries = ...\n",
" geoserie = gpd.GeoSeries(polygons(geometries), crs=\"your_crs\")\n",
"\n",
" # (optional) Get the shape of the grid\n",
" # You can use this if you have a grid given in 2d coordinates\n",
" nx, ny = ...\n",
"\n",
" # Pass the geoserie to the parent class\n",
" super().__init__(geoserie, name=\"MyModelGrid\", shape=(nx, ny))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Our example \n",
"\n",
"Let's see how this works in our example.\n"
]
},
{
"cell_type": "code",
"execution_count": 78,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"
Make this Notebook Trusted to load map: File -> Trust Notebook
"
],
"text/plain": [
""
]
},
"execution_count": 78,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"class ExampleGrid(GeoPandasGrid):\n",
" def __init__(self):\n",
" right_x = np.arange(1.0, 4.0)\n",
" left_x = right_x - 1.0\n",
" bottom_y = np.arange(5.0, 10.0)\n",
" top_y = bottom_y + 1.0\n",
"\n",
" left_bot_x, left_bot_y = np.meshgrid(left_x, bottom_y)\n",
" right_top_x, right_top_y = np.meshgrid(right_x, top_y)\n",
" left_top_x, left_top_y = np.meshgrid(left_x, top_y)\n",
" right_bot_x, right_bot_y = np.meshgrid(right_x, bottom_y)\n",
"\n",
" geometries = np.asarray(\n",
" [\n",
" np.stack(\n",
" (\n",
" left_bot_x.ravel(),\n",
" right_bot_x.ravel(),\n",
" right_top_x.ravel(),\n",
" left_top_x.ravel(),\n",
" )\n",
" ),\n",
" np.stack(\n",
" (\n",
" left_bot_y.ravel(),\n",
" right_bot_y.ravel(),\n",
" right_top_y.ravel(),\n",
" left_top_y.ravel(),\n",
" )\n",
" ),\n",
" ]\n",
" )\n",
" geoserie = gpd.GeoSeries(polygons(geometries.T), crs=\"WGS84\")\n",
"\n",
" super().__init__(geoserie, name=\"ExampleGrid\", shape=(3, 4))\n",
"\n",
"\n",
"# Finally we can create an instance of the grid\n",
"grid = ExampleGrid()\n",
"grid.gdf.explore()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Checking your grid\n",
"\n",
"As you have already seen above you can check your grid by plotting it with the explore function.\n",
"\n",
"However sometimes the grid is large and this function will be buggy.\n",
"\n",
"So we recommend selecting only a few cells to plot.\n",
"\n",
"```python\n",
"# Select the first 100 cells\n",
"grid.gdf.iloc[:100].explore()\n",
"```\n",
"\n",
"Sometimes you will see your implementation is wrong like in this example:\n",
"\n",
""
]
}
],
"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.12.2"
}
},
"nbformat": 4,
"nbformat_minor": 2
}