data/preprocessing/farm_ghg/11_Farm_emissions, nutrients and water.ipynb
{
"cells": [
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"# Exploring the use of spatial estimates of GHG emissions, nutrient application and water use"
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"Converting enviormental impact estimates associated with crop and livestock production from Halpern et al., 2022, \"The environmental footprint of global food production\" https://doi.org/10.1038/s41893-022-00965-x"
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"## Purpose"
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"Takes the impact estimates from Halpern et al., by crop group and livestock category for the year 2017 and converts these to impact intensities, which can be applied to estimate impact intensities of consumption"
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"## Methodology"
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"Take GHG emissions, excess nutrient production and water use associated with crop groups and livestock, and divide by the production to convert into an emissions intensity\n",
"\n",
"GHG emissions exclude LUC\n",
"\n",
"June 2023: Global food system pressure data was downloaded from: https://knb.ecoinformatics.org/view/doi:10.5063/F1V69H1B\n",
"\n",
"Excess N & P: \"Excess N and P nutrients are combined for the assessments in the paper. These mapped data provide excess N and P data separately for all crops (synthetic fertilizers) and reared animals (manure/excretion). Units are tonnes excess N or P (includes, leaching, runoff, volatilization), with coordinate reference system as Gall-Peters equal area with a resolution of 36km2.\". Also \"tonnes of leached/runoff/volatilized N or P from crops (synthetic fertilizers) and farmed animals (manure/excretion).\"\n",
"\n",
"Crop impacts: \"Pressure data for 26 crop categories, excluding the portion estimated to be used as animal feed. Pressures are provided per cell and include disturbance (km2eq); blue freshwater consumption (m3 water); excess nutrients (tonnes NP); and greenhouse gas emissions (tonnes CO2eq).\"\n",
"\n",
"Livestock impacts: \"Pressure data for several categories of livestock incurred at the farm site and for crop, fodder, and fish oil/meal feed. Pressures are provided per cell and include disturbance (km2eq); blue freshwater consumption (m3 water); excess nutrients (tonnes NP); and greenhouse gas emissions (tonnes CO2eq).\"\n",
"\n",
"Production data is: \"Mapped production data for crops, fisheries, livestock, mariculture. Data are in tonnes with units of production from FAO data. Coordinate reference system is Gall-Peters equal area with a resolution of 36km2.\"\n",
"\n",
"\n",
"Usage rights: Halpern et al. 2022 data is provided with the following usage rights: \"This work is dedicated to the public domain under the Creative Commons Universal 1.0 Public Domain Dedication. To view a copy of this dedication, visit https://creativecommons.org/publicdomain/zero/1.0/\""
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"## WIP - improvements"
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"Use this section only if the notebook is not final.\n",
"\n",
"Notable TODOs:\n",
"\n",
"- Todo 1;\n",
"- Todo 2;\n"
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"## Results"
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"Describe and comment the most important results."
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"## Suggested next steps"
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"State suggested next steps, based on results obtained in this notebook."
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"# Setup"
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"## Library import"
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"We import all the required PYthon libraries"
]
},
{
"cell_type": "code",
"execution_count": 129,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Requirement already satisfied: rasterio in /home/mikeha/.local/lib/python3.8/site-packages (1.3.7)\n",
"Requirement already satisfied: click>=4.0 in /usr/lib/python3/dist-packages (from rasterio) (7.0)\n",
"Requirement already satisfied: affine in /home/mikeha/.local/lib/python3.8/site-packages (from rasterio) (2.4.0)\n",
"Requirement already satisfied: numpy>=1.18 in /home/mikeha/.local/lib/python3.8/site-packages (from rasterio) (1.23.4)\n",
"Requirement already satisfied: attrs in /usr/lib/python3/dist-packages (from rasterio) (19.3.0)\n",
"Requirement already satisfied: click-plugins in /home/mikeha/.local/lib/python3.8/site-packages (from rasterio) (1.1.1)\n",
"Requirement already satisfied: certifi in /usr/lib/python3/dist-packages (from rasterio) (2019.11.28)\n",
"Requirement already satisfied: cligj>=0.5 in /home/mikeha/.local/lib/python3.8/site-packages (from rasterio) (0.7.2)\n",
"Requirement already satisfied: setuptools in /usr/lib/python3/dist-packages (from rasterio) (45.2.0)\n",
"Requirement already satisfied: snuggs>=1.4.1 in /home/mikeha/.local/lib/python3.8/site-packages (from rasterio) (1.4.7)\n",
"Requirement already satisfied: pyparsing>=2.1.6 in /home/mikeha/.local/lib/python3.8/site-packages (from snuggs>=1.4.1->rasterio) (3.0.9)\n",
"Collecting rioxarray\n",
" Downloading rioxarray-0.13.4-py3-none-any.whl (53 kB)\n",
"\u001b[K |████████████████████████████████| 53 kB 1.1 MB/s eta 0:00:01\n",
"\u001b[?25hRequirement already satisfied: rasterio>=1.1.1 in /home/mikeha/.local/lib/python3.8/site-packages (from rioxarray) (1.3.7)\n",
"Requirement already satisfied: numpy>=1.21 in /home/mikeha/.local/lib/python3.8/site-packages (from rioxarray) (1.23.4)\n",
"Requirement already satisfied: pyproj>=2.2 in /home/mikeha/.local/lib/python3.8/site-packages (from rioxarray) (3.4.0)\n",
"Collecting xarray>=0.17\n",
" Downloading xarray-2023.1.0-py3-none-any.whl (973 kB)\n",
"\u001b[K |████████████████████████████████| 973 kB 3.6 MB/s eta 0:00:01\n",
"\u001b[?25hRequirement already satisfied: packaging in /home/mikeha/.local/lib/python3.8/site-packages (from rioxarray) (21.3)\n",
"Requirement already satisfied: certifi in /usr/lib/python3/dist-packages (from rasterio>=1.1.1->rioxarray) (2019.11.28)\n",
"Requirement already satisfied: snuggs>=1.4.1 in /home/mikeha/.local/lib/python3.8/site-packages (from rasterio>=1.1.1->rioxarray) (1.4.7)\n",
"Requirement already satisfied: click>=4.0 in /usr/lib/python3/dist-packages (from rasterio>=1.1.1->rioxarray) (7.0)\n",
"Requirement already satisfied: setuptools in /usr/lib/python3/dist-packages (from rasterio>=1.1.1->rioxarray) (45.2.0)\n",
"Requirement already satisfied: affine in /home/mikeha/.local/lib/python3.8/site-packages (from rasterio>=1.1.1->rioxarray) (2.4.0)\n",
"Requirement already satisfied: cligj>=0.5 in /home/mikeha/.local/lib/python3.8/site-packages (from rasterio>=1.1.1->rioxarray) (0.7.2)\n",
"Requirement already satisfied: attrs in /usr/lib/python3/dist-packages (from rasterio>=1.1.1->rioxarray) (19.3.0)\n",
"Requirement already satisfied: click-plugins in /home/mikeha/.local/lib/python3.8/site-packages (from rasterio>=1.1.1->rioxarray) (1.1.1)\n",
"Requirement already satisfied: pandas>=1.3 in /home/mikeha/.local/lib/python3.8/site-packages (from xarray>=0.17->rioxarray) (1.5.1)\n",
"Requirement already satisfied: pyparsing!=3.0.5,>=2.0.2 in /home/mikeha/.local/lib/python3.8/site-packages (from packaging->rioxarray) (3.0.9)\n",
"Requirement already satisfied: python-dateutil>=2.8.1 in /home/mikeha/.local/lib/python3.8/site-packages (from pandas>=1.3->xarray>=0.17->rioxarray) (2.8.2)\n",
"Requirement already satisfied: pytz>=2020.1 in /home/mikeha/.local/lib/python3.8/site-packages (from pandas>=1.3->xarray>=0.17->rioxarray) (2022.5)\n",
"Requirement already satisfied: six>=1.5 in /usr/lib/python3/dist-packages (from python-dateutil>=2.8.1->pandas>=1.3->xarray>=0.17->rioxarray) (1.14.0)\n",
"Installing collected packages: xarray, rioxarray\n",
"Successfully installed rioxarray-0.13.4 xarray-2023.1.0\n"
]
}
],
"source": [
"!pip install rasterio\n",
"!pip install rioxarray"
]
},
{
"cell_type": "code",
"execution_count": 131,
"metadata": {},
"outputs": [],
"source": [
"# Data manipulation\n",
"import os\n",
"import pandas as pd\n",
"import geopandas as gpd\n",
"import numpy as np\n",
"\n",
"# Visualization\n",
"import plotly\n",
"import matplotlib\n",
"import matplotlib.pyplot as plt\n",
"import rasterio as rio\n",
"import rioxarray as riox\n",
"import xarray\n",
"\n",
"from rasterio.warp import calculate_default_transform, reproject, Resampling\n"
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"## Data import"
]
},
{
"cell_type": "code",
"execution_count": 29,
"metadata": {},
"outputs": [],
"source": [
"# Read in a list of crop super names\n",
"\n",
"root_data_path = '/mnt/c/Users/mikeha/Work/Spatial data/Environmental footprint of global food production/'\n",
"crop_table_fname = root_data_path + 'SI_SPAM_crops_tbl.csv'\n",
"crop_table = pd.read_csv(crop_table_fname,quotechar='\"')\n",
"\n",
"spam_supers = np.unique(crop_table['SPAM_super'])\n"
]
},
{
"cell_type": "code",
"execution_count": 80,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"85.47965"
]
},
"execution_count": 80,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"#Read grid cell areas for wgs84 grids (production)\n",
"dst_wgs84_grid_areas_km2 = rio.open(root_data_path + 'wgs84_grid_areas_km2.tif')\n",
"\n",
"# Read the raster data\n",
"wgs84_grid_areas_km2 = dst_wgs84_grid_areas_km2.read(1) # Assuming a single-band raster, change the index if needed\n",
"\n",
"np.nanmax(wgs84_grid_areas_km2)"
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"## Parameter definition"
]
},
{
"cell_type": "code",
"execution_count": 89,
"metadata": {},
"outputs": [],
"source": [
"# We set all relevant parameters for our notebook. (agrrements in naming convention).\n",
"\n",
"#Data paths\n",
"\n",
"prod_data_path = root_data_path + 'extra_production/'\n",
"impact_data_path = root_data_path + 'crops_food_raw/'\n",
"\n",
"impact_km2_path = root_data_path + \"crops_foods_raw_km2/\"\n",
"impact_reprojected_km2_path = root_data_path + 'crops_food_raw_km2_wgs84/'\n",
"ghg_rate_path = root_data_path + 'ghg_per_t_product_food_wgs84/'\n",
"nutrient_rate_path = root_data_path + 'nutrient_per_t_product_food_wgs84/'\n",
"\n",
"# Define the target CRS (WGS84)\n",
"target_crs = 'EPSG:4326'\n"
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"## Data processing"
]
},
{
"cell_type": "code",
"execution_count": 52,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"production no data = -3.3999999521443642e+38\n",
"ghg no data = -3.3999999521443642e+38\n",
"Nutrient no data = -3.3999999521443642e+38\n"
]
}
],
"source": [
"c = spam_supers[0]\n",
"p_src = rio.open(prod_data_path + 'crop_'+c+'_production.tif')\n",
"print('production no data = ' + str(p_src.nodata))\n",
"\n",
"ghg_src = rio.open(impact_data_path + 'gall_peter_farm_land_' + c + '_crop_produce_ghg_per_cell.tif')\n",
"print('ghg no data = ' + str(ghg_src.nodata))\n",
"\n",
"nut_src = rio.open(impact_data_path + 'gall_peter_farm_land_' + c + '_crop_produce_nutrient_per_cell.tif')\n",
"print('nutrient no data = ' + str(nut_src.nodata))"
]
},
{
"cell_type": "code",
"execution_count": 86,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"{'driver': 'GTiff',\n",
" 'dtype': 'float32',\n",
" 'nodata': -3.3999999521443642e+38,\n",
" 'width': 4736,\n",
" 'height': 3000,\n",
" 'count': 1,\n",
" 'crs': CRS.from_wkt('PROJCS[\"unknown\",GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AUTHORITY[\"EPSG\",\"4326\"]],PROJECTION[\"Cylindrical_Equal_Area\"],PARAMETER[\"standard_parallel_1\",45],PARAMETER[\"central_meridian\",0],PARAMETER[\"false_easting\",0],PARAMETER[\"false_northing\",0],UNIT[\"metre\",1,AUTHORITY[\"EPSG\",\"9001\"]],AXIS[\"Easting\",EAST],AXIS[\"Northing\",NORTH]]'),\n",
" 'transform': Affine(6000.0, 0.0, -14207430.316916062,\n",
" 0.0, -6000.0, 8999818.14485532)}"
]
},
"execution_count": 86,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"p_src.res\n",
"ghg_src.meta"
]
},
{
"cell_type": "code",
"execution_count": 127,
"metadata": {},
"outputs": [],
"source": [
"def reproject_impacts(source_raster_path,target_raster_path,impact_reprojected_km2_path,p_src,wgs84_grid_areas_km2):\n",
"\n",
" # Open the source raster\n",
" with rio.open(source_raster_path) as src:\n",
" \n",
" # convert the impact values in src to an impact per unit area (km2)\n",
" # impact pixel resolution is 6 x 6 and equal area- so divide by 36 \n",
" src_data = src.read(1, masked = True)\n",
" #src_data[src_data == src.nodata] = np.NaN\n",
" src_km2 = src_data/36.0\n",
"\n",
" #Interim output of the area rate of impact in its original projection\n",
" #with rio.open(impact_km2_path + 'test_bana.tif','w', **src.meta.copy()) as dst:\n",
" # dst.write(src_km2,indexes = 1)\n",
"\n",
" # Calculate the transformation parameters\n",
" transform, width, height = calculate_default_transform(\n",
" src_crs = src.crs,\n",
" dst_crs = target_crs,\n",
" width= src.width,height= src.height,\n",
" left=src.bounds[0], bottom=src.bounds[1], right=src.bounds[2], top=src.bounds[3],\n",
" resolution=p_src.res\n",
" )\n",
" \n",
" # Create the target raster dataset\n",
" # use p_src.width and height in order to make sure get the correct \n",
" kwargs = p_src.meta.copy()\n",
" kwargs.update({\n",
" 'crs': target_crs,\n",
" 'transform': transform,\n",
" 'width': width,\n",
" 'height': height\n",
" })\n",
"\n",
" print(kwargs)\n",
" \n",
" #create an empty array for the destination data\n",
" #dst_km2 = None # np.array(p_src.shape,np.float32)\n",
"\n",
" with rio.open(impact_reprojected_km2_path, 'w', **kwargs) as dst:\n",
"\n",
" # Reproject the source raster to the target CRS\n",
" reproject(\n",
" source=src_km2,\n",
" destination=rio.band(dst, 1),\n",
" src_transform=src.transform,\n",
" src_crs=src.crs,\n",
" dst_transform=transform,\n",
" dst_crs=target_crs,\n",
" resampling=Resampling.bilinear\n",
" )\n",
"\n",
" #convert to gridcell impact\n",
" #dst_cell = dst_km2 * wgs84_grid_areas_km2\n",
"\n",
" #dst.write(dst_cell,1)\n",
" \n",
" with rio.open(impact_reprojected_km2_path) as dst_prj_km2:\n",
" prj_data = dst_prj_km2.read(1, masked = True)\n",
" prj_data = prj_data * np.ma.masked_where(np.ma.getmask(prj_data), wgs84_grid_areas_km2[0:prj_data.shape[0],:])\n",
" with rio.open(target_raster_path, 'w', **kwargs) as dst_prj:\n",
" dst_prj.write(prj_data, 1)\n",
"\n",
" #with rio.open(impact_reprojected_km2_path + 'test_bana.tif', 'w', **kwargs) as dst:\n",
" # dst.write(dst_cell/wgs84_grid_areas_km2,indexes = 1)"
]
},
{
"cell_type": "code",
"execution_count": 157,
"metadata": {},
"outputs": [],
"source": [
"# source_raster_path = the equal area 36km2 impact data\n",
"# target_raster_path = the raster to be written out in wgs84\n",
"# match_raster_path = a raster to match the reprojection to\n",
"# impact_reprojected_km2_path = the reprojected raster in impact per km2\n",
"# wgs84_grid_areas = pixel areas (calculated in R and exported as a raster) \n",
"def reproject_impacts_riox(source_raster_path,target_raster_path,match_raster_path,impact_reprojected_km2_path,wgs84_grid_areas_km2_path):\n",
" #Open the source impact data\n",
" src_datasetreader = rio.open(source_raster_path)\n",
" src_xds = riox.open_rasterio(src_datasetreader, masked=True)\n",
" #set the crs of src file\n",
" src_xds.rio.write_crs(src_datareader.crs, inplace=True)\n",
"\n",
" #open and read the match raster file (wgs84)\n",
" match_datasetreader = rio.open(match_raster_path)\n",
" match_xds = riox.open_rasterio(match_datasetreader, masked=True)\n",
" match_xds.rio.write_crs(\"EPSG:4326\",inplace=True)\n",
"\n",
" wgs84_pixel_areas_xds = riox.open_rasterio(wgs84_grid_areas_km2_path)#,engine=\"rasterio\")\n",
" \n",
" #Convert src impact to a per km2 rate by dividing by the cell area\n",
" src_xds_km2 = src_xds/36.0\n",
"\n",
" #reproject\n",
" src_match_xds_km2 = src_xds_km2.rio.reproject_match(match_xds,resampling=Resampling.bilinear)\n",
" \n",
" #export the reprojected per km2 raster\n",
" src_match_xds_km2.rio.to_raster(impact_reprojected_km2_path)\n",
"\n",
" #convert to per cell impacts and export\n",
" src_match_xds = src_match_xds_km2 * wgs84_pixel_areas_xds\n",
" src_match_xds.rio.to_raster(target_raster_path)\n"
]
},
{
"cell_type": "code",
"execution_count": 94,
"metadata": {},
"outputs": [],
"source": [
"if os.path.exists(impact_km2_path) == False: os.mkdir(impact_km2_path)\n",
"if os.path.exists(impact_reprojected_km2_path) == False: os.mkdir(impact_reprojected_km2_path)\n",
"if os.path.exists(ghg_rate_path) == False: os.mkdir(ghg_rate_path)\n",
"if os.path.exists(nutrient_rate_path) == False: os.mkdir(nutrient_rate_path)\n"
]
},
{
"cell_type": "code",
"execution_count": 159,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"/home/mikeha/.local/lib/python3.8/site-packages/rioxarray/raster_writer.py:132: UserWarning: The nodata value (3.402823466e+38) has been automatically changed to (3.4028234663852886e+38) to match the dtype of the data.\n",
" warnings.warn(\n"
]
}
],
"source": [
"\n",
"#loop over crop types\n",
"#for c in spam_supers:\n",
"\n",
"\n",
"#read in the impact data\n",
"# Source raster path (Gall-Peters CRS)\n",
"ghg_source_raster_path = impact_data_path + 'gall_peter_farm_land_' + c + '_crop_produce_ghg_per_cell.tif'\n",
"\n",
"\n",
"# Target raster path (WGS84 CRS)\n",
"ghg_target_raster_path = ghg_rate_path + 'ghg_per_t_product_' + c + '.tif'\n",
"\n",
"\n",
"ghg_reprojected_km2_raster_path = impact_reprojected_km2_path + 'rio_ghg_per_km2_' + c + '.tif'\n",
"\n",
"# reproject_impacts(ghg_source_raster_path,\n",
"# ghg_target_raster_path,\n",
"# ghg_reprojected_km2_raster_path,\n",
"# p_src,\n",
"# wgs84_grid_areas_km2)\n",
"\n",
"\n",
"\n",
"# Target raster path (WGS84 CRS) rioxarray\n",
"rxa_ghg_target_raster_path = ghg_rate_path + 'rxa_ghg_per_t_product_' + c + '.tif'\n",
"rxa_ghg_target_km2_raster_path = impact_reprojected_km2_path + 'rxa_ghg_per_km2_' + c + '.tif'\n",
"match_raster_path = prod_data_path + 'crop_'+c+'_production.tif'\n",
"\n",
"\n",
"reproject_impacts_riox(ghg_source_raster_path,\n",
" rxa_ghg_target_raster_path,\n",
" prod_data_path + 'crop_'+c+'_production.tif',\n",
" rxa_ghg_target_km2_raster_path,\n",
" '/mnt/c/Users/mikeha/Work/Spatial data/Environmental footprint of global food production/wgs84_grid_areas_km2.tif')\n",
"\n",
"\n"
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"## References"
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"Report here relevant references:\n",
"\n",
"1. author1, article1, journal1, year1, url1\n",
"2. author2, article2, journal2, year2, url2"
]
},
{
"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.8.10"
}
},
"nbformat": 4,
"nbformat_minor": 4
}