<jupyter_start><jupyter_text>This code produces maps for a comparison of impact of potato and wheat.<jupyter_code>import os
import mypackages.myrasters as mr
import numpy as np
import fiona
import pandas as pd
import geopandas as gpd
import matplotlib
import matplotlib.pyplot as plt
from mpl_toolkits.axes_grid1 import make_axes_locatable
%matplotlib inline
data_dir = os.path.join('..', 'data')
cf_dir = os.path.join('..', 'output/characterization_factors')
out_dir = os.path.join('..', 'output/impact')
shapefile_name = 'shapefiles/countries.shp'
shapefile = os.path.join(data_dir, shapefile_name)<jupyter_output><empty_output><jupyter_text>Read in the ton-kilometers per crop:<jupyter_code>tkm = pd.read_csv(os.path.join(data_dir, 'inventory/tkm_inventory.csv'), index_col=0)<jupyter_output><empty_output><jupyter_text>Read in the characterization factor rasters:<jupyter_code>ts_raster = mr.MyRaster(os.path.join(cf_dir, 'ts_100yrs.tif'))
ms_raster = mr.MyRaster(os.path.join(cf_dir, 'ms_100yrs.tif'))
bs_raster = mr.MyRaster(os.path.join(cf_dir, 'bs_100yrs.tif'))<jupyter_output><empty_output><jupyter_text>Define a function to calculate the total impact:<jupyter_code>def total_impact_global_raster(crop):
        
    ts_array = ts_raster.get_array()
    ts = tkm[crop]['top soil'] * ts_array
    
    ts_array = None
    
    ms_array = ms_raster.get_array()
    ms = tkm[crop]['mid soil'] * ms_array
    
    ms_array = None
    
    impact = ts + ms
    
    ts = None
    ms = None
    
    bs_array = bs_raster.get_array()
    bs = tkm[crop]['bottom soil'] * bs_array
    
    impact = impact + bs
    
    impact[bs_array == bs_raster.NoDataValue] = bs_raster.NoDataValue
    
    bs_array = None
    bs = None
    
    impact = mr.mask_array(impact, bs_raster.NoDataValue)
    
    return impact<jupyter_output><empty_output><jupyter_text>Read in crops:<jupyter_code>[col for col in tkm.columns if 'Potato' in col]
[col for col in tkm.columns if 'wheat' in col]<jupyter_output><empty_output><jupyter_text>Run calculation:<jupyter_code>%%time

potato = total_impact_global_raster('Potatoes, OeLN intensive, wholesale')
wheat = total_impact_global_raster('Winter wheat TOP, OeLN intensive, wholesale')
np.percentile(potato, 98)
vmin = 0
vmax = 0.5<jupyter_output><empty_output><jupyter_text>Plot results:<jupyter_code>gdf = gpd.read_file(shapefile)
def createColourbar(lwr, upr):
    """Create a colourbar with limits of lwr and upr"""
    divider = make_axes_locatable(plt.gca())
    cax = divider.append_axes("right", size="2%", pad=0.1)
    norm = matplotlib.colors.Normalize(vmin = lwr, vmax = upr, clip = False)

    c = matplotlib.colorbar.ColorbarBase(cax, cmap='YlOrRd', norm=norm)
    return c
%%time

fig = plt.figure(figsize=(16, 9))
ax = fig.add_subplot(111)

ax.set_xlim([-180, 180])
ax.set_ylim([-56, 84])
ax.set_xticks([])
ax.set_yticks([])
ax.set_frame_on(False)

gdf.plot(ax=ax, facecolor='none', edgecolor='#000000', linewidths=0.3)

ax.imshow(potato, extent=[-180, 180, -56, 84], cmap='YlOrRd', vmin=vmin, vmax=vmax)

norm = matplotlib.colors.Normalize(vmin, vmax, clip = False)

cax = fig.add_axes([0.2, 0.25, 0.01, 0.3]) # xmin, ymin, dx, dy
matplotlib.colorbar.ColorbarBase(cax, cmap='YlOrRd', norm=norm)
cax.yaxis.set_ticks_position('left')

plt.savefig(os.path.join(out_dir, 'impact_potato.png'),
            dpi=300, bbox_inches='tight', pad_inches=0.1)
%%time

fig = plt.figure(figsize=(16, 9))
ax = fig.add_subplot(111)

ax.set_xlim([-180, 180])
ax.set_ylim([-56, 84])
ax.set_xticks([])
ax.set_yticks([])
ax.set_frame_on(False)

gdf.plot(ax=ax, facecolor='none', edgecolor='#000000', linewidths=0.3)

ax.set_xlim(-180, 180)
ax.set_ylim(-56, 84)

ax.imshow(wheat, extent=[-180, 180, -56, 84], cmap='YlOrRd', vmin=vmin, vmax=vmax)

plt.savefig(os.path.join(out_dir, 'impact_wheat.png'),
            dpi=300, bbox_inches='tight', pad_inches=0.1)
out_filename = 'impact_potato'
mr.array2geotiff_rastercopy(potato.data, os.path.join(out_dir, out_filename), bs_raster.raster)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.14359990096558553, 'fraction_numerical': 0.02921515226541223, 'mean_word_length': 5.17737003058104, 'pattern_counts': {'":': 0, '<': 29, '<?xml version=': 0, '>': 29, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '7', 'dir': 'jupyter-structured-clean-dedup', 'id': '184754', 'license': 'MIT', 'n_tokens_mistral': 1770, 'n_tokens_neox': 1628, 'n_words': 290, 'path': 'python/5b_impact_potato_vs_wheat.ipynb', 'repo_name': 'ethz-esd/compaction_stoessel_2018', 'size': '4039'}