<jupyter_start><jupyter_text>Table of Contents1&nbsp;&nbsp;MapSPAM 2010 data1.1&nbsp;&nbsp;Download data1.2&nbsp;&nbsp;Data preprocessing  MapSPAM 2010 data In this notebook we download the data that we need from [MapSPAM](http://mapspam.info/data/)The SPAM (10 x 10 km grid-cell resolution) data is specifically based on the four variables which are calculated by the model: `physical area`, `harvest area`, `production` and `yield`, for each of the 42 crops, split by `rainfed` and `irrigated` production system, as well as the total combination of both.SPAM generates a large collection of data: ~800.000 pixels X 42 crops X 4 production systems X 4 variables = ~500 million recordsPlease check the [methodology](http://mapspam.info/methodology/) section for a thorough explanation.**SPAM 2010 v1.0 Global Data (Latest)**Download pre-packaged SPAM 2010 v1.0 data and map files at the global level, as one separate file for each crop, production system, variable and format. The available formats are DBF and GeoTIFF. See the [ReadMe](https://s3.amazonaws.com/mapspam/2010/v1.0/ReadMe_v1r0_Global.txt) file first!* List of available SPAM variables:   |SPAM long name | SPAM short name |description|  |:---|:---|:---|  |PHYSICAL AREA| phys_area|Physical area is measured in hectare, and represents the actual area where a crop is grown, not counting how often production was harvested from it. Physical area is calculated for each  production system and crop, and the sum of all physical areas of the four production systems constitute the total physical area for that crop. The sum of the physical areas of all crops in a pixel may not be larger than the pixel size.|   |HARVESTED AREA|harv_area|Also measured in hectare, harvested area is at least as large as physical area, but sometimes more, since it also accounts for multiple harvests of a crop on the same plot. Like for physical area, the harvested area is calculated for each production system, and the sum of all harvested areas of all production systems in a pixel amount to the total harvested area of the pixel.|  |PRODUCTION| prod|Production, for each production system and crop, is calculated by multiplying area harvested with yield. It is measured in metric tons. The total production of a crop includes the production of all production systems of that crop.|  |YIELD| yield|Yield is a measure of productivity, the amount of production per harvested area, and is measured in kilogram/hectare. The total yield of a crop, when considering all production systems, is not the sum of the individual yields, but the weighted average of the 4 yields|  * List of available [crops](http://mapspam.info/wp-content/uploads/4-Methodology-Crops-of-SPAM-2005-2015-02-26.csv) |No. crt.| SPAM short name | SPAM long name |FAONAMES | GROUP||:---|:---|:---|:---|:---|  |1 | whea | wheat                | wheat                                | cereals                       ||2 | rice | rice                 | rice                                 | cereals                       ||3 | maiz | maize                | maize                                | cereals                       ||4 | barl | barley               | barley                               | cereals                       ||5 | pmil | pearl millet         | millet                               | cereals                       ||6 | smil | small millet         | millet                               | cereals                       ||7 | sorg | sorghum              | sorghum                              | cereals                       ||8 | ocer | other cereals        | other cereals ++                     | cereals                       ||9 | pota | potato               | potato                               | roots&tubers or starchy roots ||10| swpo | sweet potato         | sweet potato                         | roots&tubers or starchy roots ||11| yams | yams                 | yam                                  | roots&tubers or starchy roots ||12| cass | cassava              | cassava                              | roots&tubers or starchy roots ||13| orts | other roots          | yautia ++                            | roots&tubers or starchy roots ||14| bean | bean                 | beans, dry                           | pulses                        ||15| chic | chickpea             | chickpea                             | pulses                        ||16| cowp | cowpea               | cowpea                               | pulses                        ||17| pige | pigeonpea            | pigeon pea                           | pulses                        ||18| lent | lentil               | lentils                              | pulses                        ||19| opul | other pulses         | broad beans ++                       | pulses                        ||20| soyb | soybean              | soybean                              | oilcrops                      ||21| grou | groundnut            | groundnut, with shell                | oilcrops                      ||22| cnut | coconut              | coconut                              | oilcrops                      ||23| oilp | oilpalm              | palmoil                              | oilcrops                      ||24| sunf | sunflower            | sunflower seed                       | oilcrops                      ||25| rape | rapeseed             | rapeseed                             | oilcrops                      ||26| sesa | sesameseed           | sesame seed                          | oilcrops                      ||27| ooil | other oil crops      | olives ++                            | oilcrops                      ||28| sugc | sugarcane            | sugar cane                           | sugar crops                   ||29| sugb | sugarbeet            | sugarbeet                            | sugar crops                   ||30| cott | cotton               | seed cotton                          | fibres                        ||31| ofib | other fibre crops    | other fibres ++                      | fibres                        ||32| acof | arabica coffee       | coffee                               | stimulant                     ||33| rcof | robusta coffee       | coffee                               | stimulant                     ||34| coco | cocoa                | cocoa                                | stimulant                     ||35| teas | tea                  | tea                                  | stimulant                     ||36| toba | tobacco              | tobacco leaves                       | stimulant                     ||37| bana | banana               | banana                               | fruits                        ||38| plnt | plantain             | plantain                             | fruits                        ||39| trof | tropical fruit       | oranges ++                           | fruits                        ||40| temf | temperate fruit      | apples ++                            | fruits                        ||41| vege | vegetables           | cabbages and other brassicas ++      | vegetables                    ||42| rest | rest of crops        | all individual other crops           |                               |* List of available production systems|name| production systems ||:---|:---|  |Irrigated|irrigated portion of crop|  |Rainfed|rainfed portion of crop|   |Rainfed|rainfed high inputs portion of crop| |Rainfed|rainfed low inputs portion of crop| |Total | all technologies together, ie complete crop| **SPAM 2010 data structure:**The SPAM 2010 data files are structured like this:- Base url: https://s3.amazonaws.com/mapspam/2010/v1.0/dbf/- Zip files: |File name| description ||:---|:---| |spam2010v1r0_global_harv_area.dbf.zip	|SPAM area harvested, global pixels, files in 2 formats x 6 technologies, strucuture A, record type H ||spam2010v1r0_global_prod.dbf.zip		|SPAM production, global pixels, files in 2 formats x 6 technologies, strucuture A, record type P     |- File names:     All files have standard names, which allow direct identification of variable and technology:    spam2010v1r0_global_v_t.dbf        where        v = variable        t = technology        - Variables (v):        | v  | SPAM name | |:---|:---|| harvested-area | harvested area || production | production |        - Technologies (t):     | t  | description | |:---|:---||ti | irrigated portion of crop||tr | rainfed portion of crop  |<jupyter_code>import getpass
import requests
from tqdm import tqdm
import pandas as pd
import geopandas as gpd
import numpy as np
import os
import zipfile
from simpledbf import Dbf5
import fiona
import rtree
from collections import OrderedDict
from shapely.geometry import mapping, shape
import matplotlib.pyplot as plt
import shapely.wkb 
from shapely.ops import cascaded_union
from shapely.geometry import Polygon, Point, MultiPolygon<jupyter_output><empty_output><jupyter_text>Download dataWe download the `harvested area` and `production` variables for 2010<jupyter_code>baseFolder = '/home/jovyan/work/datasets'
baseUrl = 'https://s3.amazonaws.com/mapspam/2010/v1.0/dbf/'
folder = f'{baseFolder}/aqueduct/data_source/spamdata2010/'
variables = ['harv_area', 'prod']

for variable in variables:
    filename = 'spam2010v1r0_global'+'_'+variable+'.dbf.zip'
    dataUrl=baseUrl+filename
    response = requests.get(dataUrl, stream=True)
    
    file_output = folder + filename
    print(filename)
    with open(file_output, 'wb') as handle:
        for chunk in tqdm(response.iter_content(chunk_size=128)):
            handle.write(chunk)
            
    # uncompress zip file
    with zipfile.ZipFile(file_output,"r") as zip_ref:
        zip_ref.extractall(folder)
      
    # remove zip file
    os.remove(file_output)<jupyter_output><empty_output><jupyter_text>Data preprocessing Read dataWe read only the next - Variables:    - Harvested area    - Production - Production systems:    - Irrigated    - Rainfed     and merge them all together in a single table<jupyter_code>baseFolder = f'{baseFolder}/aqueduct/data_source/spamdata2010/'
variables_file = {'harv_area': 'harvested-area', 'prod': 'production'}
production_types = {'harv_area': 'area', 'prod': 'prod'}

irrigation = {'irrigated': 'ti', 'rainfed': 'tr'}
suffix = {'irrigated': 'i', 'rainfed': 'r'}

crops = {"whea": "wheat", "rice": "rice", "maiz": "maize", "barl": "barley", "pmil": "pearl millet", 
         "smil": "small millet", "sorg": "sorghum", "ocer": "other cereals", "pota": "potato", 
         "swpo": "sweet potato", "yams": "yams", "cass": "cassava", "orts": "other roots", 
         "bean": "bean", "chic": "chickpea", "cowp": "cowpea", "pige": "pigeonpea", "lent": "lentil", 
         "opul": "other pulses", "soyb": "soybean", "grou": "groundnut", "cnut": "coconut", 
         "oilp": "oilpalm", "sunf": "sunflower", "rape": "rapeseed", "sesa": "sesameseed", 
         "ooil": "other oil crops", "sugc": "sugarcane", "sugb": "sugarbeet", "cott": "cotton", 
         "ofib": "other fibre crops", "acof": "arabica coffee", "rcof": "robusta coffee", 
         "coco": "cocoa", "teas": "tea", "toba": "tobacco", "bana": "banana", "plnt": "plantain", 
         "trof": "tropical fruit", "temf": "temperate fruit", "vege": "vegetables", "rest": "rest of crops"}                           


df_prod = pd.DataFrame(columns=['iso', 'prod_level', 'cell5m', 'name_cntr', 'name_adm1', 'name_adm2', 'unit_prod', 'crop', 'prod', 'irrigation'])
df_area = pd.DataFrame(columns=['iso', 'prod_level', 'cell5m', 'name_cntr', 'name_adm1', 'name_adm2', 'unit_area', 'crop', 'area', 'irrigation'])

for itype in irrigation.keys():
    
    df_type = pd.DataFrame(columns=['cell5m', 'crop', 'irrigation', 'area', 'prod', 
                               'unit_area', 'unit_prod', 'iso', 'name_cntr',
                               'prod_level'
                              ])
    
    for variable in variables_file.keys():
        folder = 'spam2010v1r0_global'+'_'+variable+'.dbf'
        filename = 'spam2010v1r0_global'+'_'+variables_file[variable]+'_'+irrigation[itype]+'.dbf'
        
        print(filename)
        dbf = Dbf5(baseFolder+folder+'/'+filename, codec='cp1252')
        df = dbf.to_dataframe()
        df.columns = df.columns.str.lower()
        
        # Select columns
        columns = list(dict(("{}_{}".format(k,suffix[itype]),v) for k,v in crops.items()))
        columns.extend(['iso3', 'prod_level', 'cell5m', 'name_cntr', 'name_adm1', 'name_adm2', 'unit'])
        df = df[columns]

        # Change columns names
        change_names = dict(("{}_{}".format(k,suffix[itype]),v) for k,v in crops.items())
        change_names.update({'iso3': 'iso', 'unit': 'unit_'+production_types[variable]})
        df.rename(columns=change_names, inplace=True)

        # Unpivot DataFrame
        df = pd.melt(df, id_vars=['iso', 'prod_level', 'cell5m', 'name_cntr', 'name_adm1', 'name_adm2', 'unit_'+production_types[variable]], 
                          var_name="crop", value_name=production_types[variable])

        # Sort by
        df.sort_values(['iso', 'cell5m'], inplace=True)
        df.sort_values(['iso', 'cell5m'], inplace=True)
        
        # Add irrigation column
        df['irrigation'] = itype
        
        if variable == 'prod':
            df_prod = pd.concat([df_prod, df])
        else:
            df_area = pd.concat([df_area, df])         

# Merge Harvested area and Production DataFrames
df_all = pd.merge(left=df_area, right=df_prod, how='left',
                  on=['iso', 'prod_level', 'cell5m', 'name_cntr', 'name_adm1', 'name_adm2', 'crop', 'irrigation'])<jupyter_output><empty_output><jupyter_text>Filter dataWe will select all cells that have `Harvested area` bigger or equal to 10 ha.<jupyter_code>df_filtered = df_all[df_all['area'] >= 10]<jupyter_output><empty_output><jupyter_text>Save table<jupyter_code>df_filtered.to_csv(f'{baseFolder}/aqueduct/dst/spam2010/data_filtered.csv', sep=',')<jupyter_output><empty_output><jupyter_text>Intersection between the SPAM grid cells and the basinsWe will perform the clip between the `cell5M` geometries and the basins.If a `cell5M` falls into more than one basin then we will split the cell and compute the percentage of the cell covered by each basin. **Read `cell5M` geometries**<jupyter_code>cell5m = gpd.read_file(f'{baseFolder}/aqueduct/data_source/spamdata/spam_cells/April_2013_IMPACT_Cell5M.shp')
cell5m.columns = [x.lower() for x in cell5m.columns]
cell5m.rename(columns={'cell5m_id_': 'cell5m', 'iso3code': 'iso'}, inplace=True)
cell5m.sort_values(['iso', 'cell5m'], inplace=True)

# convert the 'cell5m' column to object
cell5m['cell5m'] = cell5m['cell5m'].astype('int') 
#cell5m['cell5m'] = cell5m['cell5m'].astype('object') 

cell5m.drop(columns=['iso', 'alloc_key'], inplace=True)

cell5m.head()<jupyter_output><empty_output><jupyter_text>**Read basins for baseline water risk indicators**<jupyter_code>baseline_geo = gpd.read_file(f'{baseFolder}/aqueduct/data_source/AQ_2_water_risk_atlas/Y2018M12D06_RH_Master_Shape_V01/output_V02/Y2018M12D06_RH_Master_Shape_V01.shp')
baseline_geo.drop(columns=['aqid', 'gid_1', 'pfaf_id', 'type'], inplace=True)
baseline_geo['geometry'] = baseline_geo['geometry'].apply(lambda x: x.buffer(0))
baseline_geo.to_file(f'{baseFolder}/aqueduct/dst/spam2010/baseline_geo/baseline_geo.shp')
baseline_geo.head()<jupyter_output><empty_output><jupyter_text>**Read basins for projected water risk indicators**<jupyter_code>projected_geo = gpd.read_file(f'{baseFolder}/aqueduct/data_source/AQ_2_water_risk_atlas/wri_subcatchements/wri_subcatchements.shp')
projected_geo.drop(columns=['cartodb_id', 'gu', 'country', 'basin_name', 'iso', 'area_km2', 'down_basin'], inplace=True)
projected_geo = projected_geo[projected_geo.geometry.notnull()]
projected_geo['geometry'] = projected_geo['geometry'].apply(lambda x: x.buffer(0))
projected_geo.to_file(f'{baseFolder}/aqueduct/dst/spam2010/projected_geo/projected_geo.shp')
projected_geo.head()<jupyter_output><empty_output><jupyter_text>**We use geopandas’s R-tree spatial index to find which basins lie within each grid cell** Example of an individual cells<jupyter_code>def rtree_intersect_cell(gdf, grid, cell5m):
    columns = list(gdf.columns)
    columns.append('percentage')
    gdf_new = gpd.GeoDataFrame(columns=columns)
    
    cell = grid[grid['cell5m'] == cell5m].geometry.iloc[0]
    sindex = gdf.sindex

    # basins that intersect with the cell
    possible_matches_index = list(sindex.intersection(cell.bounds))
    possible_matches = gdf.iloc[possible_matches_index]
    # intersection between the basin and the cell
    intersection = possible_matches.intersection(cell)
    # percentage of the cell covered by each basin
    percentage = intersection.apply(lambda x: x.area/cell.area)
    # precise matches where percentage > 0
    precise_matches = possible_matches[possible_matches.index.isin(percentage[percentage > 0].index)]

    for column in gdf_new.columns:
        if (column != 'geometry') and (column != 'percentage'):
            gdf_new[column]=precise_matches[column]
        else:
            if column == 'geometry':
                gdf_new[column]=intersection.loc[precise_matches.index]
            if column == 'percentage':
                gdf_new[column]=percentage.loc[precise_matches.index]
            
    fig, ax = plt.subplots(figsize=[5,5])
    ax.set_aspect('equal')
    
    possible_matches.plot(ax=ax, edgecolor='r')
    precise_matches.plot(ax=ax, color='g',  edgecolor='r')
    gdf_new.plot(ax=ax,color='r', edgecolor='w')
    gpd.GeoDataFrame(data={'geometry': [cell]}).plot(ax=ax,color='w', edgecolor='k', alpha=0.1)
    plt.xlim(cell.bounds[0]-0.2,cell.bounds[2]+0.2)
    plt.ylim(cell.bounds[1]-0.2,cell.bounds[3]+0.2);
rtree_intersect_cell(baseline_geo, cell5m, 3065782)<jupyter_output><empty_output><jupyter_text>We iterate over the cells<jupyter_code>def rtree_intersect(gdf, grid, path):
    columns = list(gdf.columns)
    columns.append('percentage')
    columns.append('cell5m')
    gdf_all = gpd.GeoDataFrame(columns=columns)
    sindex = gdf.sindex
    # we iterate over the cells
    for n, cell in enumerate(tqdm(grid.geometry)):
        gdf_new = gpd.GeoDataFrame(columns=columns)
        # basins that intersect with the cell
        possible_matches_index = list(sindex.intersection(cell.bounds))
        possible_matches = gdf.iloc[possible_matches_index]
        # intersection between the basin and the cell
        intersection = possible_matches.intersection(cell)
        # percentage of the cell covered by each basin
        percentage = intersection.apply(lambda x: x.area/cell.area)
        # precise matches where percentage > 0
        precise_matches = possible_matches[possible_matches.index.isin(percentage[percentage > 0].index)]
        
        for column in gdf_new.columns:
            if (column != 'geometry') and (column != 'percentage') and (column != 'cell5m'):
                gdf_new[column]=precise_matches[column]
            else:
                if column == 'geometry':
                    gdf_new[column]=intersection.loc[precise_matches.index]
                if column == 'percentage':
                    gdf_new[column]=percentage.loc[precise_matches.index]
                if column == 'cell5m':
                    gdf_new[column]=grid['cell5m'].iloc[n]
                    
        gdf_all = pd.concat([gdf_all, gdf_new])
        if n % 100000 == 0:
            gdf_all.to_file(path)
    
    gdf_all.to_file(path)
    return gdf_all
def rtree_intersect_fiona(gdf_path, grid, out_path):
    gdf = gpd.read_file(gdf_path)
    columns = list(gdf.columns)
    columns.append('percentage')
    columns.append('cell5m')

    properties = columns.copy()
    properties.remove('geometry')

    sindex = gdf.sindex
    
    with fiona.open(gdf_path, 'r') as layer1:
        # We copy schema and add the  new property for the new resulting shp
        crss=layer1.crs
        schema = layer1.schema.copy()
        for property in properties: 
            schema['properties'][property] = 'str:80'
    
        with fiona.open(out_path, 'w', 'ESRI Shapefile', schema, crs=crss) as layer2:
    
            # we iterate over the cells
            for n, cell in enumerate(tqdm(grid.geometry)):
                gdf_new = gpd.GeoDataFrame(columns=columns)
                # basins that intersect with the cell
                possible_matches_index = list(sindex.intersection(cell.bounds))
                possible_matches = gdf.iloc[possible_matches_index]
                # intersection between the basin and the cell
                intersection = possible_matches.intersection(cell)
                # percentage of the cell covered by each basin
                percentage = intersection.apply(lambda x: x.area/cell.area)
                # precise matches where percentage > 0
                precise_matches = possible_matches[possible_matches.index.isin(percentage[percentage > 0].index)]
        
                for column in gdf_new.columns:
                    if (column != 'geometry') and (column != 'percentage') and (column != 'cell5m'):
                        gdf_new[column]=precise_matches[column]
                    else:
                        if column == 'geometry':
                            gdf_new[column]=intersection.loc[precise_matches.index]
                        if column == 'percentage':
                            gdf_new[column]=percentage.loc[precise_matches.index]
                        if column == 'cell5m':
                            gdf_new[column]=grid['cell5m'].iloc[n]
                        
                gdf_new[properties] = gdf_new[properties].astype('object')
                for i in range(len(gdf_new)):
                    # Add the content to the right schema in the new shp
                    layer2.write({
                        'properties': OrderedDict(gdf_new[properties].iloc[i].to_dict()),
                        'geometry': mapping(gdf_new['geometry'].iloc[i])
                    });
cell5mPath = f'{baseFolder}/aqueduct/data_source/spamdata/spam_cells/April_2013_IMPACT_Cell5M.shp'
from tqdm import tqdm
import fiona
import rtree
from shapely.geometry import mapping, shape
import shapely.wkb 
from shapely.ops import cascaded_union
from shapely.geometry import Polygon, Point, MultiPolygon
from multiprocessing import Pool

def intersection(layer2, layer1, layer3, index):
    for feat2 in tqdm(layer2):
                    geom2 = shape(feat2['geometry'])
                    for fid in list(index.intersection(geom2.bounds)):
                        if fid != int(feat2['id']):
                            feat1 = layer1[fid]
                            geom1 = shape(feat1['geometry'])
                            if geom1.intersects(geom2):
                                # Add the content to the right schema in the new shp
                                layer3.write({
                                    'properties': {**feat2['properties'],**feat1['properties']},
                                    'geometry': mapping(geom1.intersection(geom2))
                                });
    
def rtree_fiona_nondiscrimination(gdf_path, grid, out_path):
    with fiona.open(gdf_path, 'r') as layer1:
        with fiona.open(grid, 'r') as layer2:
            # We copy schema and add the  new property for the new resulting shp
            crss=layer2.crs
            schema = layer2.schema.copy()
            schema['properties'] = {**schema['properties'], **layer1.schema.copy()['properties']}
            # We open a first empty shp to write new content from both others shp
            with fiona.open(out_path, 'w', 'ESRI Shapefile', schema, crs=crss) as layer3:
                index = rtree.index.Index()
                for feat1 in tqdm(layer1):
                    fid = int(feat1['id'])
                    geom1 = shape(feat1['geometry'])
                    index.insert(fid, geom1.bounds)

                for feat2 in tqdm(layer2):
                    geom2 = shape(feat2['geometry'])
                    for fid in list(index.intersection(geom2.bounds)):
                        if fid != int(feat2['id']):
                            feat1 = layer1[fid]
                            geom1 = shape(feat1['geometry'])
                            if geom1.intersects(geom2):
                                # Add the content to the right schema in the new shp
                                layer3.write({
                                    'properties': {**feat2['properties'],**feat1['properties']},
                                    'geometry': mapping(geom1.intersection(geom2))
                                });
cell5mPath = f'{baseFolder}/aqueduct/data_source/spamdata/spam_cells/April_2013_IMPACT_Cell5M.shp'
gdf_path = f'{baseFolder}/aqueduct/dst/spam2010/baseline_geo/baseline_geo.shp'
out_path = f'{baseFolder}/aqueduct/dst/spam2010/cell5m_baseline/cell5m_baseline_v2.shp'
rtree_fiona_nondiscrimination(gdf_path, cell5mPath, out_path)
gdf_path = f'{baseFolder}/aqueduct/dst/spam2010/projected_geo/projected_geo.shp'
out_path = f'{baseFolder}/data/aqueduct/dst/spam2010/cell5m_projected/cell5m_projected_v2.shp'
rtree_fiona_nondiscrimination(gdf_path, cell5mPath, out_path)<jupyter_output><empty_output><jupyter_text>**Read `cell5M` geometries intersected with baseline basins**<jupyter_code>cell5M_baseline = gpd.read_file(f'{baseFolder}/aqueduct/dst/spam2010/cell5m_baseline/cell5m_baseline_v2.shp')
cell5M_baseline['aq30_id'] = cell5M_baseline['aq30_id'].astype('int64')
cell5M_baseline['percentage'] = cell5M_baseline['percentage'].astype('float64')
cell5M_baseline['percentage'] = cell5M_baseline['percentage'].apply(lambda x: round(x,4))
cell5M_baseline['cell5m'] = cell5M_baseline['cell5m'].astype('int')
cell5M_baseline.head()<jupyter_output><empty_output><jupyter_text>**Read `cell5M` geometries intersected with projected basins**<jupyter_code>cell5M_projected = gpd.read_file(f'{baseFolder}/aqueduct/dst/spam2010/cell5m_projected/cell5m_projected_v2.shp')
cell5M_projected['basinid'] = cell5M_projected['basinid'].astype('int64')
cell5M_projected['percentage'] = cell5M_projected['percentage'].astype('float64')
cell5M_projected['percentage'] = cell5M_projected['percentage'].apply(lambda x: round(x,4))
cell5M_projected['cell5m'] = cell5M_projected['cell5m'].astype('int')
cell5M_projected.head()<jupyter_output><empty_output><jupyter_text>Add `cell5M` geometry and baisin ids to the SPAM data table<jupyter_code>df_filtered = pd.read_csv(f'{baseFolder}/aqueduct/dst/spam2010/data_filtered.csv')
df_filtered.drop(columns=['Unnamed: 0'], inplace=True)
gdf_filtered = gpd.GeoDataFrame(df_filtered)
gdf_filtered['cell5m'] = gdf_filtered['cell5m'].astype('int')
gdf_filtered.head()<jupyter_output><empty_output><jupyter_text>**Add baseline basins**<jupyter_code># Merge GeoDataFrames
filtered_cell5M_baseline = pd.merge(left=gdf_filtered, right=cell5M_baseline, how='left', on='cell5m')

# Multiply the area and the production with the percentage of each cell fraction
def f(*x):
    return x[0]*x[1]
filtered_cell5M_baseline['prod'] = filtered_cell5M_baseline[['prod','percentage']].apply(lambda x: f(*x), axis=1).round(1)
filtered_cell5M_baseline['area'] = filtered_cell5M_baseline[['area','percentage']].apply(lambda x: f(*x), axis=1).round(1)
filtered_cell5M_baseline.drop(columns=['percentage'], inplace=True)

# Drop duplicates
filtered_cell5M_baseline.drop_duplicates(subset=['cell5m', 'crop', 'irrigation', 'prod', 'aq30_id'], inplace=True)

# Save GeoDataFrame
#filtered_cell5M_baseline.to_file('/Volumes/MacBook HD/data/aqueduct/dst/spam2010/data_baseline/data_baseline.shp')<jupyter_output><empty_output><jupyter_text>**Add projected basins**<jupyter_code># Merge GeoDataFrames
filtered_cell5M_projected = pd.merge(left=gdf_filtered, right=cell5M_projected, how='left', on='cell5m')

# Multiply the area and the production with the percentage of each cell fraction
def f(*x):
    return x[0]*x[1]
filtered_cell5M_projected['prod'] = filtered_cell5M_projected[['prod','percentage']].apply(lambda x: f(*x), axis=1).round(1)
filtered_cell5M_projected['area'] = filtered_cell5M_projected[['area','percentage']].apply(lambda x: f(*x), axis=1).round(1)
filtered_cell5M_projected.drop(columns=['percentage'], inplace=True)

# Drop duplicates
filtered_cell5M_projected.drop_duplicates(subset=['cell5m', 'crop', 'irrigation', 'prod', 'basinid'], inplace=True)

# Save GeoDataFrame
#filtered_cell5M_projected.to_file('/Volumes/MacBook HD/data/aqueduct/dst/spam2010/data_projected/data_projected.shp')<jupyter_output><empty_output><jupyter_text>Ranking of the crops in each cellFor a given `cell5m`, basin and irrigation type we perform a ranking of the most productive crops, where 1 corresponds to the most productive one. **Baseline basins**<jupyter_code>rank_baseline = filtered_cell5M_baseline[['cell5m', 'crop', 'irrigation', 'prod', 'aq30_id']].copy()

rank_baseline_new = pd.DataFrame(columns=['cell5m', 'crop', 'irrigation', 'prod', 'aq30_id', 'rank'])

for irrigation in ['irrigated', 'rainfed']:
    rank = rank_baseline[rank_baseline['irrigation'] == irrigation]
    
    rank.sort_values(['cell5m'], inplace=True)
    cell5m_unique = rank['cell5m'].unique()
    cell5m_unique.sort()
    
    cell = []
    basin = []
    prod = []
    ranking = []

    for cell5m in  tqdm(cell5m_unique):
        df = rank[rank['cell5m'] == cell5m]
        for aq30_id in df['aq30_id'].unique():
            df_aq = df[df['aq30_id'] == aq30_id]
        
            cell.extend([cell5m] * len(df_aq))
            basin.extend([aq30_id] * len(df_aq))
            prod.extend(list(df_aq['prod']))
            ranking.extend(list(df_aq['prod'].rank(ascending=False)))
        
    df_rank = pd.DataFrame({'cell5m': cell, 'aq30_id': basin, 'prod': prod, 'rank': ranking}) 
    

    rank_new = pd.merge(left=rank, right=df_rank, how='inner', on=['cell5m', 'aq30_id', 'prod'])   
    
    rank_baseline_new = pd.concat([rank_baseline_new, rank_new])
    
    rank_baseline_new.to_csv('/Volumes/MacBook HD/data/aqueduct/dst/spam2010/rank_baseline.csv')
rank_baseline = pd.read_csv('/Volumes/MacBook HD/data/aqueduct/dst/spam2010/rank_baseline.csv')
rank_baseline.drop(columns='Unnamed: 0', inplace=True)
rank_baseline['rank'] = rank_baseline['rank'].astype('int64')
rank_baseline.head(1)
filtered_cell5M_baseline.head(1)
filtered_cell5M_rank_baseline = pd.merge(left=filtered_cell5M_baseline, right=rank_baseline, on=['cell5m', 'crop', 'irrigation', 'prod', 'aq30_id'])
filtered_cell5M_rank_baseline['rank'] = filtered_cell5M_rank_baseline['rank'].astype('int64')
filtered_cell5M_rank_baseline['aq30_id'] = filtered_cell5M_rank_baseline['aq30_id'].astype('int64')
filtered_cell5M_rank_baseline.drop(columns=['name_adm1', 'name_adm2'], inplace=True)
filtered_cell5M_rank_baseline.head(1)<jupyter_output><empty_output><jupyter_text>Save table<jupyter_code>filtered_cell5M_rank_baseline.to_file('/Volumes/MacBook HD/data/aqueduct/dst/spam2010/cell5m_baseline_rank/filtered_cell5M_rank_baseline.shp')<jupyter_output><empty_output><jupyter_text>**Projected basins**<jupyter_code>rank_projected = filtered_cell5M_projected[['cell5m', 'crop', 'irrigation', 'prod', 'basinid']].copy()

rank_projected_new = pd.DataFrame(columns=['cell5m', 'crop', 'irrigation', 'prod', 'basinid', 'rank'])

for irrigation in ['irrigated', 'rainfed']:
    rank = rank_projected[rank_projected['irrigation'] == irrigation]
    
    rank.sort_values(['cell5m'], inplace=True)
    cell5m_unique = rank['cell5m'].unique()
    cell5m_unique.sort()
    
    cell = []
    basin = []
    prod = []
    ranking = []

    for cell5m in  tqdm(cell5m_unique):
        df = rank[rank['cell5m'] == cell5m]
        for basinid in df['basinid'].unique():
            df_aq = df[df['basinid'] == basinid]
        
            cell.extend([cell5m] * len(df_aq))
            basin.extend([basinid] * len(df_aq))
            prod.extend(list(df_aq['prod']))
            ranking.extend(list(df_aq['prod'].rank(ascending=False)))
        
    df_rank = pd.DataFrame({'cell5m': cell, 'basinid': basin, 'prod': prod, 'rank': ranking}) 
    

    rank_new = pd.merge(left=rank, right=df_rank, how='inner', on=['cell5m', 'basinid', 'prod'])   
    
    rank_projected_new = pd.concat([rank_projected_new, rank_new])
    
    rank_projected_new.to_csv('/Volumes/MacBook HD/data/aqueduct/dst/spam2010/rank_projected.csv')
rank_projected = pd.read_csv('/Volumes/MacBook HD/data/aqueduct/dst/spam2010/rank_projected.csv')
rank_projected.drop(columns='Unnamed: 0', inplace=True)
rank_projected['rank'] = rank_projected['rank'].astype('int64')
rank_projected.head(1)
filtered_cell5M_projected.head(1)
filtered_cell5M_rank_projected = pd.merge(left=filtered_cell5M_projected, right=rank_projected, on=['cell5m', 'crop', 'irrigation', 'prod', 'basinid'])
filtered_cell5M_rank_projected['rank'] = filtered_cell5M_rank_projected['rank'].astype('int64')
filtered_cell5M_rank_projected['basinid'] = filtered_cell5M_rank_projected['basinid'].astype('int64')
filtered_cell5M_rank_projected.drop(columns=['name_adm1', 'name_adm2'], inplace=True)
filtered_cell5M_rank_projected.head(1)<jupyter_output><empty_output><jupyter_text>Save table<jupyter_code>filtered_cell5M_rank_projected.to_file('/Volumes/MacBook HD/data/aqueduct/dst/spam2010/cell5m_projected_rank/filtered_cell5M_rank_projected.shp')<jupyter_output><empty_output><jupyter_text>Split tableWe split the tables in x datasets to make it uploable to a carto account<jupyter_code>def split_table(basefile, filename, size):
    """
    this function splits a table in x datasets to make it uploable to a carto account
    """
    dataset = gpd.read_file(basefile+filename+'.shp')
    for i in range(int(np.ceil(dataset.shape[0]/size))):
        split_dataset_name = filename+'_'+str(i)
        print(split_dataset_name)
        os.mkdir(basefile+split_dataset_name)
        full_path = basefile+split_dataset_name+'/'+split_dataset_name+'.shp'
        dataset.iloc[i*size:(i+1)*size].to_file(full_path)
basefile = '/Volumes/MacBook HD/data/aqueduct/dst/spam2010/cell5m_baseline_rank/'
split_table(basefile, 'filtered_cell5M_rank_baseline', 2000000)
basefile = '/Volumes/MacBook HD/data/aqueduct/dst/spam2010/cell5m_projected_rank/'
split_table(basefile, 'filtered_cell5M_rank_projected', 2000000)<jupyter_output><empty_output><jupyter_text>Merge table in carto **Baseline basins** ```sqlSELECT the_geom, aq30_id, area, cell5m, crop, irrigation, iso, name_cntr, prod, prod_level, rank, unit_area, unit_prod from filtered_cell5m_rank_baseline_0union allSELECT the_geom, aq30_id, area, cell5m, crop, irrigation, iso, name_cntr, prod, prod_level, rank, unit_area, unit_prod from filtered_cell5m_rank_baseline_1union allSELECT the_geom, aq30_id, area, cell5m, crop, irrigation, iso, name_cntr, prod, prod_level, rank, unit_area, unit_prod from filtered_cell5m_rank_baseline_2union allSELECT the_geom, aq30_id, area, cell5m, crop, irrigation, iso, name_cntr, prod, prod_level, rank, unit_area, unit_prod from filtered_cell5m_rank_baseline_3union allSELECT the_geom, aq30_id, area, cell5m, crop, irrigation, iso, name_cntr, prod, prod_level, rank, unit_area, unit_prod from filtered_cell5m_rank_baseline_4``` **Projected basins** ```sqlSELECT the_geom, basinid, area, cell5m, crop, irrigation, iso, name_cntr, prod, prod_level, rank, unit_area, unit_prod from filtered_cell5m_rank_projected_0union allSELECT the_geom, basinid, area, cell5m, crop, irrigation, iso, name_cntr, prod, prod_level, rank, unit_area, unit_prod from filtered_cell5m_rank_projected_1union allSELECT the_geom, basinid, area, cell5m, crop, irrigation, iso, name_cntr, prod, prod_level, rank, unit_area, unit_prod from filtered_cell5m_rank_projected_2union allSELECT the_geom, basinid, area, cell5m, crop, irrigation, iso, name_cntr, prod, prod_level, rank, unit_area, unit_prod from filtered_cell5m_rank_projected_3```  Create table with each crop bucket ```sqlwith r as (SELECT sum(prod) val, cell5m, iso, crop, irrigation            FROM crops_projected group by cell5m, iso, crop, irrigation),      s as (SELECT sum(prod) val, cell5m, iso, crop            FROM crops_projected group by cell5m, iso, crop),     t as (SELECT sum(prod) val, cell5m, crop, irrigation            FROM crops_projected group by cell5m, crop, irrigation),      tt as (SELECT sum(prod) val, cell5m, crop             FROM crops_projected group by cell5m, crop),      ff as (select val, cell5m, iso, crop, irrigation from r             union all             select val, cell5m, iso, crop, 'all' as irrigation from s             union all             select val, cell5m, 'WORLD' as iso, crop, irrigation  FROM t             union all             select val, cell5m, 'WORLD' as iso, crop, 'all' as irrigation from tt) SELECT CDB_JenksBins(array_agg(val::numeric),5) as bucket, crop, irrigation, iso from ff where val is not null group by crop, irrigation, iso```<jupyter_code>!pip install carto
import requests
from carto.auth import APIKeyAuthClient
from carto.sql import BatchSQLClient
from carto.sql import SQLClient
import getpass
import time
query =("with r as (SELECT sum(prod) val, cell5m, iso, crop, irrigation FROM crops_projected group by cell5m, iso, crop, irrigation), s as (SELECT sum(prod) val, cell5m, iso, crop FROM crops_projected group by cell5m, iso, crop), t as (SELECT sum(prod) val, cell5m, crop, irrigation FROM crops_projected group by cell5m, crop, irrigation), tt as (SELECT sum(prod) val, cell5m, crop FROM crops_projected group by cell5m, crop), ff as (select * from r union all select *, 'all' as irrigation from s union all select *, 'WORLD' as iso from t union all select *, 'WORLD' as iso, 'all' as irrigation from tt) SELECT CDB_JenksBins(array_agg(val::numeric),5) as bucket, crop, irrigation, iso from ff where val is not null group by crop, irrigation, iso")

account = 'wri-rw'
urlCarto = f"https://{account}.carto.com/api/v2/sql"

sql = {"q": query}
r = requests.get(urlCarto, params=sql)

r.json()
API_KEY = getpass.getpass('RW Token_ID:')
# For API key authenticated requests, we need to create an APIKeyAuthClient instance:
USERNAME="wri-rw"
USR_BASE_URL = "https://{user}.carto.com/".format(user=USERNAME)
auth_client = APIKeyAuthClient(api_key=API_KEY, base_url=USR_BASE_URL)
def long_lasting_SQL_queries(sql_query, table_name):
    # For long lasting SQL queries we use the batch SQL API.
    LIST_OF_SQL_QUERIES = [sql_query]
    
    batchSQLClient = BatchSQLClient(auth_client)
    createJob = batchSQLClient.create(LIST_OF_SQL_QUERIES)
    
    # Check the status of a job with the job_id every 10 s
    readJob = batchSQLClient.read(createJob['job_id'])
    
    timeout = time.time() + 60*60  # 1 hour from now
    while readJob.get('status') != 'done':
        time.sleep(10)
        print(readJob.get('status'))
        if readJob.get('status') == 'failed':
            print('Job failed.')
            break
        if time.time() > timeout:
            readJob = batchSQLClient.read(createJob['job_id'])
            # Cancel a job given its job_id
            if readJob.get('status') != 'donne':
                cancelJob = batchSQLClient.cancel(createJob['job_id'])     
                print('Job cancelled after 1 hour running')
                break
            
        readJob = batchSQLClient.read(createJob['job_id'])
       
    # Read the table
    sql = SQLClient(auth_client)
    data = sql.send("select * from "+table_name)
    
    # Drop the table
    sql = SQLClient(auth_client)
    sql.send("DROP TABLE "+table_name)
        
    return data
table_name = 'job_result'

#"with r as (SELECT sum(prod) val, cell5m, iso, crop, irrigation FROM crops_projected group by cell5m, iso, crop, irrigation) SELECT CDB_JenksBins(array_agg(val::numeric),5) as bucket, crop, irrigation, iso from r where val is not null group by crop, irrigation, iso"

queries =[
    """
    with r as (SELECT sum(prod) val, cell5m, iso, crop, irrigation 
           FROM crops_projected group by cell5m, iso, crop, irrigation), 
     s as (SELECT sum(prod) val, cell5m, iso, crop 
           FROM crops_projected group by cell5m, iso, crop),
     t as (SELECT sum(prod) val, cell5m, crop, irrigation 
           FROM crops_projected group by cell5m, crop, irrigation), 
     tt as (SELECT sum(prod) val, cell5m, crop 
            FROM crops_projected group by cell5m, crop), 
     ff as (select val, cell5m, iso, crop, irrigation from r 
            union all 
            select val, cell5m, iso, crop, 'all' as irrigation from s 
            union all 
            select val, cell5m, 'WORLD' as iso, crop, irrigation  FROM t 
            union all 
            select val, cell5m, 'WORLD' as iso, crop, 'all' as irrigation from tt) 
SELECT CDB_JenksBins(array_agg(val::numeric),5) as bucket, crop, irrigation, iso 
from ff 
where val is not null 
group by crop, irrigation, iso
    """
]
query =("WITH a as ("+queries[0]+")"
        " SELECT * INTO "+table_name+" FROM a")  

data = long_lasting_SQL_queries(query, table_name)
df = pd.DataFrame(data.get("rows"))
df.info()
!ls /home/jovyan/work/datasets
df.to_csv(f'{baseFolder}/aqueduct/dst/spam2010/crop_buckets_v2.csv')<jupyter_output><empty_output><jupyter_text>```sqlselect array_agg(bucket::numeric) as bucket from (select json_array_elements(bucket::json)::text as bucket FROM "wri-rw".crop_buckets_v2 where iso='ISR' and irrigation='all' and crop = 'groundnut') s```<jupyter_code>import re
locations = "''asdfa''asdfasdf''asdfasdf'Ssdfasdf'''q"
myexpr= r"(?!'')((?<=[a-zA-Z0-9+-]|\s|[\,\.\\])'(?=[a-zA-Z0-9+-]|\s|[\,\.\\]))|(\\)|(\/)"
re.sub(myexpr,"",locations)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.11891697615909907, 'fraction_numerical': 0.01957589553132862, 'mean_word_length': 2.9954049396898332, 'pattern_counts': {'":': 43, '<': 90, '<?xml version=': 0, '>': 97, 'https://': 5, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 35, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '22', 'dir': 'jupyter-structured-clean-dedup', 'id': '247564', 'license': 'MIT', 'n_tokens_mistral': 13883, 'n_tokens_neox': 12953, 'n_words': 3078, 'path': 'Aqueduct/lab/AQ_1_V3_spam2010_data.ipynb', 'repo_name': 'resource-watch/notebooks', 'size': '41735'}