WorldPopMap.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "9c4a0c23",
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import matplotlib as mpl\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import matplotlib.ticker as ticker\n",
    "from matplotlib.colors import ListedColormap\n",
    "from tqdm import tqdm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "7397f02d-f008-4ce8-94c0-9840688adcf2",
   "metadata": {},
   "outputs": [],
   "source": [
    "import rasterio\n",
    "from rasterio.transform import Affine\n",
    "from torchvision.transforms import ToTensor\n",
    "from sklearn.neighbors import KernelDensity\n",
    "# 24sec"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "e1f7f1cd-68d7-4046-9869-a8c6b40a0728",
   "metadata": {},
   "outputs": [],
   "source": [
    "file_name = 'WorldPop2015/WorldPop2015.tif'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "98d37379-4912-4cc4-9c5a-8725d69ce02b",
   "metadata": {},
   "outputs": [],
   "source": [
    "with rasterio.open(file_name) as image:\n",
    "    image_array = image.read()\n",
    "# 30 sec\n",
    "# 50 sec "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "c815709e-b395-4ff3-bd55-565dce55a884",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "torch.Size([43202, 1, 21384])\n"
     ]
    }
   ],
   "source": [
    "torch_image = ToTensor()(image_array)\n",
    "print(torch_image.shape)\n",
    "# for night light map : [43200, 1, 17406]\n",
    "# for population map : [43202, 1, 21384]\n",
    "# 60 sec\n",
    "# 1 min 15sec\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "db6a6293-65d7-4326-95dc-535f12de1ecc",
   "metadata": {},
   "outputs": [],
   "source": [
    "# ------------------------------------imag is defined here!!!------------------------------------\n",
    "imag = torch_image.permute(1,2,0).numpy()[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "4294821c-d5de-4475-ae76-f89302eefced",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(21384, 43202)"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "imag.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "b3100f5f",
   "metadata": {},
   "outputs": [],
   "source": [
    "def normalise(imag):\n",
    "    exclude_zero = imag != 0.0\n",
    "    imag = imag[exclude_zero]\n",
    "    min = np.min(imag)\n",
    "    max = np.max(imag)\n",
    "    normalised_imag = (imag - min) / (max - min)\n",
    "    return normalised_imag\n",
    "    \n",
    "normalised_imag = normalise(imag)\n",
    "# 6 sec\n",
    "# 40 sec\n",
    "# 58 sec"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "354ed82c",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.0\n",
      "1.0\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiMAAAGsCAYAAAAPJKchAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/H5lhTAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAk50lEQVR4nO3da3BU9f3H8c8SyAY1u4I0F8JyUeQmV4PARi3axqY0paQPLEMdkzpAq4YOmqolakVFu5laWpxKQaSYtkpjsQIdQGgaGxkkVoFkGkBpKZegZoNOdTdEXTB7/g86rP9tLuQsSX7u8n7NnAc5+Z2c7x4Z9z0nJ4nDsixLAAAAhvQxPQAAALiwESMAAMAoYgQAABhFjAAAAKOIEQAAYBQxAgAAjCJGAACAUcQIAAAwihgBAABGESMAAMCouIqRnTt3avbs2Ro8eLAcDoc2bdpk6/iHH35YDoejzXbxxRf3zMAAAOCc4ipGWlpaNGnSJK1cuTKm4++55x41NjZGbePGjdPNN9/czZMCAICuiqsYmTVrlh577DF9+9vfbvfzoVBI99xzj7KysnTxxRdr+vTpqq6ujnz+kksuUUZGRmRramrSwYMHNX/+/F56BQAA4H/FVYycy6JFi1RTU6OKigr94x//0M0336yvf/3r+te//tXu+rVr12rUqFG6/vrre3lSAABwVsLESENDg5599llt2LBB119/va644grdc889uu666/Tss8+2Wf/pp5/q+eef564IAACG9TU9QHepr69Xa2urRo0aFbU/FArpsssua7N+48aNam5uVlFRUW+NCAAA2pEwMXLq1CklJSVp7969SkpKivrcJZdc0mb92rVr9c1vflPp6em9NSIAAGhHwsTIlClT1NraqpMnT57zGZCjR4/qb3/7m/785z/30nQAAKAjcRUjp06d0uHDhyMfHz16VHV1dRo4cKBGjRqlW265RYWFhVq+fLmmTJmi999/X1VVVZo4caLy8/Mjx61bt06ZmZmaNWuWiZcBAAD+H4dlWZbpIbqqurpaN954Y5v9RUVFKi8v15kzZ/TYY4/pd7/7nd59910NGjRIM2bM0COPPKIJEyZIksLhsIYNG6bCwkI9/vjjvf0SAADA/4irGAEAAIknYX60FwAAxCdiBAAAGBUXD7CGw2G99957Sk1NlcPhMD0OAADoAsuy1NzcrMGDB6tPn47vf8RFjLz33nvyeDymxwAAADE4ceKEhgwZ0uHn4yJGUlNTJf33xbhcLsPTAACArggGg/J4PJH38Y7ERYyc/daMy+UiRgAAiDPnesSCB1gBAIBRxAgAADCKGAEAAEYRIwAAwChiBAAAGEWMAAAAo4gRAABgFDECAACMIkYAAIBRxAgAADCKGAEAAEYRIwAAwChiBAAAGEWMAAAAo/qaHsC04Uu2ttl3rCzfwCQAAFyYuDMCAACMIkYAAIBRxAgAADCKGAEAAEYRIwAAwChiBAAAGEWMAAAAo4gRAABgFDECAACMIkYAAIBRxAgAADDqvGKkrKxMDodDd911V6frNmzYoDFjxiglJUUTJkzQtm3bzue0AAAggcQcI2+++aaefvppTZw4sdN1u3fv1rx58zR//nzV1taqoKBABQUF2r9/f6ynBgAACSSmGDl16pRuueUWPfPMMxowYECna5988kl9/etf17333quxY8dq2bJluvrqq/XUU0/FNDAAAEgsMcVIcXGx8vPzlZube861NTU1bdbl5eWppqamw2NCoZCCwWDUBgAAElNfuwdUVFRo3759evPNN7u03u/3Kz09PWpfenq6/H5/h8f4fD498sgjdkcDAABxyNadkRMnTmjx4sV6/vnnlZKS0lMzqbS0VIFAILKdOHGix84FAADMsnVnZO/evTp58qSuvvrqyL7W1lbt3LlTTz31lEKhkJKSkqKOycjIUFNTU9S+pqYmZWRkdHgep9Mpp9NpZzQAABCnbN0Z+epXv6r6+nrV1dVFtqlTp+qWW25RXV1dmxCRJK/Xq6qqqqh9lZWV8nq95zc5AABICLbujKSmpmr8+PFR+y6++GJddtllkf2FhYXKysqSz+eTJC1evFgzZ87U8uXLlZ+fr4qKCu3Zs0dr1qzpppcAAADiWbf/BtaGhgY1NjZGPs7JydH69eu1Zs0aTZo0SS+++KI2bdrUJmoAAMCFyWFZlmV6iHMJBoNyu90KBAJyuVzd+rWHL9naZt+xsvxuPQcAABeirr5/87dpAACAUcQIAAAwihgBAABGESMAAMAoYgQAABhFjAAAAKOIEQAAYBQxAgAAjCJGAACAUcQIAAAwihgBAABGESMAAMAoYgQAABhFjAAAAKOIEQAAYBQxAgAAjCJGAACAUcQIAAAwihgBAABGESMAAMAoYgQAABhFjAAAAKOIEQAAYBQxAgAAjCJGAACAUcQIAAAwihgBAABGESMAAMAoYgQAABhFjAAAAKOIEQAAYJStGFm1apUmTpwol8sll8slr9erl19+ucP15eXlcjgcUVtKSsp5Dw0AABJHXzuLhwwZorKyMl155ZWyLEu//e1vNWfOHNXW1uqqq65q9xiXy6VDhw5FPnY4HOc3MQAASCi2YmT27NlRHz/++ONatWqVXn/99Q5jxOFwKCMjI/YJAQBAQov5mZHW1lZVVFSopaVFXq+3w3WnTp3SsGHD5PF4NGfOHB04cOCcXzsUCikYDEZtAAAgMdmOkfr6el1yySVyOp26/fbbtXHjRo0bN67dtaNHj9a6deu0efNmPffccwqHw8rJydE777zT6Tl8Pp/cbndk83g8dscEAABxwmFZlmXngNOnT6uhoUGBQEAvvvii1q5dq1dffbXDIPn/zpw5o7Fjx2revHlatmxZh+tCoZBCoVDk42AwKI/Ho0AgIJfLZWfccxq+ZGubfcfK8rv1HAAAXIiCwaDcbvc5379tPTMiScnJyRo5cqQkKTs7W2+++aaefPJJPf300+c8tl+/fpoyZYoOHz7c6Tqn0ymn02l3NAAAEIfO+/eMhMPhqLsYnWltbVV9fb0yMzPP97QAACBB2LozUlpaqlmzZmno0KFqbm7W+vXrVV1drR07dkiSCgsLlZWVJZ/PJ0l69NFHNWPGDI0cOVIfffSRnnjiCR0/flwLFizo/lcCAADikq0YOXnypAoLC9XY2Ci3262JEydqx44duummmyRJDQ0N6tPn85stH374oRYuXCi/368BAwYoOztbu3fv7tLzJQAA4MJg+wFWE7r6AEwseIAVAICe0dX3b/42DQAAMIoYAQAARhEjAADAKGIEAAAYRYwAAACjiBEAAGAUMQIAAIwiRgAAgFHECAAAMIoYAQAARhEjAADAKGIEAAAYRYwAAACjiBEAAGAUMQIAAIwiRgAAgFHECAAAMIoYAQAARhEjAADAKGIEAAAYRYwAAACjiBEAAGAUMQIAAIwiRgAAgFHECAAAMIoYAQAARhEjAADAKGIEAAAYRYwAAACjiBEAAGAUMQIAAIyyFSOrVq3SxIkT5XK55HK55PV69fLLL3d6zIYNGzRmzBilpKRowoQJ2rZt23kNDAAAEoutGBkyZIjKysq0d+9e7dmzR1/5ylc0Z84cHThwoN31u3fv1rx58zR//nzV1taqoKBABQUF2r9/f7cMDwAA4p/DsizrfL7AwIED9cQTT2j+/PltPjd37ly1tLRoy5YtkX0zZszQ5MmTtXr16i6fIxgMyu12KxAIyOVync+4bQxfsrXNvmNl+d16DgAALkRdff+O+ZmR1tZWVVRUqKWlRV6vt901NTU1ys3NjdqXl5enmpqaTr92KBRSMBiM2gAAQGKyHSP19fW65JJL5HQ6dfvtt2vjxo0aN25cu2v9fr/S09Oj9qWnp8vv93d6Dp/PJ7fbHdk8Ho/dMQEAQJywHSOjR49WXV2d/v73v+uOO+5QUVGRDh482K1DlZaWKhAIRLYTJ05069cHAABfHH3tHpCcnKyRI0dKkrKzs/Xmm2/qySef1NNPP91mbUZGhpqamqL2NTU1KSMjo9NzOJ1OOZ1Ou6MBAIA4dN6/ZyQcDisUCrX7Oa/Xq6qqqqh9lZWVHT5jAgAALjy27oyUlpZq1qxZGjp0qJqbm7V+/XpVV1drx44dkqTCwkJlZWXJ5/NJkhYvXqyZM2dq+fLlys/PV0VFhfbs2aM1a9Z0/ysBAABxyVaMnDx5UoWFhWpsbJTb7dbEiRO1Y8cO3XTTTZKkhoYG9enz+c2WnJwcrV+/Xg8++KDuv/9+XXnlldq0aZPGjx/fva8CAADErfP+PSO9gd8zAgBA/Onx3zMCAADQHYgRAABgFDECAACMIkYAAIBRxAgAADCKGAEAAEYRIwAAwChiBAAAGEWMAAAAo4gRAABgFDECAACMIkYAAIBRxAgAADCKGAEAAEYRIwAAwChiBAAAGEWMAAAAo4gRAABgFDECAACMIkYAAIBRxAgAADCKGAEAAEYRIwAAwChiBAAAGEWMAAAAo4gRAABgFDECAACMIkYAAIBRxAgAADCKGAEAAEYRIwAAwChbMeLz+XTNNdcoNTVVaWlpKigo0KFDhzo9pry8XA6HI2pLSUk5r6EBAEDisBUjr776qoqLi/X666+rsrJSZ86c0de+9jW1tLR0epzL5VJjY2NkO378+HkNDQAAEkdfO4u3b98e9XF5ebnS0tK0d+9effnLX+7wOIfDoYyMjNgmBAAACe28nhkJBAKSpIEDB3a67tSpUxo2bJg8Ho/mzJmjAwcOdLo+FAopGAxGbQAAIDHFHCPhcFh33XWXrr32Wo0fP77DdaNHj9a6deu0efNmPffccwqHw8rJydE777zT4TE+n09utzuyeTyeWMcEAABfcA7LsqxYDrzjjjv08ssva9euXRoyZEiXjztz5ozGjh2refPmadmyZe2uCYVCCoVCkY+DwaA8Ho8CgYBcLlcs43Zo+JKtbfYdK8vv1nMAAHAhCgaDcrvd53z/tvXMyFmLFi3Sli1btHPnTlshIkn9+vXTlClTdPjw4Q7XOJ1OOZ3OWEYDAABxxta3aSzL0qJFi7Rx40a98sorGjFihO0Ttra2qr6+XpmZmbaPBQAAicfWnZHi4mKtX79emzdvVmpqqvx+vyTJ7Xarf//+kqTCwkJlZWXJ5/NJkh599FHNmDFDI0eO1EcffaQnnnhCx48f14IFC7r5pQAAgHhkK0ZWrVolSbrhhhui9j/77LP63ve+J0lqaGhQnz6f33D58MMPtXDhQvn9fg0YMEDZ2dnavXu3xo0bd36TAwCAhBDzA6y9qasPwMSCB1gBAOgZXX3/5m/TAAAAo4gRAABgFDECAACMIkYAAIBRxAgAADCKGAEAAEYRIwAAwChiBAAAGEWMAAAAo4gRAABgFDECAACMIkYAAIBRxAgAADCKGAEAAEYRIwAAwChiBAAAGEWMAAAAo4gRAABgFDECAACMIkYAAIBRxAgAADCKGAEAAEYRIwAAwChiBAAAGEWMAAAAo4gRAABgFDECAACMIkYAAIBRxAgAADCKGAEAAEYRIwAAwChbMeLz+XTNNdcoNTVVaWlpKigo0KFDh8553IYNGzRmzBilpKRowoQJ2rZtW8wDAwCAxGIrRl599VUVFxfr9ddfV2Vlpc6cOaOvfe1ramlp6fCY3bt3a968eZo/f75qa2tVUFCggoIC7d+//7yHBwAA8c9hWZYV68Hvv/++0tLS9Oqrr+rLX/5yu2vmzp2rlpYWbdmyJbJvxowZmjx5slavXt2l8wSDQbndbgUCAblcrljHbdfwJVvb7DtWlt+t5wAA4ELU1ffv83pmJBAISJIGDhzY4Zqamhrl5uZG7cvLy1NNTU2Hx4RCIQWDwagNAAAkpphjJBwO66677tK1116r8ePHd7jO7/crPT09al96err8fn+Hx/h8Prnd7sjm8XhiHRMAAHzBxRwjxcXF2r9/vyoqKrpzHklSaWmpAoFAZDtx4kS3nwMAAHwx9I3loEWLFmnLli3auXOnhgwZ0unajIwMNTU1Re1rampSRkZGh8c4nU45nc5YRgMAAHHG1p0Ry7K0aNEibdy4Ua+88opGjBhxzmO8Xq+qqqqi9lVWVsrr9dqbFAAAJCRbd0aKi4u1fv16bd68WampqZHnPtxut/r37y9JKiwsVFZWlnw+nyRp8eLFmjlzppYvX678/HxVVFRoz549WrNmTTe/FAAAEI9s3RlZtWqVAoGAbrjhBmVmZka2F154IbKmoaFBjY2NkY9zcnK0fv16rVmzRpMmTdKLL76oTZs2dfrQKwAAuHDYujPSlV9JUl1d3WbfzTffrJtvvtnOqQAAwAWCv00DAACMIkYAAIBRxAgAADCKGAEAAEYRIwAAwChiBAAAGEWMAAAAo4gRAABgFDECAACMIkYAAIBRxAgAADCKGAEAAEYRIwAAwChiBAAAGEWMAAAAo4gRAABgFDECAACMIkYAAIBRxAgAADCKGAEAAEYRIwAAwChiBAAAGEWMAAAAo4gRAABgFDECAACMIkYAAIBRxAgAADCKGAEAAEYRIwAAwChiBAAAGEWMAAAAo2zHyM6dOzV79mwNHjxYDodDmzZt6nR9dXW1HA5Hm83v98c6MwAASCC2Y6SlpUWTJk3SypUrbR136NAhNTY2Rra0tDS7pwYAAAmor90DZs2apVmzZtk+UVpami699FLbxwEAgMTWa8+MTJ48WZmZmbrpppv02muvdbo2FAopGAxGbQAAIDH1eIxkZmZq9erV+tOf/qQ//elP8ng8uuGGG7Rv374Oj/H5fHK73ZHN4/H09JgAAMAQh2VZVswHOxzauHGjCgoKbB03c+ZMDR06VL///e/b/XwoFFIoFIp8HAwG5fF4FAgE5HK5Yh23XcOXbG2z71hZfreeAwCAC1EwGJTb7T7n+7ftZ0a6w7Rp07Rr164OP+90OuV0OntxIgAAYIqR3zNSV1enzMxME6cGAABfMLbvjJw6dUqHDx+OfHz06FHV1dVp4MCBGjp0qEpLS/Xuu+/qd7/7nSRpxYoVGjFihK666ip9+umnWrt2rV555RX95S9/6b5XAQAA4pbtGNmzZ49uvPHGyMclJSWSpKKiIpWXl6uxsVENDQ2Rz58+fVo/+tGP9O677+qiiy7SxIkT9de//jXqawAAgAvXeT3A2lu6+gBMLHiAFQCAntHV92/+Ng0AADCKGAEAAEYRIwAAwChiBAAAGEWMAAAAo4gRAABgFDECAACMIkYAAIBRxAgAADCKGAEAAEYRIwAAwChiBAAAGEWMAAAAo4gRAABgFDECAACMIkYAAIBRxAgAADCKGAEAAEYRIwAAwChiBAAAGEWMAAAAo4gRAABgFDECAACMIkYAAIBRxAgAADCKGAEAAEYRIwAAwChiBAAAGEWMAAAAo4gRAABgFDECAACMsh0jO3fu1OzZszV48GA5HA5t2rTpnMdUV1fr6quvltPp1MiRI1VeXh7DqAAAIBHZjpGWlhZNmjRJK1eu7NL6o0ePKj8/XzfeeKPq6up01113acGCBdqxY4ftYQEAQOLpa/eAWbNmadasWV1ev3r1ao0YMULLly+XJI0dO1a7du3SL3/5S+Xl5dk9PQAASDA9/sxITU2NcnNzo/bl5eWppqamw2NCoZCCwWDUBgAAElOPx4jf71d6enrUvvT0dAWDQX3yySftHuPz+eR2uyObx+Pp6TEBAIAhX8ifpiktLVUgEIhsJ06cMD0SAADoIbafGbErIyNDTU1NUfuamprkcrnUv3//do9xOp1yOp09PRoAAPgC6PE7I16vV1VVVVH7Kisr5fV6e/rUAAAgDtiOkVOnTqmurk51dXWS/vuju3V1dWpoaJD032+xFBYWRtbffvvtOnLkiO677z69/fbb+vWvf60//vGPuvvuu7vnFQAAgLhmO0b27NmjKVOmaMqUKZKkkpISTZkyRQ899JAkqbGxMRImkjRixAht3bpVlZWVmjRpkpYvX661a9fyY70AAECS5LAsyzI9xLkEg0G53W4FAgG5XK5u/drDl2xts+9YWX63ngMAgAtRV9+/v5A/TQMAAC4cxAgAADCKGAEAAEYRIwAAwChiBAAAGEWMAAAAo4gRAABgFDECAACMIkYAAIBRxAgAADCKGAEAAEYRIwAAwChiBAAAGEWMAAAAo4gRAABgFDECAACMIkYAAIBRxAgAADCKGAEAAEYRIwAAwChiBAAAGEWMAAAAo4gRAABgFDECAACMIkYAAIBRxAgAADCKGAEAAEYRIwAAwChiBAAAGEWMAAAAo4gRAABgVEwxsnLlSg0fPlwpKSmaPn263njjjQ7XlpeXy+FwRG0pKSkxDwwAABKL7Rh54YUXVFJSoqVLl2rfvn2aNGmS8vLydPLkyQ6PcblcamxsjGzHjx8/r6EBAEDisB0jv/jFL7Rw4ULddtttGjdunFavXq2LLrpI69at6/AYh8OhjIyMyJaenn5eQwMAgMRhK0ZOnz6tvXv3Kjc39/Mv0KePcnNzVVNT0+Fxp06d0rBhw+TxeDRnzhwdOHCg0/OEQiEFg8GoDQAAJCZbMfLBBx+otbW1zZ2N9PR0+f3+do8ZPXq01q1bp82bN+u5555TOBxWTk6O3nnnnQ7P4/P55Ha7I5vH47EzJgAAiCM9/tM0Xq9XhYWFmjx5smbOnKmXXnpJX/rSl/T00093eExpaakCgUBkO3HiRE+PCQAADOlrZ/GgQYOUlJSkpqamqP1NTU3KyMjo0tfo16+fpkyZosOHD3e4xul0yul02hkNAADEKVt3RpKTk5Wdna2qqqrIvnA4rKqqKnm93i59jdbWVtXX1yszM9PepAAAICHZujMiSSUlJSoqKtLUqVM1bdo0rVixQi0tLbrtttskSYWFhcrKypLP55MkPfroo5oxY4ZGjhypjz76SE888YSOHz+uBQsWdO8rAQAAccl2jMydO1fvv/++HnroIfn9fk2ePFnbt2+PPNTa0NCgPn0+v+Hy4YcfauHChfL7/RowYICys7O1e/dujRs3rvteBQAAiFsOy7Is00OcSzAYlNvtViAQkMvl6tavPXzJ1jb7jpXld+s5AAC4EHX1/Zu/TQMAAIwiRgAAgFHECAAAMIoYAQAARhEjAADAKGIEAAAYRYwAAACjiBEAAGAUMQIAAIwiRgAAgFHECAAAMIoYAQAARhEjAADAKGIEAAAYRYwAAACjiBEAAGAUMQIAAIwiRgAAgFHECAAAMIoYAQAARhEjAADAKGIEAAAYRYwAAACjiBEAAGAUMQIAAIwiRgAAgFHECAAAMIoYAQAARhEjAADAKGIEAAAYRYwAAACjYoqRlStXavjw4UpJSdH06dP1xhtvdLp+w4YNGjNmjFJSUjRhwgRt27YtpmEBAEDisR0jL7zwgkpKSrR06VLt27dPkyZNUl5enk6ePNnu+t27d2vevHmaP3++amtrVVBQoIKCAu3fv/+8hwcAAPHPYVmWZeeA6dOn65prrtFTTz0lSQqHw/J4PPrhD3+oJUuWtFk/d+5ctbS0aMuWLZF9M2bM0OTJk7V69eounTMYDMrtdisQCMjlctkZ95yGL9l6zjXHyvK79ZwAAFwIuvr+3dfOFz19+rT27t2r0tLSyL4+ffooNzdXNTU17R5TU1OjkpKSqH15eXnatGlTh+cJhUIKhUKRjwOBgKT/vqjuFg59fM41Q+/e0Gbf/kfyun0WAAASydn37XPd97AVIx988IFaW1uVnp4etT89PV1vv/12u8f4/f521/v9/g7P4/P59Mgjj7TZ7/F47Izbo9wrTE8AAEB8aG5ultvt7vDztmKkt5SWlkbdTQmHw/rPf/6jyy67TA6Ho9vOEwwG5fF4dOLEiW7/9g8+x3XuPVzr3sF17h1c597Rk9fZsiw1Nzdr8ODBna6zFSODBg1SUlKSmpqaovY3NTUpIyOj3WMyMjJsrZckp9Mpp9MZte/SSy+1M6otLpeLf+i9gOvce7jWvYPr3Du4zr2jp65zZ3dEzrL10zTJycnKzs5WVVVVZF84HFZVVZW8Xm+7x3i93qj1klRZWdnhegAAcGGx/W2akpISFRUVaerUqZo2bZpWrFihlpYW3XbbbZKkwsJCZWVlyefzSZIWL16smTNnavny5crPz1dFRYX27NmjNWvWdO8rAQAAccl2jMydO1fvv/++HnroIfn9fk2ePFnbt2+PPKTa0NCgPn0+v+GSk5Oj9evX68EHH9T999+vK6+8Ups2bdL48eO771XEyOl0aunSpW2+JYTuxXXuPVzr3sF17h1c597xRbjOtn/PCAAAQHfib9MAAACjiBEAAGAUMQIAAIwiRgAAgFEJHyMrV67U8OHDlZKSounTp+uNN97odP2GDRs0ZswYpaSkaMKECdq2bVsvTRrf7FznZ555Rtdff70GDBigAQMGKDc395z/XfA5u/+mz6qoqJDD4VBBQUHPDpgg7F7njz76SMXFxcrMzJTT6dSoUaP4/0cX2L3OK1as0OjRo9W/f395PB7dfffd+vTTT3tp2vi0c+dOzZ49W4MHD5bD4ej0b8OdVV1drauvvlpOp1MjR45UeXl5zw5pJbCKigorOTnZWrdunXXgwAFr4cKF1qWXXmo1NTW1u/61116zkpKSrJ/97GfWwYMHrQcffNDq16+fVV9f38uTxxe71/m73/2utXLlSqu2ttZ66623rO9973uW2+223nnnnV6ePP7YvdZnHT161MrKyrKuv/56a86cOb0zbByze51DoZA1depU6xvf+Ia1a9cu6+jRo1Z1dbVVV1fXy5PHF7vX+fnnn7ecTqf1/PPPW0ePHrV27NhhZWZmWnfffXcvTx5ftm3bZj3wwAPWSy+9ZEmyNm7c2On6I0eOWBdddJFVUlJiHTx40PrVr35lJSUlWdu3b++xGRM6RqZNm2YVFxdHPm5tbbUGDx5s+Xy+dtd/5zvfsfLz86P2TZ8+3frBD37Qo3PGO7vX+X999tlnVmpqqvXb3/62p0ZMGLFc688++8zKycmx1q5daxUVFREjXWD3Oq9atcq6/PLLrdOnT/fWiAnB7nUuLi62vvKVr0TtKykpsa699toenTORdCVG7rvvPuuqq66K2jd37lwrLy+vx+ZK2G/TnD59Wnv37lVubm5kX58+fZSbm6uampp2j6mpqYlaL0l5eXkdrkds1/l/ffzxxzpz5owGDhzYU2MmhFiv9aOPPqq0tDTNnz+/N8aMe7Fc5z//+c/yer0qLi5Wenq6xo8fr5/+9KdqbW3trbHjTizXOScnR3v37o18K+fIkSPatm2bvvGNb/TKzBcKE++FX8i/2tsdPvjgA7W2tkZ+M+xZ6enpevvtt9s9xu/3t7ve7/f32JzxLpbr/L9+/OMfa/DgwW3+8SNaLNd6165d+s1vfqO6urpemDAxxHKdjxw5oldeeUW33HKLtm3bpsOHD+vOO+/UmTNntHTp0t4YO+7Ecp2/+93v6oMPPtB1110ny7L02Wef6fbbb9f999/fGyNfMDp6LwwGg/rkk0/Uv3//bj9nwt4ZQXwoKytTRUWFNm7cqJSUFNPjJJTm5mbdeuuteuaZZzRo0CDT4yS0cDistLQ0rVmzRtnZ2Zo7d64eeOABrV692vRoCaW6ulo//elP9etf/1r79u3TSy+9pK1bt2rZsmWmR8N5Stg7I4MGDVJSUpKampqi9jc1NSkjI6PdYzIyMmytR2zX+ayf//znKisr01//+ldNnDixJ8dMCHav9b///W8dO3ZMs2fPjuwLh8OSpL59++rQoUO64oorenboOBTLv+nMzEz169dPSUlJkX1jx46V3+/X6dOnlZyc3KMzx6NYrvNPfvIT3XrrrVqwYIEkacKECWppadH3v/99PfDAA1F/Fw2x6+i90OVy9chdESmB74wkJycrOztbVVVVkX3hcFhVVVXyer3tHuP1eqPWS1JlZWWH6xHbdZakn/3sZ1q2bJm2b9+uqVOn9saocc/utR4zZozq6+tVV1cX2b71rW/pxhtvVF1dnTweT2+OHzdi+Td97bXX6vDhw5HYk6R//vOfyszMJEQ6EMt1/vjjj9sEx9kAtPgza93GyHthjz0a+wVQUVFhOZ1Oq7y83Dp48KD1/e9/37r00kstv99vWZZl3XrrrdaSJUsi61977TWrb9++1s9//nPrrbfespYuXcqP9naB3etcVlZmJScnWy+++KLV2NgY2Zqbm029hLhh91r/L36apmvsXueGhgYrNTXVWrRokXXo0CFry5YtVlpamvXYY4+Zeglxwe51Xrp0qZWammr94Q9/sI4cOWL95S9/sa644grrO9/5jqmXEBeam5ut2tpaq7a21pJk/eIXv7Bqa2ut48ePW5ZlWUuWLLFuvfXWyPqzP9p77733Wm+99Za1cuVKfrT3fP3qV7+yhg4daiUnJ1vTpk2zXn/99cjnZs6caRUVFUWt/+Mf/2iNGjXKSk5Otq666ipr69atvTxxfLJznYcNG2ZJarMtXbq09wePQ3b/Tf9/xEjX2b3Ou3fvtqZPn245nU7r8ssvtx5//HHrs88+6+Wp44+d63zmzBnr4Ycftq644gorJSXF8ng81p133ml9+OGHvT94HPnb3/7W7v9zz17boqIia+bMmW2OmTx5spWcnGxdfvnl1rPPPtujMzosi3tbAADAnIR9ZgQAAMQHYgQAABhFjAAAAKOIEQAAYBQxAgAAjCJGAACAUcQIAAAwihgBAABGESMAAMAoYgQAABhFjAAAAKOIEQAAYNT/AdSw95C5gfz9AAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "print(normalised_imag.min())\n",
    "print(normalised_imag.max())\n",
    "def plotHistNorm(imag):\n",
    "    plt.hist(imag.ravel(), bins=100)\n",
    "    plt.show()\n",
    "plotHistNorm(normalised_imag)\n",
    "# auto : "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "bc729a00",
   "metadata": {},
   "outputs": [],
   "source": [
    "def saveImgToCsv(imag, filename, lon_min=-180, lon_max=180, lat_min=-90, lat_max=+90):\n",
    "    # Generate lat/lon values for each point\n",
    "    latitudes = np.linspace(lat_min, lat_max, imag.shape[0])\n",
    "    longitudes = np.linspace(lon_min, lon_max, imag.shape[1])\n",
    "\n",
    "    # Create a grid of coordinates\n",
    "    lon_grid, lat_grid = np.meshgrid(longitudes, latitudes)\n",
    "\n",
    "    # Flatten the arrays to create columns\n",
    "    lat_flat = lat_grid.flatten()\n",
    "    lon_flat = lon_grid.flatten()\n",
    "    value_flat = np.ma.filled(imag.flatten(), np.nan)  # Convert masked values to NaN\n",
    "\n",
    "    # Create a DataFrame from the flattened arrays\n",
    "    data = {\n",
    "        'Latitude': lat_flat,\n",
    "        'Longitude': lon_flat,\n",
    "        'Value': value_flat\n",
    "    }\n",
    "\n",
    "    df = pd.DataFrame(data)\n",
    "\n",
    "    df_clean = df.dropna()\n",
    "    df_clean.to_csv(filename, index=False)\n",
    "    # np.savetxt(filename, imag, delimiter=\",\")\n",
    "\n",
    "\n",
    "# saveImgToCsv(imag, 'World.csv')\n",
    "# 12 min 1.2GB\n",
    "# 15 min 1.44GB\n",
    "# 18 min 1.76GB\n",
    "# 21 min 2.12GB\n",
    "# 28 min 2.88GB\n",
    "# \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "29a1b258",
   "metadata": {},
   "outputs": [],
   "source": [
    "# saveImgToCsv(imag, 'World.csv')\n",
    "# 12 min 1.2GB\n",
    "# "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "2c983846",
   "metadata": {},
   "outputs": [],
   "source": [
    "def toCsvToImag(filename, lon_min=-180, lon_max=180, lat_min=-90, lat_max=+90):\n",
    "    df = pd.read_csv(filename)\n",
    "    lon_step = 360 / 43200\n",
    "    lat_step = 180 / 17406\n",
    "    latitudes = np.arange(lat_min, lat_max, lat_step)\n",
    "    longitudes = np.arange(lon_min, lon_max, lon_step)\n",
    "    imag = np.zeros((len(latitudes), len(longitudes)))\n",
    "    for i in range(len(df)):\n",
    "        lat = df['Latitude'][i]\n",
    "        lon = df['Longitude'][i]\n",
    "        value = df['Value'][i]\n",
    "        lat_index = np.where(latitudes == lat)\n",
    "        lon_index = np.where(longitudes == lon)\n",
    "        imag[lat_index, lon_index] = value\n",
    "    return imag\n",
    "# World = toCsvToImag('World.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "d04dbada-ff06-4ac4-be96-8eaa3e4397c7",
   "metadata": {},
   "outputs": [],
   "source": [
    "#plt.imshow(np.log(imag[:17406/2, :17406/2]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "79a790d1-c025-4bb2-a9aa-aa273e9760db",
   "metadata": {},
   "outputs": [],
   "source": [
    "# plt.imshow(np.log(imag))\n",
    "# plt.colorbar()\n",
    "# 1 min"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "d687de69",
   "metadata": {},
   "outputs": [],
   "source": [
    "def getTransform(file_name):\n",
    "    with rasterio.open(file_name) as image:\n",
    "        # Extract the transformation and use it to compute the geographical coordinates\n",
    "        transform = image.transform\n",
    "    return transform\n",
    "transform = getTransform(file_name)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "3dfea615-ef2b-459e-b6cf-204397034443",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/var/folders/5g/9cr0t2sj5mj8w4wcjj4by3300000gq/T/ipykernel_60610/496434615.py:15: RuntimeWarning: divide by zero encountered in log\n",
      "  imag_log = np.log(imag[y_pixels_min:y_pixels_max, x_pixels_min:x_pixels_max])\n"
     ]
    }
   ],
   "source": [
    "def get_log_imag(imag, transform, lon_min = -180, lat_min = -90, lon_max = 180, lat_max = 90):\n",
    "    inv_transform = ~transform\n",
    "    x_pixels_min, y_pixels_min = inv_transform * (lon_min, lat_max)\n",
    "    x_pixels_max, y_pixels_max = inv_transform * (lon_max, lat_min)\n",
    "\n",
    "    # Round and convert to integer\n",
    "    x_pixels_min, y_pixels_min = int(np.floor(x_pixels_min)), int(np.floor(y_pixels_min))\n",
    "    x_pixels_max, y_pixels_max = int(np.ceil(x_pixels_max)), int(np.ceil(y_pixels_max))\n",
    "\n",
    "    # Ensure the slicing indices are within the image bounds\n",
    "    x_pixels_min, y_pixels_min = max(x_pixels_min, 0), max(y_pixels_min, 0)\n",
    "    x_pixels_max, y_pixels_max = min(x_pixels_max, imag.shape[1]), min(y_pixels_max, imag.shape[0])\n",
    "\n",
    "    # Log transformation\n",
    "    imag_log = np.log(imag[y_pixels_min:y_pixels_max, x_pixels_min:x_pixels_max])\n",
    "    imag_log = np.ma.masked_invalid(imag_log)\n",
    "    boundary = np.array([lon_min, lon_max, lat_min, lat_max])\n",
    "    return imag_log, boundary\n",
    "\n",
    "World, bWorld = get_log_imag(imag, transform)\n",
    "Europe, bEurope = get_log_imag(imag, transform, -25, 35, 45, 70)\n",
    "Africa, bAfrica = get_log_imag(imag, transform, -25, -35, 45, 35)\n",
    "South_America, bSouth_America = get_log_imag(imag, transform, -90, -60, -30, 15)\n",
    "North_America, bNorth_America = get_log_imag(imag, transform, -180, 15, -30, 90)\n",
    "Australia, bAustralia = get_log_imag(imag, transform, 110, -45, 160, -10)\n",
    "Asia, bAsia = get_log_imag(imag, transform, 40, 0, 180, 70)\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "478590b8",
   "metadata": {},
   "outputs": [
    {
     "ename": "KeyboardInterrupt",
     "evalue": "",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mKeyboardInterrupt\u001b[0m                         Traceback (most recent call last)",
      "Cell \u001b[0;32mIn[17], line 22\u001b[0m\n\u001b[1;32m     20\u001b[0m plotlogimag(Europe, bEurope, \u001b[38;5;124m'\u001b[39m\u001b[38;5;124mEurope\u001b[39m\u001b[38;5;124m'\u001b[39m)\n\u001b[1;32m     21\u001b[0m plotlogimag(Africa, bAfrica, \u001b[38;5;124m'\u001b[39m\u001b[38;5;124mAfrica\u001b[39m\u001b[38;5;124m'\u001b[39m)\n\u001b[0;32m---> 22\u001b[0m \u001b[43mplotlogimag\u001b[49m\u001b[43m(\u001b[49m\u001b[43mSouth_America\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mbSouth_America\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;124;43m'\u001b[39;49m\u001b[38;5;124;43mSouth America\u001b[39;49m\u001b[38;5;124;43m'\u001b[39;49m\u001b[43m)\u001b[49m\n\u001b[1;32m     23\u001b[0m plotlogimag(North_America, bNorth_America, \u001b[38;5;124m'\u001b[39m\u001b[38;5;124mNorth America\u001b[39m\u001b[38;5;124m'\u001b[39m)\n\u001b[1;32m     24\u001b[0m plotlogimag(Australia, bAustralia, \u001b[38;5;124m'\u001b[39m\u001b[38;5;124mAustralia\u001b[39m\u001b[38;5;124m'\u001b[39m)\n",
      "Cell \u001b[0;32mIn[17], line 18\u001b[0m, in \u001b[0;36mplotlogimag\u001b[0;34m(imag_log, boundary, title)\u001b[0m\n\u001b[1;32m     16\u001b[0m ax\u001b[38;5;241m.\u001b[39mgrid()\n\u001b[1;32m     17\u001b[0m fig\u001b[38;5;241m.\u001b[39mcolorbar(cax, ax\u001b[38;5;241m=\u001b[39max, label\u001b[38;5;241m=\u001b[39m\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mLog Brightness\u001b[39m\u001b[38;5;124m'\u001b[39m)\n\u001b[0;32m---> 18\u001b[0m \u001b[43mfig\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43msavefig\u001b[49m\u001b[43m(\u001b[49m\u001b[38;5;124;43m'\u001b[39;49m\u001b[38;5;124;43mWorldPop2015/\u001b[39;49m\u001b[38;5;124;43m'\u001b[39;49m\u001b[38;5;241;43m+\u001b[39;49m\u001b[43mtitle\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mdpi\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43m \u001b[49m\u001b[38;5;241;43m200\u001b[39;49m\u001b[43m)\u001b[49m\n",
      "File \u001b[0;32m/usr/local/lib/python3.9/site-packages/matplotlib/figure.py:3390\u001b[0m, in \u001b[0;36mFigure.savefig\u001b[0;34m(self, fname, transparent, **kwargs)\u001b[0m\n\u001b[1;32m   3388\u001b[0m     \u001b[38;5;28;01mfor\u001b[39;00m ax \u001b[38;5;129;01min\u001b[39;00m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39maxes:\n\u001b[1;32m   3389\u001b[0m         _recursively_make_axes_transparent(stack, ax)\n\u001b[0;32m-> 3390\u001b[0m \u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mcanvas\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mprint_figure\u001b[49m\u001b[43m(\u001b[49m\u001b[43mfname\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[43mkwargs\u001b[49m\u001b[43m)\u001b[49m\n",
      "File \u001b[0;32m/usr/local/lib/python3.9/site-packages/matplotlib/backend_bases.py:2193\u001b[0m, in \u001b[0;36mFigureCanvasBase.print_figure\u001b[0;34m(self, filename, dpi, facecolor, edgecolor, orientation, format, bbox_inches, pad_inches, bbox_extra_artists, backend, **kwargs)\u001b[0m\n\u001b[1;32m   2189\u001b[0m \u001b[38;5;28;01mtry\u001b[39;00m:\n\u001b[1;32m   2190\u001b[0m     \u001b[38;5;66;03m# _get_renderer may change the figure dpi (as vector formats\u001b[39;00m\n\u001b[1;32m   2191\u001b[0m     \u001b[38;5;66;03m# force the figure dpi to 72), so we need to set it again here.\u001b[39;00m\n\u001b[1;32m   2192\u001b[0m     \u001b[38;5;28;01mwith\u001b[39;00m cbook\u001b[38;5;241m.\u001b[39m_setattr_cm(\u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mfigure, dpi\u001b[38;5;241m=\u001b[39mdpi):\n\u001b[0;32m-> 2193\u001b[0m         result \u001b[38;5;241m=\u001b[39m \u001b[43mprint_method\u001b[49m\u001b[43m(\u001b[49m\n\u001b[1;32m   2194\u001b[0m \u001b[43m            \u001b[49m\u001b[43mfilename\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m   2195\u001b[0m \u001b[43m            \u001b[49m\u001b[43mfacecolor\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mfacecolor\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m   2196\u001b[0m \u001b[43m            \u001b[49m\u001b[43medgecolor\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43medgecolor\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m   2197\u001b[0m \u001b[43m            \u001b[49m\u001b[43morientation\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43morientation\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m   2198\u001b[0m \u001b[43m            \u001b[49m\u001b[43mbbox_inches_restore\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43m_bbox_inches_restore\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m   2199\u001b[0m \u001b[43m            \u001b[49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[43mkwargs\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m   2200\u001b[0m \u001b[38;5;28;01mfinally\u001b[39;00m:\n\u001b[1;32m   2201\u001b[0m     \u001b[38;5;28;01mif\u001b[39;00m bbox_inches \u001b[38;5;129;01mand\u001b[39;00m restore_bbox:\n",
      "File \u001b[0;32m/usr/local/lib/python3.9/site-packages/matplotlib/backend_bases.py:2043\u001b[0m, in \u001b[0;36mFigureCanvasBase._switch_canvas_and_return_print_method.<locals>.<lambda>\u001b[0;34m(*args, **kwargs)\u001b[0m\n\u001b[1;32m   2039\u001b[0m     optional_kws \u001b[38;5;241m=\u001b[39m {  \u001b[38;5;66;03m# Passed by print_figure for other renderers.\u001b[39;00m\n\u001b[1;32m   2040\u001b[0m         \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mdpi\u001b[39m\u001b[38;5;124m\"\u001b[39m, \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mfacecolor\u001b[39m\u001b[38;5;124m\"\u001b[39m, \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124medgecolor\u001b[39m\u001b[38;5;124m\"\u001b[39m, \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124morientation\u001b[39m\u001b[38;5;124m\"\u001b[39m,\n\u001b[1;32m   2041\u001b[0m         \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mbbox_inches_restore\u001b[39m\u001b[38;5;124m\"\u001b[39m}\n\u001b[1;32m   2042\u001b[0m     skip \u001b[38;5;241m=\u001b[39m optional_kws \u001b[38;5;241m-\u001b[39m {\u001b[38;5;241m*\u001b[39minspect\u001b[38;5;241m.\u001b[39msignature(meth)\u001b[38;5;241m.\u001b[39mparameters}\n\u001b[0;32m-> 2043\u001b[0m     print_method \u001b[38;5;241m=\u001b[39m functools\u001b[38;5;241m.\u001b[39mwraps(meth)(\u001b[38;5;28;01mlambda\u001b[39;00m \u001b[38;5;241m*\u001b[39margs, \u001b[38;5;241m*\u001b[39m\u001b[38;5;241m*\u001b[39mkwargs: \u001b[43mmeth\u001b[49m\u001b[43m(\u001b[49m\n\u001b[1;32m   2044\u001b[0m \u001b[43m        \u001b[49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[43margs\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[43m{\u001b[49m\u001b[43mk\u001b[49m\u001b[43m:\u001b[49m\u001b[43m \u001b[49m\u001b[43mv\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;28;43;01mfor\u001b[39;49;00m\u001b[43m \u001b[49m\u001b[43mk\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mv\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;129;43;01min\u001b[39;49;00m\u001b[43m \u001b[49m\u001b[43mkwargs\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mitems\u001b[49m\u001b[43m(\u001b[49m\u001b[43m)\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;28;43;01mif\u001b[39;49;00m\u001b[43m \u001b[49m\u001b[43mk\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;129;43;01mnot\u001b[39;49;00m\u001b[43m \u001b[49m\u001b[38;5;129;43;01min\u001b[39;49;00m\u001b[43m \u001b[49m\u001b[43mskip\u001b[49m\u001b[43m}\u001b[49m\u001b[43m)\u001b[49m)\n\u001b[1;32m   2045\u001b[0m \u001b[38;5;28;01melse\u001b[39;00m:  \u001b[38;5;66;03m# Let third-parties do as they see fit.\u001b[39;00m\n\u001b[1;32m   2046\u001b[0m     print_method \u001b[38;5;241m=\u001b[39m meth\n",
      "File \u001b[0;32m/usr/local/lib/python3.9/site-packages/matplotlib/backends/backend_agg.py:497\u001b[0m, in \u001b[0;36mFigureCanvasAgg.print_png\u001b[0;34m(self, filename_or_obj, metadata, pil_kwargs)\u001b[0m\n\u001b[1;32m    450\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21mprint_png\u001b[39m(\u001b[38;5;28mself\u001b[39m, filename_or_obj, \u001b[38;5;241m*\u001b[39m, metadata\u001b[38;5;241m=\u001b[39m\u001b[38;5;28;01mNone\u001b[39;00m, pil_kwargs\u001b[38;5;241m=\u001b[39m\u001b[38;5;28;01mNone\u001b[39;00m):\n\u001b[1;32m    451\u001b[0m \u001b[38;5;250m    \u001b[39m\u001b[38;5;124;03m\"\"\"\u001b[39;00m\n\u001b[1;32m    452\u001b[0m \u001b[38;5;124;03m    Write the figure to a PNG file.\u001b[39;00m\n\u001b[1;32m    453\u001b[0m \n\u001b[0;32m   (...)\u001b[0m\n\u001b[1;32m    495\u001b[0m \u001b[38;5;124;03m        *metadata*, including the default 'Software' key.\u001b[39;00m\n\u001b[1;32m    496\u001b[0m \u001b[38;5;124;03m    \"\"\"\u001b[39;00m\n\u001b[0;32m--> 497\u001b[0m     \u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43m_print_pil\u001b[49m\u001b[43m(\u001b[49m\u001b[43mfilename_or_obj\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;124;43m\"\u001b[39;49m\u001b[38;5;124;43mpng\u001b[39;49m\u001b[38;5;124;43m\"\u001b[39;49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mpil_kwargs\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mmetadata\u001b[49m\u001b[43m)\u001b[49m\n",
      "File \u001b[0;32m/usr/local/lib/python3.9/site-packages/matplotlib/backends/backend_agg.py:445\u001b[0m, in \u001b[0;36mFigureCanvasAgg._print_pil\u001b[0;34m(self, filename_or_obj, fmt, pil_kwargs, metadata)\u001b[0m\n\u001b[1;32m    440\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21m_print_pil\u001b[39m(\u001b[38;5;28mself\u001b[39m, filename_or_obj, fmt, pil_kwargs, metadata\u001b[38;5;241m=\u001b[39m\u001b[38;5;28;01mNone\u001b[39;00m):\n\u001b[1;32m    441\u001b[0m \u001b[38;5;250m    \u001b[39m\u001b[38;5;124;03m\"\"\"\u001b[39;00m\n\u001b[1;32m    442\u001b[0m \u001b[38;5;124;03m    Draw the canvas, then save it using `.image.imsave` (to which\u001b[39;00m\n\u001b[1;32m    443\u001b[0m \u001b[38;5;124;03m    *pil_kwargs* and *metadata* are forwarded).\u001b[39;00m\n\u001b[1;32m    444\u001b[0m \u001b[38;5;124;03m    \"\"\"\u001b[39;00m\n\u001b[0;32m--> 445\u001b[0m     \u001b[43mFigureCanvasAgg\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mdraw\u001b[49m\u001b[43m(\u001b[49m\u001b[38;5;28;43mself\u001b[39;49m\u001b[43m)\u001b[49m\n\u001b[1;32m    446\u001b[0m     mpl\u001b[38;5;241m.\u001b[39mimage\u001b[38;5;241m.\u001b[39mimsave(\n\u001b[1;32m    447\u001b[0m         filename_or_obj, \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mbuffer_rgba(), \u001b[38;5;28mformat\u001b[39m\u001b[38;5;241m=\u001b[39mfmt, origin\u001b[38;5;241m=\u001b[39m\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mupper\u001b[39m\u001b[38;5;124m\"\u001b[39m,\n\u001b[1;32m    448\u001b[0m         dpi\u001b[38;5;241m=\u001b[39m\u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mfigure\u001b[38;5;241m.\u001b[39mdpi, metadata\u001b[38;5;241m=\u001b[39mmetadata, pil_kwargs\u001b[38;5;241m=\u001b[39mpil_kwargs)\n",
      "File \u001b[0;32m/usr/local/lib/python3.9/site-packages/matplotlib/backends/backend_agg.py:388\u001b[0m, in \u001b[0;36mFigureCanvasAgg.draw\u001b[0;34m(self)\u001b[0m\n\u001b[1;32m    385\u001b[0m \u001b[38;5;66;03m# Acquire a lock on the shared font cache.\u001b[39;00m\n\u001b[1;32m    386\u001b[0m \u001b[38;5;28;01mwith\u001b[39;00m (\u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mtoolbar\u001b[38;5;241m.\u001b[39m_wait_cursor_for_draw_cm() \u001b[38;5;28;01mif\u001b[39;00m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mtoolbar\n\u001b[1;32m    387\u001b[0m       \u001b[38;5;28;01melse\u001b[39;00m nullcontext()):\n\u001b[0;32m--> 388\u001b[0m     \u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mfigure\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mdraw\u001b[49m\u001b[43m(\u001b[49m\u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mrenderer\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m    389\u001b[0m     \u001b[38;5;66;03m# A GUI class may be need to update a window using this draw, so\u001b[39;00m\n\u001b[1;32m    390\u001b[0m     \u001b[38;5;66;03m# don't forget to call the superclass.\u001b[39;00m\n\u001b[1;32m    391\u001b[0m     \u001b[38;5;28msuper\u001b[39m()\u001b[38;5;241m.\u001b[39mdraw()\n",
      "File \u001b[0;32m/usr/local/lib/python3.9/site-packages/matplotlib/artist.py:95\u001b[0m, in \u001b[0;36m_finalize_rasterization.<locals>.draw_wrapper\u001b[0;34m(artist, renderer, *args, **kwargs)\u001b[0m\n\u001b[1;32m     93\u001b[0m \u001b[38;5;129m@wraps\u001b[39m(draw)\n\u001b[1;32m     94\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21mdraw_wrapper\u001b[39m(artist, renderer, \u001b[38;5;241m*\u001b[39margs, \u001b[38;5;241m*\u001b[39m\u001b[38;5;241m*\u001b[39mkwargs):\n\u001b[0;32m---> 95\u001b[0m     result \u001b[38;5;241m=\u001b[39m \u001b[43mdraw\u001b[49m\u001b[43m(\u001b[49m\u001b[43martist\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mrenderer\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[43margs\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[43mkwargs\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m     96\u001b[0m     \u001b[38;5;28;01mif\u001b[39;00m renderer\u001b[38;5;241m.\u001b[39m_rasterizing:\n\u001b[1;32m     97\u001b[0m         renderer\u001b[38;5;241m.\u001b[39mstop_rasterizing()\n",
      "File \u001b[0;32m/usr/local/lib/python3.9/site-packages/matplotlib/artist.py:72\u001b[0m, in \u001b[0;36mallow_rasterization.<locals>.draw_wrapper\u001b[0;34m(artist, renderer)\u001b[0m\n\u001b[1;32m     69\u001b[0m     \u001b[38;5;28;01mif\u001b[39;00m artist\u001b[38;5;241m.\u001b[39mget_agg_filter() \u001b[38;5;129;01mis\u001b[39;00m \u001b[38;5;129;01mnot\u001b[39;00m \u001b[38;5;28;01mNone\u001b[39;00m:\n\u001b[1;32m     70\u001b[0m         renderer\u001b[38;5;241m.\u001b[39mstart_filter()\n\u001b[0;32m---> 72\u001b[0m     \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[43mdraw\u001b[49m\u001b[43m(\u001b[49m\u001b[43martist\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mrenderer\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m     73\u001b[0m \u001b[38;5;28;01mfinally\u001b[39;00m:\n\u001b[1;32m     74\u001b[0m     \u001b[38;5;28;01mif\u001b[39;00m artist\u001b[38;5;241m.\u001b[39mget_agg_filter() \u001b[38;5;129;01mis\u001b[39;00m \u001b[38;5;129;01mnot\u001b[39;00m \u001b[38;5;28;01mNone\u001b[39;00m:\n",
      "File \u001b[0;32m/usr/local/lib/python3.9/site-packages/matplotlib/figure.py:3154\u001b[0m, in \u001b[0;36mFigure.draw\u001b[0;34m(self, renderer)\u001b[0m\n\u001b[1;32m   3151\u001b[0m         \u001b[38;5;66;03m# ValueError can occur when resizing a window.\u001b[39;00m\n\u001b[1;32m   3153\u001b[0m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mpatch\u001b[38;5;241m.\u001b[39mdraw(renderer)\n\u001b[0;32m-> 3154\u001b[0m \u001b[43mmimage\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43m_draw_list_compositing_images\u001b[49m\u001b[43m(\u001b[49m\n\u001b[1;32m   3155\u001b[0m \u001b[43m    \u001b[49m\u001b[43mrenderer\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;28;43mself\u001b[39;49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43martists\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43msuppressComposite\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m   3157\u001b[0m \u001b[38;5;28;01mfor\u001b[39;00m sfig \u001b[38;5;129;01min\u001b[39;00m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39msubfigs:\n\u001b[1;32m   3158\u001b[0m     sfig\u001b[38;5;241m.\u001b[39mdraw(renderer)\n",
      "File \u001b[0;32m/usr/local/lib/python3.9/site-packages/matplotlib/image.py:132\u001b[0m, in \u001b[0;36m_draw_list_compositing_images\u001b[0;34m(renderer, parent, artists, suppress_composite)\u001b[0m\n\u001b[1;32m    130\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m not_composite \u001b[38;5;129;01mor\u001b[39;00m \u001b[38;5;129;01mnot\u001b[39;00m has_images:\n\u001b[1;32m    131\u001b[0m     \u001b[38;5;28;01mfor\u001b[39;00m a \u001b[38;5;129;01min\u001b[39;00m artists:\n\u001b[0;32m--> 132\u001b[0m         \u001b[43ma\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mdraw\u001b[49m\u001b[43m(\u001b[49m\u001b[43mrenderer\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m    133\u001b[0m \u001b[38;5;28;01melse\u001b[39;00m:\n\u001b[1;32m    134\u001b[0m     \u001b[38;5;66;03m# Composite any adjacent images together\u001b[39;00m\n\u001b[1;32m    135\u001b[0m     image_group \u001b[38;5;241m=\u001b[39m []\n",
      "File \u001b[0;32m/usr/local/lib/python3.9/site-packages/matplotlib/artist.py:72\u001b[0m, in \u001b[0;36mallow_rasterization.<locals>.draw_wrapper\u001b[0;34m(artist, renderer)\u001b[0m\n\u001b[1;32m     69\u001b[0m     \u001b[38;5;28;01mif\u001b[39;00m artist\u001b[38;5;241m.\u001b[39mget_agg_filter() \u001b[38;5;129;01mis\u001b[39;00m \u001b[38;5;129;01mnot\u001b[39;00m \u001b[38;5;28;01mNone\u001b[39;00m:\n\u001b[1;32m     70\u001b[0m         renderer\u001b[38;5;241m.\u001b[39mstart_filter()\n\u001b[0;32m---> 72\u001b[0m     \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[43mdraw\u001b[49m\u001b[43m(\u001b[49m\u001b[43martist\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mrenderer\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m     73\u001b[0m \u001b[38;5;28;01mfinally\u001b[39;00m:\n\u001b[1;32m     74\u001b[0m     \u001b[38;5;28;01mif\u001b[39;00m artist\u001b[38;5;241m.\u001b[39mget_agg_filter() \u001b[38;5;129;01mis\u001b[39;00m \u001b[38;5;129;01mnot\u001b[39;00m \u001b[38;5;28;01mNone\u001b[39;00m:\n",
      "File \u001b[0;32m/usr/local/lib/python3.9/site-packages/matplotlib/axes/_base.py:3070\u001b[0m, in \u001b[0;36m_AxesBase.draw\u001b[0;34m(self, renderer)\u001b[0m\n\u001b[1;32m   3067\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m artists_rasterized:\n\u001b[1;32m   3068\u001b[0m     _draw_rasterized(\u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mfigure, artists_rasterized, renderer)\n\u001b[0;32m-> 3070\u001b[0m \u001b[43mmimage\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43m_draw_list_compositing_images\u001b[49m\u001b[43m(\u001b[49m\n\u001b[1;32m   3071\u001b[0m \u001b[43m    \u001b[49m\u001b[43mrenderer\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;28;43mself\u001b[39;49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43martists\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mfigure\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43msuppressComposite\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m   3073\u001b[0m renderer\u001b[38;5;241m.\u001b[39mclose_group(\u001b[38;5;124m'\u001b[39m\u001b[38;5;124maxes\u001b[39m\u001b[38;5;124m'\u001b[39m)\n\u001b[1;32m   3074\u001b[0m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mstale \u001b[38;5;241m=\u001b[39m \u001b[38;5;28;01mFalse\u001b[39;00m\n",
      "File \u001b[0;32m/usr/local/lib/python3.9/site-packages/matplotlib/image.py:132\u001b[0m, in \u001b[0;36m_draw_list_compositing_images\u001b[0;34m(renderer, parent, artists, suppress_composite)\u001b[0m\n\u001b[1;32m    130\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m not_composite \u001b[38;5;129;01mor\u001b[39;00m \u001b[38;5;129;01mnot\u001b[39;00m has_images:\n\u001b[1;32m    131\u001b[0m     \u001b[38;5;28;01mfor\u001b[39;00m a \u001b[38;5;129;01min\u001b[39;00m artists:\n\u001b[0;32m--> 132\u001b[0m         \u001b[43ma\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mdraw\u001b[49m\u001b[43m(\u001b[49m\u001b[43mrenderer\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m    133\u001b[0m \u001b[38;5;28;01melse\u001b[39;00m:\n\u001b[1;32m    134\u001b[0m     \u001b[38;5;66;03m# Composite any adjacent images together\u001b[39;00m\n\u001b[1;32m    135\u001b[0m     image_group \u001b[38;5;241m=\u001b[39m []\n",
      "File \u001b[0;32m/usr/local/lib/python3.9/site-packages/matplotlib/artist.py:72\u001b[0m, in \u001b[0;36mallow_rasterization.<locals>.draw_wrapper\u001b[0;34m(artist, renderer)\u001b[0m\n\u001b[1;32m     69\u001b[0m     \u001b[38;5;28;01mif\u001b[39;00m artist\u001b[38;5;241m.\u001b[39mget_agg_filter() \u001b[38;5;129;01mis\u001b[39;00m \u001b[38;5;129;01mnot\u001b[39;00m \u001b[38;5;28;01mNone\u001b[39;00m:\n\u001b[1;32m     70\u001b[0m         renderer\u001b[38;5;241m.\u001b[39mstart_filter()\n\u001b[0;32m---> 72\u001b[0m     \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[43mdraw\u001b[49m\u001b[43m(\u001b[49m\u001b[43martist\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mrenderer\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m     73\u001b[0m \u001b[38;5;28;01mfinally\u001b[39;00m:\n\u001b[1;32m     74\u001b[0m     \u001b[38;5;28;01mif\u001b[39;00m artist\u001b[38;5;241m.\u001b[39mget_agg_filter() \u001b[38;5;129;01mis\u001b[39;00m \u001b[38;5;129;01mnot\u001b[39;00m \u001b[38;5;28;01mNone\u001b[39;00m:\n",
      "File \u001b[0;32m/usr/local/lib/python3.9/site-packages/matplotlib/image.py:649\u001b[0m, in \u001b[0;36m_ImageBase.draw\u001b[0;34m(self, renderer, *args, **kwargs)\u001b[0m\n\u001b[1;32m    647\u001b[0m         renderer\u001b[38;5;241m.\u001b[39mdraw_image(gc, l, b, im, trans)\n\u001b[1;32m    648\u001b[0m \u001b[38;5;28;01melse\u001b[39;00m:\n\u001b[0;32m--> 649\u001b[0m     im, l, b, trans \u001b[38;5;241m=\u001b[39m \u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mmake_image\u001b[49m\u001b[43m(\u001b[49m\n\u001b[1;32m    650\u001b[0m \u001b[43m        \u001b[49m\u001b[43mrenderer\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mrenderer\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mget_image_magnification\u001b[49m\u001b[43m(\u001b[49m\u001b[43m)\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m    651\u001b[0m     \u001b[38;5;28;01mif\u001b[39;00m im \u001b[38;5;129;01mis\u001b[39;00m \u001b[38;5;129;01mnot\u001b[39;00m \u001b[38;5;28;01mNone\u001b[39;00m:\n\u001b[1;32m    652\u001b[0m         renderer\u001b[38;5;241m.\u001b[39mdraw_image(gc, l, b, im)\n",
      "File \u001b[0;32m/usr/local/lib/python3.9/site-packages/matplotlib/image.py:939\u001b[0m, in \u001b[0;36mAxesImage.make_image\u001b[0;34m(self, renderer, magnification, unsampled)\u001b[0m\n\u001b[1;32m    936\u001b[0m transformed_bbox \u001b[38;5;241m=\u001b[39m TransformedBbox(bbox, trans)\n\u001b[1;32m    937\u001b[0m clip \u001b[38;5;241m=\u001b[39m ((\u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mget_clip_box() \u001b[38;5;129;01mor\u001b[39;00m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39maxes\u001b[38;5;241m.\u001b[39mbbox) \u001b[38;5;28;01mif\u001b[39;00m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mget_clip_on()\n\u001b[1;32m    938\u001b[0m         \u001b[38;5;28;01melse\u001b[39;00m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mfigure\u001b[38;5;241m.\u001b[39mbbox)\n\u001b[0;32m--> 939\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43m_make_image\u001b[49m\u001b[43m(\u001b[49m\u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43m_A\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mbbox\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mtransformed_bbox\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mclip\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m    940\u001b[0m \u001b[43m                        \u001b[49m\u001b[43mmagnification\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43munsampled\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43munsampled\u001b[49m\u001b[43m)\u001b[49m\n",
      "File \u001b[0;32m/usr/local/lib/python3.9/site-packages/matplotlib/image.py:526\u001b[0m, in \u001b[0;36m_ImageBase._make_image\u001b[0;34m(self, A, in_bbox, out_bbox, clip_bbox, magnification, unsampled, round_to_pixel_border)\u001b[0m\n\u001b[1;32m    521\u001b[0m mask \u001b[38;5;241m=\u001b[39m (np\u001b[38;5;241m.\u001b[39mwhere(A\u001b[38;5;241m.\u001b[39mmask, np\u001b[38;5;241m.\u001b[39mfloat32(np\u001b[38;5;241m.\u001b[39mnan), np\u001b[38;5;241m.\u001b[39mfloat32(\u001b[38;5;241m1\u001b[39m))\n\u001b[1;32m    522\u001b[0m         \u001b[38;5;28;01mif\u001b[39;00m A\u001b[38;5;241m.\u001b[39mmask\u001b[38;5;241m.\u001b[39mshape \u001b[38;5;241m==\u001b[39m A\u001b[38;5;241m.\u001b[39mshape  \u001b[38;5;66;03m# nontrivial mask\u001b[39;00m\n\u001b[1;32m    523\u001b[0m         \u001b[38;5;28;01melse\u001b[39;00m np\u001b[38;5;241m.\u001b[39mones_like(A, np\u001b[38;5;241m.\u001b[39mfloat32))\n\u001b[1;32m    524\u001b[0m \u001b[38;5;66;03m# we always have to interpolate the mask to account for\u001b[39;00m\n\u001b[1;32m    525\u001b[0m \u001b[38;5;66;03m# non-affine transformations\u001b[39;00m\n\u001b[0;32m--> 526\u001b[0m out_alpha \u001b[38;5;241m=\u001b[39m \u001b[43m_resample\u001b[49m\u001b[43m(\u001b[49m\u001b[38;5;28;43mself\u001b[39;49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mmask\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mout_shape\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mt\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mresample\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[38;5;28;43;01mTrue\u001b[39;49;00m\u001b[43m)\u001b[49m\n\u001b[1;32m    527\u001b[0m \u001b[38;5;28;01mdel\u001b[39;00m mask  \u001b[38;5;66;03m# Make sure we don't use mask anymore!\u001b[39;00m\n\u001b[1;32m    528\u001b[0m \u001b[38;5;66;03m# Agg updates out_alpha in place.  If the pixel has no image\u001b[39;00m\n\u001b[1;32m    529\u001b[0m \u001b[38;5;66;03m# data it will not be updated (and still be 0 as we initialized\u001b[39;00m\n\u001b[1;32m    530\u001b[0m \u001b[38;5;66;03m# it), if input data that would go into that output pixel than\u001b[39;00m\n\u001b[1;32m    531\u001b[0m \u001b[38;5;66;03m# it will be `nan`, if all the input data for a pixel is good\u001b[39;00m\n\u001b[1;32m    532\u001b[0m \u001b[38;5;66;03m# it will be 1, and if there is _some_ good data in that output\u001b[39;00m\n\u001b[1;32m    533\u001b[0m \u001b[38;5;66;03m# pixel it will be between [0, 1] (such as a rotated image).\u001b[39;00m\n",
      "File \u001b[0;32m/usr/local/lib/python3.9/site-packages/matplotlib/image.py:208\u001b[0m, in \u001b[0;36m_resample\u001b[0;34m(image_obj, data, out_shape, transform, resample, alpha)\u001b[0m\n\u001b[1;32m    206\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m resample \u001b[38;5;129;01mis\u001b[39;00m \u001b[38;5;28;01mNone\u001b[39;00m:\n\u001b[1;32m    207\u001b[0m     resample \u001b[38;5;241m=\u001b[39m image_obj\u001b[38;5;241m.\u001b[39mget_resample()\n\u001b[0;32m--> 208\u001b[0m \u001b[43m_image\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mresample\u001b[49m\u001b[43m(\u001b[49m\u001b[43mdata\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mout\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mtransform\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m    209\u001b[0m \u001b[43m                \u001b[49m\u001b[43m_interpd_\u001b[49m\u001b[43m[\u001b[49m\u001b[43minterpolation\u001b[49m\u001b[43m]\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m    210\u001b[0m \u001b[43m                \u001b[49m\u001b[43mresample\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m    211\u001b[0m \u001b[43m                \u001b[49m\u001b[43malpha\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m    212\u001b[0m \u001b[43m                \u001b[49m\u001b[43mimage_obj\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mget_filternorm\u001b[49m\u001b[43m(\u001b[49m\u001b[43m)\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m    213\u001b[0m \u001b[43m                \u001b[49m\u001b[43mimage_obj\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mget_filterrad\u001b[49m\u001b[43m(\u001b[49m\u001b[43m)\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m    214\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m out\n",
      "\u001b[0;31mKeyboardInterrupt\u001b[0m: "
     ]
    }
   ],
   "source": [
    "def plotlogimag(imag_log, boundary, title):\n",
    "    lon_min = boundary[0]\n",
    "    lon_max = boundary[1]\n",
    "    lat_min = boundary[2]\n",
    "    lat_max = boundary[3]\n",
    "\n",
    "    # Plot the data with geographic coordinates\n",
    "    fig, ax = plt.subplots(figsize=(20, 11), dpi = 200)\n",
    "    cax = ax.imshow(imag_log, extent=[lon_min, lon_max, lat_min, lat_max], origin='upper', cmap = 'jet')\n",
    "\n",
    "    ax.set_xlabel('Longitude')\n",
    "    ax.set_ylabel('Latitude')\n",
    "    ax.set_title(f'Artificial Sky Brightness of {title}')\n",
    "    ax.set_xticks(np.arange(lon_min, lon_max, 10))\n",
    "    ax.set_yticks(np.arange(lat_min, lat_max, 10))\n",
    "    ax.grid()\n",
    "    fig.colorbar(cax, ax=ax, label='Log Brightness')\n",
    "    fig.savefig('WorldPop2015/'+title, dpi = 200)\n",
    "plotlogimag(World, bWorld, 'the World')\n",
    "plotlogimag(Europe, bEurope, 'Europe')\n",
    "plotlogimag(Africa, bAfrica, 'Africa')\n",
    "plotlogimag(South_America, bSouth_America, 'South America')\n",
    "plotlogimag(North_America, bNorth_America, 'North America')\n",
    "plotlogimag(Australia, bAustralia, 'Australia')\n",
    "plotlogimag(Asia, bAsia, 'Asia')\n",
    "\n",
    "# "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c023e33f",
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_H(data_log):\n",
    "    std = np.std(data_log)\n",
    "    N = data_log.shape[0] * data_log.shape[1]\n",
    "    H = 1.06 * std * N ** (-1/5)\n",
    "    return H\n",
    "h = get_H(World)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "60589964",
   "metadata": {},
   "outputs": [],
   "source": [
    "# def Epinechkov_kernel(x, y):\n",
    "#     \"\"\"\n",
    "#     Compute the Epinechnikov kernel for a pair of points x and y\n",
    "#     \"\"\"\n",
    "#     # Compute the distance between the points\n",
    "#     d = np.sqrt((x[0] - y[0])**2 + (x[1] - y[1])**2)\n",
    "#     # Compute the kernel\n",
    "#     if d < h:\n",
    "#         return 3/4 * (1 - d**2 / h**2)\n",
    "#     else:\n",
    "#         return 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0585a237",
   "metadata": {},
   "outputs": [],
   "source": [
    "# def Gaussian_kernel(x, y, sigma=1):\n",
    "#     \"\"\"\n",
    "#     Compute the Gaussian kernel for a pair of points x and y\n",
    "#     \"\"\"\n",
    "#     # Compute the distance between the points\n",
    "#     d = np.sqrt((x[0] - y[0])**2 + (x[1] - y[1])**2)\n",
    "#     # Compute the kernel\n",
    "#     return np.exp(-d**2 / (2 * sigma**2))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "703c9eea",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[4.2614938e-05 4.2823671e-05 4.3040294e-05 ... 0.0000000e+00\n",
      "  0.0000000e+00 0.0000000e+00]\n",
      " [3.9259321e-05 3.9470444e-05 3.9689068e-05 ... 0.0000000e+00\n",
      "  0.0000000e+00 0.0000000e+00]\n",
      " [3.5898120e-05 3.6049503e-05 3.6437134e-05 ... 0.0000000e+00\n",
      "  0.0000000e+00 0.0000000e+00]\n",
      " ...\n",
      " [0.0000000e+00 0.0000000e+00 0.0000000e+00 ... 0.0000000e+00\n",
      "  0.0000000e+00 0.0000000e+00]\n",
      " [0.0000000e+00 0.0000000e+00 0.0000000e+00 ... 0.0000000e+00\n",
      "  0.0000000e+00 0.0000000e+00]\n",
      " [0.0000000e+00 0.0000000e+00 0.0000000e+00 ... 0.0000000e+00\n",
      "  0.0000000e+00 0.0000000e+00]]\n"
     ]
    }
   ],
   "source": [
    "def KDEcore(data):\n",
    "    kde = KernelDensity(bandwidth = h, kernel='epanechnikov')\n",
    "    kde.fit(data)\n",
    "    \n",
    "    \n",
    "    \n",
    "    longitudes = np.linspace(-180, 180, 43200)\n",
    "    latitudes = np.linspace(-90, 90, 17406)\n",
    "\n",
    "\n",
    "# print(imag[2000:2024])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5f2c632d",
   "metadata": {},
   "outputs": [],
   "source": [
    "def twoDimensionalEpanechnikovKDE(data):\n",
    "    X = data[['lon', 'lat']].values\n",
    "    kde = KernelDensity(bandwidth=0.01, kernel='epanechnikov')\n",
    "    kde.fit(X)\n",
    "    longitudes = np.linspace(data['lon'].min(), data['lon'].max(), 100)\n",
    "    latitudes = np.linspace(data['lat'].min(), data['lat'].max(), 100)\n",
    "    total_area = (longitudes.max() - longitudes.min()) * (latitudes.max() - latitudes.min())\n",
    "\n",
    "    X_grid, Y_grid = np.meshgrid(longitudes, latitudes)\n",
    "    grid_points = np.vstack([X_grid.ravel(), Y_grid.ravel()]).T\n",
    "\n",
    "    # Calculate the log density (log PDF) values for the grid\n",
    "    log_pdf_values = kde.score_samples(grid_points)\n",
    "    PDF_values = np.exp(log_pdf_values).reshape(X_grid.shape)\n",
    "\n",
    "    fig, ax = plt.subplots(figsize=(12, 9))\n",
    "    contourf_plot = ax.contourf(X_grid, Y_grid, PDF_values, cmap = 'cool')\n",
    "#     ax.scatter(data['lon'], data['lat'], s=5, color = COLOUR[2], alpha=0.5, label='Crash Incidents')\n",
    "    fig.colorbar(contourf_plot, ax=ax, label='KDE Value')\n",
    "    ax.set_xlabel('Longitude')\n",
    "    ax.set_ylabel('Latitude')\n",
    "    ax.set_title('Map of Epanechnikov KDE')\n",
    "\n",
    "    ## box volume proportion\n",
    "    box_longitude = np.array([-78.76, -78.72])\n",
    "    box_latitude = np.array([35.74, 35.78])\n",
    "    box_area = (box_longitude[1] - box_longitude[0]) * (box_latitude[1] - box_latitude[0])\n",
    "    \n",
    "    # plot\n",
    "    ax.plot([box_longitude[0], box_longitude[1], box_longitude[1], box_longitude[0], box_longitude[0]], \n",
    "            [box_latitude[0], box_latitude[0], box_latitude[1], box_latitude[1], box_latitude[0]], \n",
    "            color=COLOUR[2], linestyle='--', label='Region of Interest')\n",
    "    proportion = getProportion(PDF_values, longitudes, latitudes, box_longitude, box_latitude)\n",
    "    print(f\"Proportion of the region within the total area: {proportion:.4f}\")\n",
    "    ax.legend()\n",
    "\n",
    "    # Add text to the plot\n",
    "    d = {'prob. proportion[box/total]': f'{proportion:.2%}', \n",
    "         'area proportion[box/total]': f'{box_area / total_area:.2%}',\n",
    "         'PDF total(normalisation)' : f'{normalisation(PDF_values, longitudes, latitudes):.4f}',\n",
    "         'max PDF value' : f'{PDF_values.max():.4f}',\n",
    "         }\n",
    "    add_text_to_ax(0.05, 0.30, nice_string_output(d), ax, fontsize=12)\n",
    "\n",
    "\n",
    "    fig_mag, ax_mag = plt.subplots(figsize=(12, 9))\n",
    "    contourf_plot_mag = ax_mag.contourf(X_grid, Y_grid, PDF_values, cmap = 'cool')\n",
    "    fig_mag.colorbar(contourf_plot_mag, ax=ax_mag, label='KDE Value')\n",
    "    ax_mag.scatter(data['lon'], data['lat'], s=5, color = COLOUR[1], alpha=0.5, label='Crash Incidents')\n",
    "    ax_mag.set_xlabel('Longitude')\n",
    "    ax_mag.set_ylabel('Latitude')\n",
    "    ax_mag.set_title('Map of Epanechnikov KDE(magnified)')\n",
    "    ax_mag.set_xlim(-78.9, -78.7)\n",
    "    ax_mag.set_ylim(35.725, 35.825)\n",
    "    ax_mag.plot([box_longitude[0], box_longitude[1], box_longitude[1], box_longitude[0], box_longitude[0]], \n",
    "            [box_latitude[0], box_latitude[0], box_latitude[1], box_latitude[1], box_latitude[0]], \n",
    "            color=COLOUR[2], linestyle='--', label='Region of Interest')\n",
    "    ax_mag.legend(loc = 'upper right')\n",
    "    add_text_to_ax(0.05, 0.32, nice_string_output(d), ax_mag, fontsize=12)\n",
    "    \n",
    "twoDimensionalEpanechnikovKDE(data)"
   ]
  }
 ],
 "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.18"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}