Jupyter notebook to add area to calculate country geojson

country_area.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ded1878e-3e4b-45ce-a73e-894ff9853562",
   "metadata": {},
   "outputs": [],
   "source": [
    "import json\n",
    "from shapely.geometry import Polygon, MultiPolygon, shape, mapping\n",
    "from shapely.affinity import translate\n",
    "from shapely.ops import transform\n",
    "import pyproj\n",
    "from functools import partial"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "aaffd500-d3af-4c92-afab-22c43b3d474c",
   "metadata": {},
   "outputs": [],
   "source": [
    "input_geojson = \"/path/to/countries.geojson\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e5c418e9-93e3-4e20-9fa9-f58e66f0f0de",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "with open(input_geojson, \"r\") as f:\n",
    "    contents = json.load(f)\n",
    "    \n",
    "features = contents[\"features\"]\n",
    "\n",
    "countries = []\n",
    "for feature in features:\n",
    "    mpoly = shape(feature[\"geometry\"])\n",
    "    props = feature[\"properties\"]\n",
    "    \n",
    "    new_geoms = []\n",
    "    for i, part in enumerate(mpoly):\n",
    "        if part.centroid.x < 0:\n",
    "            new_geoms.append(translate(part, xoff=360.0, yoff=0.0, zoff=0.0))\n",
    "        else:\n",
    "            new_geoms.append(part)\n",
    "    new_mpoly = MultiPolygon(new_geoms)\n",
    "    \n",
    "    geom_aea = transform(\n",
    "        partial(\n",
    "            pyproj.transform,\n",
    "            pyproj.Proj(init='EPSG:4326'),\n",
    "            pyproj.Proj(\n",
    "                proj='aea',\n",
    "                lat_1=new_mpoly.bounds[1],\n",
    "                lat_2=new_mpoly.bounds[3]\n",
    "            )\n",
    "        ),\n",
    "        new_mpoly\n",
    "    )\n",
    "    \n",
    "    props[\"area_deg\"] = new_mpoly.area\n",
    "    props[\"area_aea\"] = geom_aea.area\n",
    "    \n",
    "    countries.append({\n",
    "        \"type\": \"Feature\",\n",
    "        \"properties\": props,\n",
    "        \"geometry\": mapping(new_mpoly)\n",
    "    })\n",
    "    \n",
    "sorted_features = sorted(countries, key=lambda d: d[\"properties\"][\"area_aea\"])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "22e6f037-1d34-4072-bd2f-da15ff7d45cb",
   "metadata": {},
   "outputs": [],
   "source": [
    "output_geojson = \"/path/to/output/area_sorted_countries.geojson\"\n",
    "with open(output, 'w') as f:\n",
    "    json.dump(\n",
    "        {\n",
    "            \"type\": \"FeatureCollection\",\n",
    "            \"features\": sorted_features\n",
    "        }, f\n",
    "    )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "73982c19-4886-4c6a-b7d9-33baf2e6042c",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.8.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}