<jupyter_start><jupyter_text>Critical loads: Soil workflow (high-resolution method; 2018 onwards)In Spring 2018, the worflow for calculating critical loads was refined to make use of new, higher-resolution input datasets. During November 2018, data handling for the Critical Loads project was also redesigned, with the ultimate aim of centralising all key datasets on NIVA's new cloud platform. This notebook deals with data processing for the soil exceedance calculations, using the new 0.1 degree deposition grid and the "old" soil CLs. Exceedances are calculated as $$E_{soil} = S_{dep} − CL_{soil}$$(See e-mail from Kari received 01/11/2017 at 13.47 for details).**Note:** Deposition data supplied prior to 2017-18 use a different spatial grid (the "BLR" grid). **The workflow here only applies to data supplied from 2017-18 onwards**.   1. Organising soil data on the DSToolkitNIVA's JupyterHub includes a PostGIS database server capable of storing relational and vector geospatial datasets. I have created a database named `critical_loads` and, within this, a schema named `soil`. This schema contains the following table: * **s_critical_load:** Non-spatial table specifying the soil critical load for around 660 BLR cells with "productive forest" (see e-mail from Kari received 16.11.2017 at 12.52 for details). This dataset was originally exported from `RESA2.TALEGREN_VALUES`  2. Define data series and resolution of interestChoose the dataset you whish to work with (see notebook 01 for a list of available `series_ids`.**Note:** If you want to work with the 30 m resolution data, you will need to **log-in on one of the high memory machines**, as the grids involved are large.<jupyter_code># Series of interest
ser_id = 28
short_name = '1216'

# Choose cell size (30m or 60m)
cell_size = 120<jupyter_output><empty_output><jupyter_text>3. Create raster of soil critical loadsThe code below does not need to be re-run unless the soil critical loads are updated.<jupyter_code># Snap raster
snap_tif = f'/home/jovyan/shared/critical_loads/raster/sat_veg_{cell_size}m_norut.tif'

# Soil CL raster to create
cl_tif = f'/home/jovyan/shared/critical_loads/raster/soil_cl_mgspm2pyr_{cell_size}m.tif'
# Get all CLs for S
sql = ("SELECT * FROM soil.s_critical_load")
cl_df = pd.read_sql(sql, eng)

# Read BLR grid
gdf = nivapy.da.read_postgis('deposition', 'dep_grid_blr', eng)

# Join
gdf = gdf.merge(cl_df, on='cell_id')

# Rasterise
gdf.to_file('temp_soil_cl.geojson', driver='GeoJSON')
cl.vec_to_ras('temp_soil_cl.geojson', cl_tif, snap_tif, 'cl_mgSpm2', -1, 'Int16')
os.remove('temp_soil_cl.geojson')<jupyter_output><empty_output><jupyter_text>4. Calculate exceedances for soil<jupyter_code># Raster of S deposition (created in notebook 01)
sdep_tif = f'/home/jovyan/shared/critical_loads/raster/deposition/sdep_mgpm2pyr_{short_name}_{cell_size}m.tif'

# Exceedance .tif to create
ex_tif = f'/home/jovyan/shared/critical_loads/raster/exceedance/soil_ex_mgpm2pyr_{short_name}_{cell_size}m.tif'
# Read tifs
dep_s, ndv, epsg, extent = nivapy.spatial.read_raster(sdep_tif)
dep_s = dep_s.astype(np.float32)
dep_s[dep_s==ndv] = np.nan

cl_soil, ndv, epsg, extent = nivapy.spatial.read_raster(cl_tif)
cl_soil = cl_soil.astype(np.float32)
cl_soil[cl_soil==ndv] = np.nan

# Exceedance
ex_soil = dep_s - cl_soil
ex_soil[ex_soil<0] = 0

# Get exceeded area
ex_area = np.count_nonzero(ex_soil > 0)*cell_size*cell_size/1.E6
nor_area = np.count_nonzero(~np.isnan(ex_soil))*cell_size*cell_size/1.E6
ex_pct = 100*ex_area/nor_area

ex_df = pd.DataFrame({'series_id':ser_id,
                      'total_area_km2':nor_area,
                      'exceeded_area_km2':ex_area,
                      'exceeded_area_pct':ex_pct},
                     index=[0])
ex_df = ex_df.round(0).astype(int)
ex_df['medium'] = 'soil'
ex_df = ex_df[['series_id', 'medium', 'total_area_km2', 'exceeded_area_km2', 
               'exceeded_area_pct']]

# Set NaN to -1
ex_soil[np.isnan(ex_soil)] = -1

# Save geotiff
cl.write_geotiff(ex_soil, ex_tif, snap_tif, -1, gdal.GDT_Int16)

ex_df<jupyter_output><empty_output><jupyter_text>**If you connected to the database with `admin=True`**, these results can be saved back to the database.<jupyter_code>## Write summary data to db
#ex_df.to_sql('national_summary', 
#             eng,
#             'summaries',
#             if_exists='append',
#             index=False)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.10326086956521739, 'fraction_numerical': 0.024682971014492752, 'mean_word_length': 4.819499341238472, 'pattern_counts': {'":': 0, '<': 18, '<?xml version=': 0, '>': 18, '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': '4', 'dir': 'jupyter-structured-clean-dedup', 'id': '241498', 'license': 'MIT', 'n_tokens_mistral': 1677, 'n_tokens_neox': 1541, 'n_words': 437, 'path': 'notebooks/05_soil_new_grid.ipynb', 'repo_name': 'JamesSample/critical_loads_2', 'size': '4416'}