<jupyter_start><jupyter_text>Tutorial: Examining Protected Areas and Their Benthic and Geomorphic Distribution PurposeIn this tutorial, we'll walk through using the [World Database on Protected Areas](https://www.protectedplanet.net/en/thematic-areas/wdpa) (WDPA) data and the Allen Coral Atlas (ACA) to determine the amount of ACA mapped area that intersects protected areas in Tonga SetupWe recommend using a Conda environment with the packages listed below. Instructions to add a Conda environment to a jupyter notebook can be found [here](https://medium.com/@nrk25693/how-to-add-your-conda-environment-to-your-jupyter-notebook-in-just-4-steps-abeab8b8d084), and Conda documentation can be found [here](https://docs.conda.io/en/latest/).  Contents  1. Obtaining Data from WDPA 2. Visualizing Protected Areas 3. Determime What ACA Mapped Area Intersects Protected Areas 4. Visualizing Unprotected Areas  Considerations - for ipyLeaflet maps to display, you may need to locally run the following commands:         conda install -n your-environment-name -c conda-forge widgetsnbextension        conda install -n your-environment-name -c conda-forge ipywidgets    Packages to Import<jupyter_code>from ipyleaflet import (Map, GeoData, basemaps, WidgetControl, GeoJSON,
    LayersControl, Icon, Marker,basemap_to_tiles, Choropleth,
    MarkerCluster, Heatmap,SearchControl, 
    FullScreenControl)
from ipywidgets import Text, HTML
from branca.colormap import linear

import geopandas as gpd
import json
import requests

from osgeo import ogr

from owslib.wfs import WebFeatureService
from owslib.wms import WebMapService<jupyter_output><empty_output><jupyter_text>1. Obtaining data from WDPAWe'll be looking at Tonga again for protected areas. The data can be obtained from [here](https://www.protectedplanet.net/country/TON).We'll be using the shapefiles which is listed as 'SHP' on the WDPA site.When downloaded, the Tonga shapefiles have three parts and are located in the following folders: WDPA_WDOECM_TON_shp0, WDPA_WDOECM_TON_shp1, WDPA_WDOECM_TON_shp2.Using geopandas, we'll work with the shapefiles as geodataframes for easier data analysis.Let's start off by opening our downloaded shapefiles as geodataframes for all of the protected areas for Tonga(Note: file path may need to be changed to local location of downloaded Tonga shapefiles):<jupyter_code>tonga1 = gpd.read_file('data/WDPA_WDOECM_Feb2021_Public_TON_shp/WDPA_WDOECM_Feb2021_Public_TON_shp_0/WDPA_WDOECM_Feb2021_Public_TON_shp-polygons.shp')
tonga2 = gpd.read_file('data/WDPA_WDOECM_Feb2021_Public_TON_shp/WDPA_WDOECM_Feb2021_Public_TON_shp_1/WDPA_WDOECM_Feb2021_Public_TON_shp-polygons.shp')
tonga3 = gpd.read_file('data/WDPA_WDOECM_Feb2021_Public_TON_shp/WDPA_WDOECM_Feb2021_Public_TON_shp_2/WDPA_WDOECM_Feb2021_Public_TON_shp-polygons.shp')<jupyter_output><empty_output><jupyter_text>Here, we have combined all of the shapefiles into one dataframe, `tonga`.<jupyter_code>tonga_protected = tonga1.append(tonga2)
tonga_protected = tonga_protected.append(tonga3)<jupyter_output><empty_output><jupyter_text>Now let's take a look at the data! The main pieces of data we're interested in are the protected area names, designation, managing authority, and the geometry of the area.<jupyter_code>tonga_protected<jupyter_output><empty_output><jupyter_text>2. Visualizing Protected RegionsFor visualization, we'll be using IPyLeaflet. Further details on IPyLeaflet library can be found in its official documentation [here](https://ipyleaflet.readthedocs.io/en/latest/) and in this tutorial [here](https://towardsdatascience.com/ipyleaflet-interactive-mapping-in-jupyter-notebook-994f19611e79).We'll first create a basemap zoomed in with Tonga as the center.<jupyter_code>center = [-21.248,-175.188]
zoom = 10


tonga_map = Map(basemap=basemaps.OpenStreetMap.Mapnik, center=center, zoom=zoom)
tonga_map<jupyter_output><empty_output><jupyter_text>Next, we'll add a layer which displays all the protected areas from the geometries in our `tonga_protected` dataframe.<jupyter_code>protected_areas = GeoData(geo_dataframe = tonga_protected,
                  style={'color': 'black', 
                         'fillColor': '#E0D071',
                         'opacity':0.03, 
                         'weight':1.9, 'dashArray': '2', 
                         'fillOpacity':0.6},
                  hover_style={'fillColor': '#b08a3e', 
                               'fillOpacity': 0.8},
                  name = 'Tonga1')
tonga_map.add_layer(protected_areas)
tonga_map<jupyter_output><empty_output><jupyter_text>We'll now add a guide which will display the name of the protected region, designation, and managing authority when hovering over a particular protected area.<jupyter_code>html = HTML('''Hover Over Protected Areas''')
html.layout.margin = '0px 20px 20px 20px'
control = WidgetControl(widget=html, position='topright')
tonga_map.add_control(control)

def update_html(feature, **kwargs):
    html.value = '''
    <h3>{}</h3>
    <h4>Designation type: {}</h4>
    <h4>Managing authority: {}</h4>
    '''.format(feature['properties']['ORIG_NAME'],
              feature['properties']['DESIG'],
              feature['properties']['MANG_AUTH'])

protected_areas.on_hover(update_html)
tonga_map<jupyter_output><empty_output><jupyter_text>3. Determine What ACA Mapped Areas Are Protected The EEZ data is obtained from MarineRegions and the benthic and geomorphic cover data will be from the Allen Coral Atlas.We'll repeat steps from our tutorial that obtains benthic and geomorphic cover from a given EEZ:First, we'll create a dictionary with the keys being all the areas within the ACA and the values being the associated Marine Regions ID. We'll need the Marine Regions ID in the second section in order to obtain the EEZ of a specific country/area. To do so, we'll grab the Marine Regions EEZ map and restrict the output with a bounding box that is representative of the ACA area.Note that this request can take up to 3 minutes.<jupyter_code>wfs = WebFeatureService(url='https://geo.vliz.be/geoserver/MarineRegions/wfs', version='1.1.0')
response = wfs.getfeature(typename='MarineRegions:eez', 
                          bbox=(-180,-23.5,180,20), 
                          srsname='urn:x-ogc:def:crs:EPSG:4326', 
                          outputFormat='application/json')
eez_data = json.load(response)<jupyter_output><empty_output><jupyter_text>Here, we create a dictionary and ensure that disputed and/or shared territories are included with the following filtering parameters.<jupyter_code>mrgid_dict = {}
for country in eez_data['features']:
    if country['properties']['pol_type'] == '200NM':
        mrgid_dict[country['properties']['territory1']] = country['properties']['mrgid']
    else:
        mrgid_dict[country['properties']['geoname']] = country['properties']['mrgid']<jupyter_output><empty_output><jupyter_text>We'll also use our `get_eez_map` function from the first tutorial to get Tonga's EEZ again.<jupyter_code>def get_eez_map(country_mrgid):
    '''Given a country mrgid(Marine region country code), returns a geojson of a country\'s EEZ'''
    url = 'https://geo.vliz.be/geoserver/MarineRegions/wfs?service=WFS&version=1.1.0&request=GetFeature&typeName=eez&cql_filter=mrgid=' + str(country_mrgid) + '&outputformat=application/json'
    response = requests.get(url)
    eez = response.json()
    return eez
tonga_eez_json = get_eez_map(mrgid_dict['Tonga'])<jupyter_output><empty_output><jupyter_text>We'll now take the geojson of Tonga's EEZ and transform it to a GeoDataFrame.<jupyter_code>tonga_eez = gpd.GeoDataFrame.from_features(tonga_eez_json)<jupyter_output><empty_output><jupyter_text>Next we will (1) make a bounding box, (2) get the benthic and geomorphic layers, and (3) obtain the intersection between the ACA and EEZ layers.Here is a summary of the functions to achieve these steps:- `make_bounding_box(geojson)`: Transforms a geojson to a bounding box- `get_aca_layer(layer, bounding_box)`: Given a bounding box, a request is made and downloads either the benthic and geomorphic layers from the Allen Coral Atlas<jupyter_code>def make_bounding_box(geojson):
    '''From a geojson, creates a bounding box which is a list of 4 coordinates'''
    geom_text = json.dumps(geojson["features"][0]["geometry"])
    geom = ogr.CreateGeometryFromJson(geom_text)
    env = geom.GetEnvelope()

    MAX_LONG = 180.0
    MIN_LONG = -180.0
    MAX_LAT = 90.0
    MIN_LAT = -90.0

    return (max(min(MAX_LONG, env[0]), MIN_LONG),
            max(min(MAX_LAT, env[2]), MIN_LAT),
            max(min(MAX_LONG, env[1]), MIN_LONG),
            max(min(MAX_LAT, env[3]), MIN_LAT))
def get_aca_layer(layer, bounding_box):
    '''Gets either geomorphic or benthic layer in geojson from the ACA and returns is a geopandas GeoDataFrame'''
    # layer options: ['reef_polygons_benthic_expanded'] or ['reef_polygons_geomorphic_expanded']
    wms = WebMapService('https://allencoralatlas.org/geoserver/coral-atlas/wms', version='1.3.0')
    country_aca_layer = wms.getmap(layers=layer, 
                                   srs='EPSG:4326', 
                                   bbox=bounding_box, 
                                   size=(256,256),
                                   timeout=60,
                                   format='application/json;type=geojson')
    try:
        country_aca_layer_geojson = json.load(country_aca_layer)
        return gpd.GeoDataFrame.from_features(country_aca_layer_geojson)
    except ValueError:
        if country_aca_layer._response.status_code == 429:
            print('Too many requests have been made in a given amount of time, please wait another minute to send another request.')
        else:
            print('Error acquiring layer, response code: ' + country_aca_layer._response.status_code)<jupyter_output><empty_output><jupyter_text>Let's first get the bounding box.<jupyter_code>tonga_bbox = make_bounding_box(tonga_eez_json)<jupyter_output><empty_output><jupyter_text>We'll now take our `tongabbox` and obtain the benthic and geomorphic layers using our previously defined `get_aca_layer` function we defined earlier. Currently, to request the benthic layer, ['benthic_data_verbose'], is used and to request the geomorphic layer, ['geomorphic_data_verbose'], is used.<jupyter_code>benthic_tonga = get_aca_layer(['benthic_data_verbose'], tonga_bbox)
geomorphic_tonga = get_aca_layer(['geomorphic_data_verbose'], tonga_bbox)<jupyter_output><empty_output><jupyter_text>We'll now obtain the intersection between the ACA layers and the Marine Regions EEZ. The intersection will be a GeoDataFrame which will make data analysis easier.<jupyter_code>benthic_tonga_eez = gpd.overlay(benthic_tonga,
                                tonga_eez,
                                how="intersection").set_crs('EPSG:4326')
geomorphic_tonga_eez = gpd.overlay(geomorphic_tonga,
                                   tonga_eez,
                                   how="intersection").set_crs('EPSG:4326')<jupyter_output><empty_output><jupyter_text>We'll now determine the benthic and geomorphic cover within Tonga's current protected areas.<jupyter_code>benthic_tonga_eez_protected = gpd.overlay(benthic_tonga_eez,
                                          tonga_protected,
                                          how="intersection").set_crs('EPSG:4326')
geomorphic_tonga_eez_protected = gpd.overlay(geomorphic_tonga_eez,
                                             tonga_protected,
                                             how="intersection").set_crs('EPSG:4326')<jupyter_output><empty_output><jupyter_text>We can now use IPyLeaflet to visualize the ACA mapped area within the protected areas.<jupyter_code>benthic = GeoData(geo_dataframe = benthic_tonga_eez_protected,
                  style={'color': 'orange', 
                         'fillColor': '#b0553e',
                         'opacity':0.03, 
                         'weight':1.9, 'dashArray': '2', 
                         'fillOpacity':0.8},
                  hover_style={'fillColor': '#b0553e', 
                               'fillOpacity': 0.9},
                  name = 'benthic')
tonga_map.add_layer(benthic)
tonga_map<jupyter_output><empty_output><jupyter_text>Let's now calculate the numerical breakdown of the amount of benthic and geomorphic cover that is protected. We'll first reproject our data to the Mollweide projection which is an equal area projection. This will give a more accurate estimate of the protected areas in square meters.We'll also add a 'total_area' column to the dataframe to simplify a calculation.<jupyter_code>benthic_tonga_eez = benthic_tonga_eez.to_crs('ESRI:54009')
geomorphic_tonga_eez = geomorphic_tonga_eez.to_crs('ESRI:54009')
benthic_tonga_eez_protected = benthic_tonga_eez_protected.to_crs('ESRI:54009')
geomorphic_tonga_eez_protected = geomorphic_tonga_eez_protected.to_crs('ESRI:54009')<jupyter_output><empty_output><jupyter_text>Using Pandas, we'll now query the `benthic_tonga_eez_protected` and `geomorphic_tonga_eez_protected` dataframes to calculate Tonga's percent distribution of subtypes for benthic and geomorphic alongside the percentage that is protected.First we'll add a total_area column to simplify the calculations:<jupyter_code>benthic_tonga_eez['total_area'] = benthic_tonga_eez['geometry'].map(lambda p: p.area / 10**6)
geomorphic_tonga_eez['total_area'] = geomorphic_tonga_eez['geometry'].map(lambda p: p.area / 10**6)
benthic_tonga_eez_protected['total_area'] = benthic_tonga_eez_protected['geometry'].map(lambda p: p.area / 10**6)
geomorphic_tonga_eez_protected['total_area'] = geomorphic_tonga_eez_protected['geometry'].map(lambda p: p.area / 10**6)
benthic_tonga_dist = benthic_tonga_eez.groupby('class_name')['total_area'].sum().to_frame()
benthic_tonga_dist['protected_area'] = benthic_tonga_eez_protected.groupby('class_name')['total_area'].sum().to_frame()
benthic_tonga_dist['percent_protected'] = benthic_tonga_dist['protected_area']/benthic_tonga_dist['total_area'] * 100
benthic_tonga_dist<jupyter_output><empty_output><jupyter_text>Now let's do the same for the geomorphic data.<jupyter_code>geomorphic_tonga_dist = geomorphic_tonga_eez.groupby('class_name')['total_area'].sum().to_frame()
geomorphic_tonga_dist['protected_area'] = geomorphic_tonga_eez_protected.groupby('class_name')['total_area'].sum().to_frame()
geomorphic_tonga_dist['percent_protected'] = geomorphic_tonga_dist['protected_area']/geomorphic_tonga_dist['total_area'] * 100
geomorphic_tonga_dist<jupyter_output><empty_output><jupyter_text>4. Visualizing Unprotected Areas<jupyter_code>unprotected_benthic_tonga = gpd.overlay(benthic_tonga_eez.to_crs('EPSG:4326'), tonga_protected, how='difference')
center = [-21.248,-175.188]
zoom = 10

tonga_unprotected_areas_map = Map(basemap=basemaps.OpenStreetMap.Mapnik, center=center, zoom=zoom)
unprotected_benthic = GeoData(geo_dataframe = unprotected_benthic_tonga,
                  style={'color': 'orange', 
                         'fillColor': '#b0553e',
                         'opacity':0.03, 
                         'weight':1.9, 'dashArray': '2', 
                         'fillOpacity':0.8},
                  hover_style={'fillColor': '#b0553e', 
                               'fillOpacity': 0.9},
                  name = 'benthic')
tonga_unprotected_areas_map.add_layer(unprotected_benthic)
tonga_unprotected_areas_map<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.09493588913353193, 'fraction_numerical': 0.01910374303846652, 'mean_word_length': 4.2706484641638225, 'pattern_counts': {'":': 0, '<': 91, '<?xml version=': 0, '>': 91, 'https://': 9, 'lorem ipsum': 0, 'www.': 2, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 19, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '21', 'dir': 'jupyter-structured-clean-dedup', 'id': '546170', 'license': 'MIT', 'n_tokens_mistral': 5318, 'n_tokens_neox': 4971, 'n_words': 1151, 'path': 'WMS/protected_areas_tutorial.ipynb', 'repo_name': 'CoralMapping/tutorials', 'size': '15442'}