Land_Clusters.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "precise-limitation",
   "metadata": {},
   "outputs": [],
   "source": [
    "#In this program, we calculate the largest Land cluster on Earth in the presence of different levels of water\n",
    "#Libraries\n",
    "from scipy.ndimage.filters import gaussian_filter\n",
    "import matplotlib.pyplot as plt\n",
    "from netCDF4 import Dataset\n",
    "import numpy as np\n",
    "import os\n",
    "from tqdm import tqdm\n",
    "from scipy.ndimage import measurements\n",
    "import time\n",
    "import sympy\n",
    "import concurrent.futures"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "subject-customs",
   "metadata": {},
   "outputs": [],
   "source": [
    "#This function gives the different resolution of the earth map. we use it to extract the orginal resolution.\n",
    "#the code has been written by\n",
    "def Etopo(lon_area, lat_area, resolution):\n",
    "    ### Input\n",
    "    # resolution: resolution of topography for both of longitude and latitude [deg]\n",
    "    # (Original resolution is 0.0167 deg)\n",
    "    # lon_area and lat_area: the region of the map which you want like [100, 130], [20, 25]\n",
    "    ###\n",
    "    ### Output\n",
    "    # Mesh type longitude, latitude, and topography data\n",
    "    ###\n",
    "\n",
    "    # Read NetCDF data\n",
    "    data = Dataset(\"ETOPO1_Bed_c_gdal.grd\", \"r\")\n",
    "    # Get data\n",
    "    lon_range = data.variables['x_range'][:]\n",
    "    lat_range = data.variables['y_range'][:]\n",
    "    topo_range = data.variables['z_range'][:]\n",
    "    spacing = data.variables['spacing'][:]\n",
    "    dimension = data.variables['dimension'][:]\n",
    "    z = data.variables['z'][:]\n",
    "    lon_num = dimension[0]\n",
    "    lat_num = dimension[1]\n",
    "  \n",
    "    # Prepare array\n",
    "    lon_input = np.zeros(lon_num); lat_input = np.zeros(lat_num)\n",
    "    for i in range(lon_num):\n",
    "        lon_input[i] = lon_range[0] + i * spacing[0]\n",
    "    for i in range(lat_num):\n",
    "        lat_input[i] = lat_range[0] + i * spacing[1]\n",
    "\n",
    "    # Create 2D array\n",
    "    lon, lat = np.meshgrid(lon_input, lat_input)\n",
    "\n",
    "    # Convert 2D array from 1D array for z value\n",
    "    topo = np.reshape(z, (lat_num, lon_num))\n",
    "  \n",
    "    # Skip the data for resolution\n",
    "    if ((resolution < spacing[0]) | (resolution < spacing[1])):\n",
    "        print('Set the highest resolution')\n",
    "    else:\n",
    "        skip = int(resolution/spacing[0])\n",
    "        lon = lon[::skip,::skip]\n",
    "        lat = lat[::skip,::skip]\n",
    "        topo = topo[::skip,::skip]\n",
    "\n",
    "    topo = topo[::-1]\n",
    "  \n",
    "    # Select the range of map\n",
    "    range1 = np.where((lon>=lon_area[0]) & (lon<=lon_area[1]))\n",
    "    lon = lon[range1]; lat = lat[range1]; topo = topo[range1]\n",
    "    range2 = np.where((lat>=lat_area[0]) & (lat<=lat_area[1]))\n",
    "    lon = lon[range2]; lat = lat[range2]; topo = topo[range2]\n",
    "  \n",
    "    # Convert 2D again\n",
    "    lon_num = len(np.unique(lon))\n",
    "    lat_num = len(np.unique(lat))\n",
    "    lon = np.reshape(lon, (lat_num, lon_num))\n",
    "    lat = np.reshape(lat, (lat_num, lon_num))\n",
    "    topo = np.reshape(topo, (lat_num, lon_num))\n",
    "\n",
    "    return lon, lat, topo"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "neural-ordering",
   "metadata": {},
   "outputs": [],
   "source": [
    "# convert degrees to radians\n",
    "def degree2radians(degree):\n",
    "    return degree*np.pi/180"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "foreign-ferry",
   "metadata": {},
   "outputs": [],
   "source": [
    "#defining the orginal resoluton of the data\n",
    "resolution = 0.0167\n",
    "# resolution = 1\n",
    "\n",
    "lon_area = [-180., 180.]\n",
    "lat_area = [-90., 90.]\n",
    "# Get mesh-shape topography data\n",
    "lon_topo, lat_topo, topo = Etopo(lon_area, lat_area, resolution)\n",
    "lon_topo = degree2radians(lon_topo)\n",
    "lat_topo = degree2radians(lat_topo)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "saved-journal",
   "metadata": {},
   "outputs": [],
   "source": [
    "#The Course graining part. we define the of the course graining window with N and course grain it with thw function rebin\n",
    "N = 10\n",
    "def rebin(a, shape):\n",
    "    sh = shape[0],a.shape[0]//shape[0],shape[1],a.shape[1]//shape[1]\n",
    "    return a.reshape(sh).mean(-1).mean(1)\n",
    "\n",
    "Width = int(10800/N)\n",
    "Length = int(21600/N)\n",
    "topo = rebin(topo,(Width, Length))\n",
    "lat_topo = rebin(lat_topo,(Width, Length))\n",
    "lon_topo = rebin(lon_topo,(Width, Length))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "honest-meter",
   "metadata": {},
   "outputs": [],
   "source": [
    "#we extract all the elevations and sort then in an ascending order\n",
    "indexes = np.argsort(topo.flatten())\n",
    "hieght = np.sort(topo.flatten())\n",
    "#defining the parallelization variables\n",
    "l = len(hieght)\n",
    "d = sympy.divisors(l)\n",
    "num_iter = d[int(len(d)/2)]\n",
    "num_proc = d[int(len(d)/2) - 1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "protective-isaac",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Land Clusters\n",
    "\n",
    "reversed_indexes = indexes[::-1]\n",
    "total_landmass = np.zeros(len(hieght))\n",
    "big_cluster = np.zeros(len(hieght))\n",
    "total_lat = np.sum(np.abs(np.cos(lat_topo)))\n",
    "all_pos = np.shape(lat_topo)[0]*np.shape(lat_topo)[1]\n",
    "\n",
    "\n",
    "squares1 = np.zeros((num_proc, num_iter))\n",
    "squares2 = np.zeros((num_proc, num_iter))\n",
    "\n",
    "\n",
    "def land_function(ii):\n",
    "    total_landmass = np.zeros(num_iter)\n",
    "    big_cluster = np.zeros(num_iter)\n",
    "    iso_hight = np.zeros(num_iter)\n",
    "    iterr = 0\n",
    "    for i in range(ii*num_iter, ii*num_iter + num_iter):\n",
    "        sphere = 0\n",
    "        new_topo=np.zeros(hieght.shape, bool)\n",
    "        new_topo[reversed_indexes[:(i+1)]] = True\n",
    "        new_topo = new_topo.reshape((Width, Length))\n",
    "        pp = (i+1)/(all_pos)\n",
    "        lw, num = measurements.label(new_topo)\n",
    "        for l in range(len(lw)):\n",
    "            if(lw[l][0] != 0 and lw[l][-1] != 0 and lw[l][-1] != lw[l][0]):\n",
    "                lw[lw == lw[l][0]] = lw[l][-1]\n",
    "\n",
    "        biggest_size = 0\n",
    "        elements = np.unique(lw)\n",
    "        elements = elements[elements != 0]\n",
    "        unique_clustersize = np.zeros(len(elements))\n",
    "        iterat = 0\n",
    "        for j in elements:\n",
    "            unique_clustersize[iterat] = np.sum(lw == j)\n",
    "            iterat +=1\n",
    "\n",
    "        elements = elements[unique_clustersize.argsort()]\n",
    "\n",
    "        a = 0\n",
    "        if len(elements)>2:\n",
    "            for n in elements[[-1,-2,-3]]:\n",
    "                mask = (lw == n)\n",
    "                a = np.sum(np.abs(np.cos((mask*lat_topo)[(mask*lat_topo) != 0])))\n",
    "                if a>biggest_size:\n",
    "                    biggest_size = a\n",
    "\n",
    "        elif len(elements) == 2:\n",
    "            for n in elements[[-1,-2]]:\n",
    "                mask = (lw == n)\n",
    "                a = np.sum(np.abs(np.cos((mask*lat_topo)[(mask*lat_topo) != 0])))\n",
    "                if a>biggest_size:\n",
    "                    biggest_size = a\n",
    "\n",
    "        else:\n",
    "            mask = (lw == elements[0])\n",
    "            biggest_size = np.sum(np.abs(np.cos((mask*lat_topo)[(mask*lat_topo) != 0])))\n",
    "\n",
    "\n",
    "        total_landmass[iterr] = pp\n",
    "        big_cluster[iterr] = float(biggest_size / total_lat)\n",
    "        iterr += 1    \n",
    "    \n",
    "    \n",
    "        \n",
    "    return ii, big_cluster, total_landmass, iso_hight\n",
    "\n",
    "with concurrent.futures.ProcessPoolExecutor() as executor: \n",
    "    for row, result1, result2 in executor.map(land_function, range(num_proc)):\n",
    "        print(row)\n",
    "        squares1[row] = result1\n",
    "        squares2[row] = result2\n",
    "\n",
    "\n",
    "np.save(\"/home/complex/c++/Earth/big_cluster_land_60.npy\", squares1.flatten())\n",
    "np.save(\"/home/complex/c++/Earth/total_landmass_land_60.npy\", squares2.flatten())"
   ]
  }
 ],
 "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.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}