WDPA_WDOECM_Pre_Processing_Python_Script_V2.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "## WDPA and WD-OECM pre-processing script - public and restricted data (Python)\n",
    "\n",
    "# Import necessary Libraries\n",
    "from dateutil.relativedelta import relativedelta\n",
    "import datetime\n",
    "from arcgis.gis import GIS\n",
    "import arcpy\n",
    "import math\n",
    "import os\n",
    "from arcgis.features import FeatureLayer\n",
    "from arcgis.features import FeatureLayerCollection\n",
    "import pandas as pd\n",
    "\n",
    "#Itentify the month and year corresponding to the version of the WDPA and WD-OECM used in the Protected Planet Report 2024 (August 2024)\n",
    "date = datetime.datetime(2024, 8, 1)\n",
    "month = date.strftime(\"%b\")  # \"Aug\"\n",
    "Month = date.strftime(\"%B\")   # \"August\"\n",
    "year = date.strftime(\"%Y\")    # \"2024\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Set default file geodatabase for intermediate outputs\n",
    "arcpy.env.workspace = r\"path_to_file_geodatabase_for_intermediate_outputs\"\n",
    "\n",
    "#Define paths to WDPA and WDOECM points and polygons \n",
    "\n",
    "# Public WDPA paths\n",
    "wdpa_path = r'path_to_all_public_wdpa_point_polygon_file_geodatabase'\n",
    "wdpa_poly_path = os.path.join(wdpa_path,'WDPA_poly_{0}{1}'.format(month,year))\n",
    "wdpa_point_path = os.path.join(wdpa_path,'WDPA_point_{0}{1}'.format(month,year))\n",
    "print(wdpa_path)\n",
    "\n",
    "# Public WDOECM paths\n",
    "wdoecm_path = r'path_to_all_public_wdoecm_point_polygon_file_geodatabase'\n",
    "wdoecm_poly_path = os.path.join(wdoecm_path,'WDOECM_poly_{0}{1}'.format(month,year)) \n",
    "wdoecm_point_path = os.path.join(wdoecm_path,'WDOECM_point_{0}{1}'.format(month,year))\n",
    "\n",
    "#Restricted WDPA paths \n",
    "restricted_wdpa_path = r'path_to_all_restricted_wdpa_point_polygon_file_geodatabase'\n",
    "restricted_wdpa_poly_path = os.path.join(restricted_wdpa_path,'wdpa_poly_restricted')\n",
    "restricted_wdpa_point_path = os.path.join(restricted_wdpa_path,'wdpa_point_restricted')\n",
    "\n",
    "Datasets = [\"wdpa\", \"wdoecm\"]\n",
    "\n",
    "#Read in and subset WDPA and WDOECM\n",
    "\n",
    "# Public WDPA points\n",
    "## Export features where REP_AREA isn't equaly to \"0\", DESIG_ENG is not equal to \"UNESCO-MAB Biosphere Reserve\", and STATUS is not equal to either \"Not Reported\" or \"Proposed\"\n",
    "arcpy.conversion.ExportFeatures(wdpa_point_path,r\"memory\\wdpa_point_mp\", \"REP_AREA <> 0 And DESIG_ENG <> 'UNESCO-MAB Biosphere Reserve' And STATUS <> 'Not Reported' And STATUS <> 'Proposed'\")\n",
    "## Convert to single part\n",
    "arcpy.management.MultipartToSinglepart(\"wdpa_point_mp\", \"wdpa_point\")\n",
    "\n",
    "# Public WDPA polygons\n",
    "## Export features where DESIG_ENG is not equal to \"UNESCO-MAB Biosphere Reserve\", and STATUS is not equal to either \"Not Reported\" or \"Proposed\"\n",
    "arcpy.conversion.ExportFeatures(wdpa_poly_path,r\"memory\\wdpa_poly_mp\", \"DESIG_ENG <> 'UNESCO-MAB Biosphere Reserve' And STATUS <> 'Not Reported' And STATUS <> 'Proposed'\")\n",
    "## Convert to single part\n",
    "arcpy.management.MultipartToSinglepart(\"wdpa_poly_mp\", \"wdpa_poly\")\n",
    "\n",
    "# Restricted WDPA points\n",
    "## Export features where REP_AREA isn't equaly to \"0\", DESIG_ENG is not equal to \"UNESCO-MAB Biosphere Reserve\", and STATUS is not equal to either \"Not Reported\" or \"Proposed\"\n",
    "arcpy.conversion.ExportFeatures(restricted_wdpa_point_path,r\"memory\\restricted_wdpa_point_mp\", \"REP_AREA <> 0 And DESIG_ENG <> 'UNESCO-MAB Biosphere Reserve' And STATUS <> 'Not Reported' And STATUS <> 'Proposed'\")\n",
    "## Convert to single part\n",
    "arcpy.management.MultipartToSinglepart(\"restricted_wdpa_point_mp\", \"restricted_wdpa_point\")\n",
    "\n",
    "# Restricted WDPA polygons\n",
    "## Export features where DESIG_ENG is not equal to \"UNESCO-MAB Biosphere Reserve\", and STATUS is not equal to either \"Not Reported\" or \"Proposed\"\n",
    "arcpy.conversion.ExportFeatures(restricted_wdpa_poly_path,r\"memory\\restricted_wdpa_poly_mp\", \"DESIG_ENG <> 'UNESCO-MAB Biosphere Reserve' And STATUS <> 'Not Reported' And STATUS <> 'Proposed'\")\n",
    "## Convert to single part\n",
    "arcpy.management.MultipartToSinglepart(\"restricted_wdpa_poly_mp\", \"restricted_wdpa_poly\")\n",
    "\n",
    "# Public WDOECM points\n",
    "## Export features where REP_AREA is not equal to \"0\" and STATUS is not equal to either \"Not Reported\" or \"Proposed\"\n",
    "arcpy.conversion.ExportFeatures(wdoecm_point_path,r\"memory\\wdoecm_point_mp\", \"\")\n",
    "## Convert to single part\n",
    "arcpy.management.MultipartToSinglepart(\"wdoecm_point_mp\", \"wdoecm_point\")\n",
    "\n",
    "## Export features where STATUS is not equal to either \"Not Reported\" or \"Proposed\"\n",
    "arcpy.conversion.ExportFeatures(wdoecm_poly_path,r\"memory\\wdoecm_poly_mp\", \"\")\n",
    "## Convert to single part\n",
    "arcpy.management.MultipartToSinglepart(\"wdoecm_poly_mp\", \"wdoecm_poly\")\n",
    "\n",
    "##Buffer WDPA, WDOECM, and restricted WDPA point shapefiles to REP_AREA value; converting to polygon shapefiles. \n",
    "Datasets = [\"wdpa\", \"wdoecm\", \"restricted_wdpa\"]\n",
    "\n",
    "for i in Datasets:\n",
    "    point_fc = f\"{i}_point\"\n",
    "    \n",
    "    # Add buffer_radius field\n",
    "    arcpy.management.AddField(point_fc, \"buffer_radius\", \"LONG\", field_length=None, field_alias=\"\", field_is_nullable=\"NULLABLE\", field_is_required=\"NON_REQUIRED\")\n",
    "    \n",
    "    # Calculate the buffer radius based on the REP_AREA field value\n",
    "    arcpy.management.CalculateField(point_fc, \"buffer_radius\", \"math.sqrt(!REP_AREA!/3.1415)\", \"PYTHON3\", code_block=\"\")\n",
    "    \n",
    "    # Create a buffer field combining the radius with a unit string\n",
    "    arcpy.management.CalculateField(point_fc, \"buffer_field\", 'str(!buffer_radius!) + \" Kilometers\"', \"PYTHON3\", code_block=\"\")\n",
    "    \n",
    "    # Perform the buffering analysis\n",
    "    buffer_fc = f\"{i}_point_buffered\"\n",
    "    arcpy.analysis.PairwiseBuffer(point_fc, buffer_fc, \"buffer_field\", \"None\", None, \"GEODESIC\", \"0 Meters\")\n",
    "\n",
    "## Repair geometries of buffered points layers.\n",
    "arcpy.management.RepairGeometry(\"wdpa_point_buffered\", \"KEEP_NULL\", \"OGC\")\n",
    "arcpy.management.RepairGeometry(\"wdoecm_point_buffered\", \"KEEP_NULL\", \"OGC\")\n",
    "arcpy.management.RepairGeometry(\"restricted_wdpa_point_buffered\", \"KEEP_NULL\", \"OGC\")\n",
    "\n",
    "## Merge WDPA points + polys (all), WDOECM points + polys (all), and erase protected area extent from OECM extent to produce a\n",
    "# layer show additional coverage provided by OECMs only. \n",
    "arcpy.management.Merge([\"wdpa_point_buffered\", \"wdpa_poly\", \"restricted_wdpa_point_buffered\", \"restricted_wdpa_poly\"], \"wdpa_merged_public_restricted\")\n",
    "arcpy.management.Merge([\"wdoecm_point_buffered\", \"wdoecm_poly\"], \"wdoecm_merged\")\n",
    "arcpy.analysis.PairwiseErase(\"wdoecm_merged\", \"wdpa_merged_public_restricted\", \"wdoecm_merged_pa_erase\", \"\")\n",
    "\n",
    "## Repair geometries of newly created layers\n",
    "arcpy.management.RepairGeometry(\"wdpa_merged_public_restricted\", \"KEEP_NULL\", \"OGC\")\n",
    "arcpy.management.RepairGeometry(\"wdoecm_merged\", \"KEEP_NULL\", \"OGC\")\n",
    "arcpy.management.RepairGeometry(\"wdoecm_merged_pa_erase\", \"KEEP_NULL\", \"OGC\")\n",
    "\n",
    "# Write intermediate WDPA and WDOECM to local file geodatabases before clearing cache to free memory. \n",
    "#WDPA\n",
    "##Set working environment for output files\n",
    "wdpa_temp_path = r'path_to_file_geodatabase_for_intermediate_outputs'\n",
    "wdpa_temp_poly_point_path = os.path.join(wdpa_temp_path,\"wdpa_merged_public_restricted\")\n",
    "wdpa_output_path = r'path_to_file_geodatabase_for_preprocessing_outputs'\n",
    "wdpa_output_poly_point_path = os.path.join(wdpa_output_path,\"wdpa_poly_buffpoint_{0}{1}\".format(month, year))\n",
    "arcpy.conversion.ExportFeatures(wdpa_temp_poly_point_path, wdpa_output_poly_point_path)\n",
    "\n",
    "#WDOECM\n",
    "wdoecm_temp_path = r'path_to_file_geodatabase_for_intermediate_outputs'\n",
    "wdoecm_temp_poly_point_path = os.path.join(wdoecm_temp_path,\"wdoecm_merged\")\n",
    "wdoecm_output_path = r'path_to_file_geodatabase_for_preprocessing_outputs'\n",
    "wdoecm_output_poly_point_path = os.path.join(wdoecm_output_path,\"wdoecm_poly_buffpoint_{0}{1}\".format(month, year))\n",
    "arcpy.conversion.ExportFeatures(wdoecm_temp_poly_point_path, wdoecm_output_poly_point_path)\n",
    "\n",
    "wdoecm_pa_erase_temp_path = r'path_to_file_geodatabase_for_intermediate_outputs'\n",
    "wdoecm_pa_erase_temp_poly_point_path = os.path.join(wdoecm_pa_erase_temp_path,\"wdoecm_merged_pa_erase\")\n",
    "wdoecm_pa_erase_output_path = r'path_to_file_geodatabase_for_preprocessing_outputs'\n",
    "wdoecm_pa_erase_output_poly_point_path = os.path.join(wdoecm_pa_erase_output_path,\"wdoecm_pa_erase_poly_buffpoint_{0}{1}\".format(month, year))\n",
    "arcpy.conversion.ExportFeatures(wdoecm_pa_erase_temp_poly_point_path, wdoecm_pa_erase_output_poly_point_path)\n",
    "\n",
    "# Define the path to the geodatabase and the feature class\n",
    "basemap_path = r'path_to_file_geodatabase_containing_basemap'\n",
    "basemap_name = 'basemap_name'\n",
    "basemap_path_name = os.path.join(basemap_path, basemap_name)\n",
    "basemap = arcpy.management.MakeFeatureLayer(basemap_path_name, basemap_name)\n",
    "\n",
    "#Pairwise intersect WDPA and WD-OECM in memory with basemap\n",
    "##WDPA\n",
    "#Intersect (pairwise) WDPA (polygon + point data) with global basemap\n",
    "arcpy.analysis.PairwiseIntersect([wdpa_output_poly_point_path, basemap], \"wdpa_poly_point_basemap_intersect\", \"ALL\", \"\", \"INPUT\")\n",
    "#Repair interesected WDPA\n",
    "arcpy.management.RepairGeometry(\"wdpa_poly_point_basemap_intersect\", \"KEEP_NULL\", \"OGC\")\n",
    "#Export intersected WDPA from memory to output file geodatabase\n",
    "arcpy.conversion.FeatureClassToFeatureClass(\"wdpa_poly_point_basemap_intersect\", wdpa_output_path, \"wdpa_poly_point_public_restricted_basemap_intersect_{0}{1}\".format(month, year))\n",
    "\n",
    "##WDOECM (verion with PA extent erased)\n",
    "#Intersect (pairwise) WDOECM (polygon + point data) with global basemap\n",
    "arcpy.analysis.PairwiseIntersect([wdoecm_pa_erase_output_poly_point_path, basemap], \"wdoecm_basemap_intersect\", \"ALL\", \"\", \"INPUT\")\n",
    "#Repair intesected WDOECM\n",
    "arcpy.management.RepairGeometry(\"wdoecm_basemap_intersect\", \"KEEP_NULL\", \"OGC\")\n",
    "#Export intersected WDOECM from memory to output file geodatabase\n",
    "arcpy.conversion.FeatureClassToFeatureClass(\"wdoecm_basemap_intersect\", wdoecm_output_path, \"wdoecm_basemap_intersect_{0}{1}\".format(month, year))\n",
    "\n",
    "#Split WDPA and WDOECM by STATUS_YR field to assess change in coverage over time. \n",
    "wdpa_public_restricted_basemap_intersect_output_path = os.path.join(wdpa_output_path,\"wdpa_poly_point_public_restricted_basemap_intersect_{0}{1}\".format(month, year))\n",
    "\n",
    "#Define output paths for WDPA and WDOECM subset by STATUS_YR\n",
    "wdpa_status_yr_2020_and_less_path = os.path.join(wdpa_output_path,\"wdpa_poly_point_public_restricted_basemap_intersect_statusyr_2020_and_less_{0}{1}\".format(month, year))\n",
    "wdpa_status_yr_2022_and_less_path = os.path.join(wdpa_output_path,\"wdpa_poly_point_public_restricted_basemap_intersect_statusyr_2022_and_less_{0}{1}\".format(month, year))\n",
    "wdpa_status_yr_2024_and_less_path = os.path.join(wdpa_output_path,\"wdpa_poly_point_public_restricted_basemap_intersect_statusyr_2024_and_less_{0}{1}\".format(month, year))\n",
    "wdoecm_status_yr_2020_and_less_path = os.path.join(wdoecm_output_path,\"wdoecm_poly_point_basemap_intersect_statusyr_2020_and_less_{0}{1}\".format(month, year))\n",
    "wdoecm_status_yr_2022_and_less_path = os.path.join(wdoecm_output_path,\"wdoecm_poly_point_basemap_intersect_statusyr_2022_and_less_{0}{1}\".format(month, year))\n",
    "wdoecm_status_yr_2024_and_less_path = os.path.join(wdoecm_output_path,\"wdoecm_poly_point_basemap_intersect_statusyr_2024_and_less_{0}{1}\".format(month, year))\n",
    "\n",
    "#Subset pre-processed WDPA and WDOECM by STATUS_YR\n",
    "arcpy.conversion.ExportFeatures(wdpa_public_restricted_basemap_intersect_output_path, wdpa_status_yr_2020_and_less_path, \"STATUS_YR <= 2020\")\n",
    "arcpy.conversion.ExportFeatures(wdpa_public_restricted_basemap_intersect_output_path, wdpa_status_yr_2022_and_less_path, \"STATUS_YR <= 2022\")\n",
    "arcpy.conversion.ExportFeatures(wdpa_public_restricted_basemap_intersect_output_path, wdpa_status_yr_2024_and_less_path, \"STATUS_YR <= 2024\")\n",
    "arcpy.conversion.ExportFeatures(wdoecm_basemap_intersect_output_path, wdoecm_status_yr_2020_and_less_path, \"STATUS_YR <= 2020\")\n",
    "arcpy.conversion.ExportFeatures(wdoecm_basemap_intersect_output_path, wdoecm_status_yr_2022_and_less_path, \"STATUS_YR <= 2022\")\n",
    "arcpy.conversion.ExportFeatures(wdoecm_basemap_intersect_output_path, wdoecm_status_yr_2024_and_less_path, \"STATUS_YR <= 2024\")\n",
    "\n",
    "#Define fields to flatten (pairwise dissolve) WDPA and WDOECM (subset by different STATUS_YR values)\n",
    "## PA_DEF = WDPA/WDOECM field, binary value (0 = OECM, 1 = protected area)\n",
    "## Sub_region_NAME = Name of UNSD M49 sub-region in global basemap\n",
    "## ISO3 = WDPA/WDOECM field, containing International Standards Organization (ISO) 3 digit code corresponding to a country or territory in the WDPA/WDOECM.\n",
    "\n",
    "flattenGroups = [\"PA_DEF\", \"Sub_region_Name\", \"ISO3\"]\n",
    "\n",
    "#Create flattened WDPA for STATUS_YR <= 2024 for each field name in flattenGroups\n",
    "for i in flattenGroups:\n",
    "    output_path = f\"{wdpa_status_yr_2024_and_less_path}_{i}\"\n",
    "    arcpy.analysis.PairwiseDissolve(wdpa_status_yr_2024_and_less_path, output_path, i,None, \"SINGLE_PART\")\n",
    "    arcpy.management.RepairGeometry(output_path, \"DELETE_NULL\", \"OGC\") \n",
    "\n",
    "#Create flattened WDPA for STATUS_YR <= 2022 for each field name in flattenGroups\n",
    "for i in flattenGroups:\n",
    "    output_path = f\"{wdpa_status_yr_2022_and_less_path}_{i}\"\n",
    "    arcpy.analysis.PairwiseDissolve(wdpa_status_yr_2022_and_less_path, output_path, i,None, \"SINGLE_PART\")\n",
    "    arcpy.management.RepairGeometry(output_path, \"DELETE_NULL\", \"OGC\")\n",
    "\n",
    "#Create flattened WDPA for STATUS_YR <= 2020 for each field name in flattenGroups\n",
    "for i in flattenGroups:\n",
    "    output_path = f\"{wdpa_status_yr_2020_and_less_path}_{i}\"\n",
    "    arcpy.analysis.PairwiseDissolve(wdpa_status_yr_2020_and_less_path, output_path, i,None, \"SINGLE_PART\")\n",
    "    arcpy.management.RepairGeometry(output_path, \"DELETE_NULL\", \"OGC\")\n",
    "\n",
    "#Create flattened WDOECM for STATUS_YR <= 2024 for each field name in flattenGroups    \n",
    "for i in flattenGroups:\n",
    "    output_path = f\"{wdoecm_status_yr_2024_and_less_path}_{i}\"\n",
    "    arcpy.analysis.PairwiseDissolve(wdoecm_status_yr_2024_and_less_path, output_path, i,None, \"SINGLE_PART\")\n",
    "    arcpy.management.RepairGeometry(output_path, \"DELETE_NULL\", \"OGC\")\n",
    "    \n",
    "#Create flattened WDOECM for STATUS_YR <= 2022 for each field name in flattenGroups\n",
    "for i in flattenGroups:\n",
    "    output_path = f\"{wdoecm_status_yr_2022_and_less_path}_{i}\"\n",
    "    arcpy.analysis.PairwiseDissolve(wdoecm_status_yr_2022_and_less_path, output_path, i,None, \"SINGLE_PART\")\n",
    "    arcpy.management.RepairGeometry(output_path, \"DELETE_NULL\", \"OGC\")\n",
    "\n",
    "#Create flattened WDOECM for STATUS_YR <= 2020 for each field name in flattenGroups\n",
    "for i in flattenGroups:\n",
    "    output_path = f\"{wdoecm_status_yr_2020_and_less_path}_{i}\"\n",
    "    arcpy.analysis.PairwiseDissolve(wdoecm_status_yr_2020_and_less_path, output_path, i,None, \"SINGLE_PART\")\n",
    "    arcpy.management.RepairGeometry(output_path, \"DELETE_NULL\", \"OGC\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Buffer Points "
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.9.16"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}