Equal_Dist_Grid_Points 06 R0 all.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "#\n",
    "# This program defines a grid of points on the surface of a geoid.\n",
    "# From a given point (P0) a number of \"columns\" and \"rows\" are defined.\n",
    "# The size of the grid is specified by introducing:\n",
    "# - the number of columns to the \"right\" and \"left\" of P0, \n",
    "# - the number of rows \"below\" and \"above\" P0.\n",
    "# The points are created so that the geographic distance between them (row-row and col-col) is a given value.\n",
    "# All points in the same row than P0 have the same latitude as P0.\n",
    "# All points in the same column than P0 share with it the longitude.\n",
    "# The distances are calculated with a geodesic formula.\n",
    "# Units: \n",
    "# - distance: km.\n",
    "# - coordinates: decimal degrees.\n",
    "# The coordinates are used as (lon, lat) NOT as (lat, lon).\n",
    "# The results are exporte as:\n",
    "# - a csv file with the coordinates of the points,\n",
    "# - a geojson file with the polygons that can be defined from each 4-set of adjacent points.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Version log.\n",
    "# \n",
    "# R0 (20200401)\n",
    "# All parts tested and put together.\n",
    "# Seems to work well.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Files and directories.\n",
    "\n",
    "# Directories:\n",
    "RootDir = 'D:/0 DOWN/00 PY RG/GIS/GRIDS/'\n",
    "\n",
    "# I/O files:\n",
    "FileNameOut1 = RootDir + 'Data_out 06 pnt R0.csv'\n",
    "FileNameOut2 = RootDir + 'Data_out 06 pol R0.geojson'\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import modules.\n",
    "from math import radians, degrees, cos, copysign\n",
    "from geopy.distance import geodesic, lonlat\n",
    "import csv, json\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Functions.\n",
    "def f_dlon(lat, dist, R = 6371):\n",
    "    \"\"\"\n",
    "    Calculate the increase in longitude required to cover a given distance along a parallel.\n",
    "    Asumes a sphere.\n",
    "    Data in degrees and km.\n",
    "    \"\"\"\n",
    "    return degrees(dist / (R * cos(radians(lat))))\n",
    "\n",
    "def f_dlat(dist, R = 6371):\n",
    "    \"\"\"\n",
    "    Calculate the increase in latitude required to cover a given distance along a meridian.\n",
    "    Asumes a sphere.\n",
    "    Data in degrees and km.\n",
    "    \"\"\"\n",
    "    return degrees(dist / R)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Input data.\n",
    "# Starting point:\n",
    "P0 = [6., 45.]\n",
    "# Geographic distance between points of the grid (x and y-wise):\n",
    "D0 = 100.\n",
    "# Grid size, measured as number or rows/columns :\n",
    "nx_1 = 16 # decreasing longitude, columns - excluding the reference \n",
    "nx_2 = 21 # increasing longitude, columns - excluding the reference\n",
    "ny_1 = 11 # decreasing latitude, rows - excluding the reference\n",
    "ny_2 = 18 # increasing latitude, rows - excluding the reference\n",
    "\n",
    "# Iteration parameters:\n",
    "# Accepted absolute error:\n",
    "abs_err = D0 / 1000.\n",
    "# Maximum number of iterations:\n",
    "iter_max = 1000\n",
    "\n",
    "# Associated variables:\n",
    "# Grid size (total number of rows/columns):\n",
    "NX = nx_1 + nx_2 + 1 # longitude\n",
    "NY = ny_1 + ny_2 + 1 # latitude\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Check suitability of the data.\n",
    "# For extreme latitudes:\n",
    "lat_min, lat_max = -80., 80. # Arbitrarily close to problematic latitudes\n",
    "lat_1 = P0[1] - f_dlat((ny_1)*D0)\n",
    "lat_2 = P0[1] + f_dlat((ny_2)*D0)\n",
    "if (lat_1 < lat_min) or (lat_2 > lat_max):\n",
    "    print (' lat_min = {0:4.2f}, lat_max = {1:4.2f}'.format(lat_1, lat_2))\n",
    "    raise Exception (' Warning: extreme latitudes exceeded. Program stopped.')\n",
    "\n",
    "# For extreme longitudes:\n",
    "lat_aux = max(abs(lat_1), abs(lat_2))\n",
    "lon_aux = f_dlon(lat_aux, (NX-1)*D0)\n",
    "lon_max = 180. # Arbitrarily chosen.\n",
    "if (lon_aux > lon_max):\n",
    "    print (' lon_max = {0:4.2f}'.format(lon_aux))\n",
    "    raise Exception (' Warning: extreme longitude exceeded. Program stopped.')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Prepare the container of the points, a list of lists.\n",
    "aux = float('nan')\n",
    "\n",
    "# The access is first from x and then from y!\n",
    "data = []\n",
    "for j in range (0, NX, 1):\n",
    "    aux_j = []\n",
    "    for i in range (0, NY, 1):\n",
    "        aux_j.append([aux, aux])\n",
    "    data.append(aux_j)\n",
    "    \n",
    "# Enter P00:\n",
    "data[nx_1 ][ny_1] = P0\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Generate the central row, at the middle of the grid.\n",
    "# Right side:\n",
    "for x in range(nx_1, NX-1, 1):\n",
    "    # P1 is defined:\n",
    "    P1      = data[x][ny_1]\n",
    "    # P2 is on the same parallel than P1:\n",
    "    P2      = [P1[0] + f_dlon(P1[1], D0), P1[1]]\n",
    "    D12     = geodesic(lonlat(*P1), lonlat(*P2)).km\n",
    "    D12_err = D12 - D0\n",
    "    iter_count   = 0\n",
    "    while (abs(D12_err) > abs_err) and (iter_count <= iter_max):\n",
    "        lon_aux = f_dlon(P2[1], abs(D12_err/2))\n",
    "        P2      = [P2[0] - copysign(lon_aux, D12_err), P2[1]]\n",
    "        D12     = geodesic(lonlat(*P1), lonlat(*P2)).km\n",
    "        D12_err = D12 - D0\n",
    "        iter_count  += 1\n",
    "        if (iter_count == iter_max):\n",
    "            print (\" Warning at P2: iter_max exceeded. D_err = \", abs(D12_err)) \n",
    "\n",
    "    # Save to data_list:\n",
    "    data[x+1][ny_1] = P2\n",
    "\n",
    "# Left side:\n",
    "for x in range(nx_1, 0, -1):\n",
    "    # P2 is defined:\n",
    "    P2      = data[x][ny_1]\n",
    "    # P1 is on the same parallel than P2:\n",
    "    P1      = [P2[0] - f_dlon(P2[1], D0), P2[1]]\n",
    "    D12     = geodesic(lonlat(*P1), lonlat(*P2)).km\n",
    "    D12_err = D12 - D0\n",
    "    iter_count   = 0\n",
    "    while (abs(D12_err) > abs_err) and (iter_count < iter_max):\n",
    "        lon_aux = f_dlon(P1[1], abs(D12_err/2))\n",
    "        P1      = [P1[0] + copysign(lon_aux, D12_err), P1[1]]\n",
    "        D12     = geodesic(lonlat(*P1), lonlat(*P2)).km\n",
    "        D12_err = D12 - D0\n",
    "        iter_count  += 1\n",
    "        if (iter_count == iter_max):\n",
    "            print (\" Warning at P1: iter_max exceeded. D_err = \", abs(D12_err))\n",
    "\n",
    "    # Save to data_list:\n",
    "    data[x-1][ny_1] = P1\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Generate the central column, at the center of the grid.\n",
    "# Upper side:\n",
    "for y in range(ny_1, NY-1, 1):\n",
    "    # P1 is defined:\n",
    "    P1      = data[nx_1][y]\n",
    "    # P4 is on the same meridian than P1:\n",
    "    P4      = [P1[0], P1[1] + f_dlat(D0)]\n",
    "    D14     = geodesic(lonlat(*P1), lonlat(*P4)).km\n",
    "    D14_err = D14 - D0\n",
    "    iter_count   = 0\n",
    "    while (abs(D14_err) > abs_err) and (iter_count <= iter_max):\n",
    "        lat_aux = f_dlat(abs(D14_err/2))\n",
    "        P4      = [P4[0], P4[1] - copysign(lat_aux, D14_err)]\n",
    "        D14     = geodesic(lonlat(*P1), lonlat(*P4)).km\n",
    "        D14_err = D14 - D0\n",
    "        iter_count  += 1\n",
    "        if (iter_count == iter_max):\n",
    "            print (\" Warning at P4: iter_max exceeded. D_err = \", abs(D14_err))\n",
    "    \n",
    "    # Save to data_list:\n",
    "    data[nx_1][y+1] = P4\n",
    "\n",
    "# Lower side:\n",
    "for y in range(ny_1, 0, -1):\n",
    "    # P4 is defined:\n",
    "    P4      = data[nx_1][y]\n",
    "    # P1 is on the same meridian than P4:\n",
    "    P1      = [P4[0], P4[1] - f_dlat(D0)]\n",
    "    D14     = geodesic(lonlat(*P1), lonlat(*P4)).km\n",
    "    D14_err = D14 - D0\n",
    "    iter_count   = 0\n",
    "    while (abs(D14_err) > abs_err) and (iter_count <= iter_max):\n",
    "        lat_aux = f_dlat(abs(D14_err/2))\n",
    "        P1      = [P1[0], P1[1] + copysign(lat_aux, D14_err)]\n",
    "        D14     = geodesic(lonlat(*P1), lonlat(*P4)).km\n",
    "        D14_err = D14 - D0\n",
    "        iter_count  += 1\n",
    "        if (iter_count == iter_max):\n",
    "            print (\" Warning at P1: iter_max exceeded. D_err = \", abs(D14_err))\n",
    "    \n",
    "    # Save to data_list:\n",
    "    data[nx_1][y-1] = P1\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Complete the grid from the central row and column.\n",
    "# First quadrant:\n",
    "for x in range(nx_1, NX-1, 1):\n",
    "    for y in range(ny_1, NY-1, 1):\n",
    "        # P1, P2 and P4 are defined:\n",
    "        P1 = data[x  ][y  ]\n",
    "        P2 = data[x+1][y  ]\n",
    "        P4 = data[x  ][y+1]\n",
    "        # P3 is approx. on the same longitude that P2, and approx. on the same latitude than P4:\n",
    "        P3      = [P2[0], P4[1]]\n",
    "        D23     = geodesic(lonlat(*P2), lonlat(*P3)).km\n",
    "        D34     = geodesic(lonlat(*P3), lonlat(*P4)).km\n",
    "        D23_err = D23 - D0\n",
    "        D34_err = D34 - D0\n",
    "        D_err   = max(abs(D23_err), abs(D34_err))\n",
    "        iter_count   = 0\n",
    "        while (abs(D_err) > abs_err) and (iter_count <= iter_max):\n",
    "            lat_aux = f_dlat(abs(D23_err/2))\n",
    "            lon_aux = f_dlon(P3[1] - copysign(lat_aux, D23_err), abs(D34_err/2))\n",
    "            P3      = [P3[0] - copysign(lon_aux, D34_err), P3[1] - copysign(lat_aux, D23_err)]\n",
    "            D23     = geodesic(lonlat(*P2), lonlat(*P3)).km\n",
    "            D34     = geodesic(lonlat(*P3), lonlat(*P4)).km\n",
    "            D23_err = D23 - D0\n",
    "            D34_err = D34 - D0\n",
    "            D_err   = max(abs(D23_err), abs(D34_err))\n",
    "            iter_count  += 1\n",
    "            if (iter_count == iter_max):\n",
    "                print (\" Warning at P3: iter_max exceeded. D_err = \", D_err) \n",
    "    \n",
    "        # Save to data_list:\n",
    "        data[x+1][y+1] = P3\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Second quadrant:\n",
    "for x in range(nx_1, 0, -1):\n",
    "    for y in range(ny_1, NY-1, 1):\n",
    "        # P1, P2 and P3 are defined:\n",
    "        P1 = data[x-1][y  ]\n",
    "        P2 = data[x  ][y  ]\n",
    "        P3 = data[x  ][y+1]\n",
    "        # P4 is approx. on the same longitude that P1, and approx. on the same latitude than P3:\n",
    "        P4      = [P1[0], P3[1]]\n",
    "        D14     = geodesic(lonlat(*P1), lonlat(*P4)).km\n",
    "        D34     = geodesic(lonlat(*P3), lonlat(*P4)).km\n",
    "        D14_err = D14 - D0\n",
    "        D34_err = D34 - D0\n",
    "        D_err   = max(abs(D14_err), abs(D34_err))\n",
    "        iter_count   = 0\n",
    "        while (abs(D_err) > abs_err) and (iter_count <= iter_max):\n",
    "            lat_aux = f_dlat(abs(D14_err/2))\n",
    "            lon_aux = f_dlon(P4[1] - copysign(lat_aux, D14_err), abs(D34_err/2))\n",
    "            P4      = [P4[0] + copysign(lon_aux, D34_err), P4[1] - copysign(lat_aux, D14_err)]\n",
    "            D14     = geodesic(lonlat(*P1), lonlat(*P4)).km\n",
    "            D34     = geodesic(lonlat(*P3), lonlat(*P4)).km\n",
    "            D14_err = D14 - D0\n",
    "            D34_err = D34 - D0\n",
    "            D_err   = max(abs(D14_err), abs(D34_err))\n",
    "            iter_count  += 1\n",
    "            if (iter_count == iter_max):\n",
    "                print (\" Warning at P4: iter_max exceeded. D_err = \", D_err) \n",
    "    \n",
    "        # Save to data_list:\n",
    "        data[x-1][y+1] = P4"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Third quadrant:\n",
    "for x in range(nx_1, 0, -1):\n",
    "    for y in range(ny_1, 0, -1):\n",
    "        # P2, P3 and P4 are defined:\n",
    "        P2 = data[x  ][y-1]\n",
    "        P3 = data[x  ][y  ]\n",
    "        P4 = data[x-1][y  ]\n",
    "        # P1 is approx. on the same longitude that P4, and approx. on the same latitude than P2:\n",
    "        P1      = [P4[0], P2[1]]\n",
    "        D12     = geodesic(lonlat(*P1), lonlat(*P2)).km\n",
    "        D14     = geodesic(lonlat(*P1), lonlat(*P4)).km\n",
    "        D12_err = D12 - D0\n",
    "        D14_err = D14 - D0\n",
    "        D_err   = max(abs(D12_err), abs(D14_err))\n",
    "        iter_count   = 0\n",
    "        while (abs(D_err) > abs_err) and (iter_count <= iter_max):\n",
    "            lat_aux = f_dlat(abs(D14_err/2))\n",
    "            lon_aux = f_dlon(P1[1] + copysign(lat_aux, D14_err), abs(D12_err/2))\n",
    "            P1      = [P1[0] + copysign(lon_aux, D12_err), P1[1] + copysign(lat_aux, D14_err)]\n",
    "            D12     = geodesic(lonlat(*P1), lonlat(*P2)).km\n",
    "            D14     = geodesic(lonlat(*P1), lonlat(*P4)).km\n",
    "            D12_err = D12 - D0\n",
    "            D14_err = D14 - D0\n",
    "            D_err   = max(abs(D12_err), abs(D14_err))\n",
    "            iter_count  += 1\n",
    "            if (iter_count == iter_max):\n",
    "                print (\" Warning at P1: iter_max exceeded. D_err = \", D_err) \n",
    "    \n",
    "        # Save to data_list:\n",
    "        data[x-1][y-1] = P1\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Fourth quadrant:\n",
    "for x in range(nx_1, NX-1, 1):\n",
    "    for y in range(ny_1, 0, -1):\n",
    "        # P1, P3 and P4 are defined:\n",
    "        P1 = data[x  ][y-1]\n",
    "        P3 = data[x+1][y  ]\n",
    "        P4 = data[x  ][y  ]\n",
    "        # P2 is approx. on the same longitude that P3, and approx. on the same latitude than P1:\n",
    "        P2      = [P3[0], P1[1]]\n",
    "        D12     = geodesic(lonlat(*P1), lonlat(*P2)).km\n",
    "        D23     = geodesic(lonlat(*P2), lonlat(*P3)).km\n",
    "        D12_err = D12 - D0\n",
    "        D23_err = D23 - D0\n",
    "        D_err   = max(abs(D12_err), abs(D23_err))\n",
    "        iter_count   = 0\n",
    "        while (abs(D_err) > abs_err) and (iter_count <= iter_max):\n",
    "            lat_aux = f_dlat(abs(D23_err/2))\n",
    "            lon_aux = f_dlon(P2[1] + copysign(lat_aux, D23_err), abs(D12_err/2))\n",
    "            P2      = [P2[0] - copysign(lon_aux, D12_err), P2[1] + copysign(lat_aux, D23_err)]\n",
    "            D12     = geodesic(lonlat(*P1), lonlat(*P2)).km\n",
    "            D23     = geodesic(lonlat(*P2), lonlat(*P3)).km\n",
    "            D12_err = D12 - D0\n",
    "            D23_err = D14 - D0\n",
    "            D_err   = max(abs(D12_err), abs(D23_err))\n",
    "            iter_count  += 1\n",
    "            if (iter_count == iter_max):\n",
    "                print (\" Warning at P2: iter_max exceeded. D_err = \", D_err) \n",
    "                print (x, y)\n",
    "    \n",
    "        # Save to data_list:\n",
    "        data[x+1][y-1] = P2\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Check distances.\n",
    "# Vertical:\n",
    "for x in range(0, NX, 1):\n",
    "    for y in range(0, NY-1, 1):\n",
    "        P1 = data[x][y  ]\n",
    "        P4 = data[x][y+1]\n",
    "        print(geodesic(lonlat(*P1), lonlat(*P4)))\n",
    "\n",
    "# Horizontal:\n",
    "for y in range(0, NY, 1):\n",
    "    for x in range(0, NX-1, 1):\n",
    "        P1 = data[x  ][y]\n",
    "        P2 = data[x+1][y]\n",
    "        print(geodesic(lonlat(*P1), lonlat(*P2)))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Save the grid.\n",
    "# Points:\n",
    "with open(FileNameOut1, 'w', newline = '') as fileCSVout:\n",
    "    writeCSV = csv.writer(fileCSVout, delimiter = ',')\n",
    "    for x in range(0, NX, 1):\n",
    "        writeCSV.writerows(data[x])\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create the dictionary for the geojson FeatureCollection.\n",
    "fc = {\"type\": \"FeatureCollection\", \"features\": []}\n",
    "\n",
    "# Add the polygons:\n",
    "i=0\n",
    "for x in range(0, NX-1, 1):\n",
    "    for y in range(0, NY-1, 1):\n",
    "        aux = [(p) for p in [data[x][y], data[x+1][y], data[x+1][y+1], data[x][y+1]]]\n",
    "        f ={\"type\": \"Feature\", \n",
    "            \"geometry\": {\"type\": \"Polygon\", \"coordinates\": [aux]},\n",
    "            \"properties\": {\"pol\": \"pol_from_grid\", \"id\": i}\n",
    "            }\n",
    "        fc[\"features\"].append(f)\n",
    "        i += 1\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Save as geojson.\n",
    "# json conversion is needed because the python dictionary uses '' while geojson requires \"\".\n",
    "with open(FileNameOut2, 'w') as output_file:\n",
    "    output_file.write(json.dumps(fc))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "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.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}