{
  "nbformat_minor": 0,
  "cells": [
    {
      "execution_count": null,
      "outputs": [],
      "cell_type": "code",
      "metadata": {
        "collapsed": false
      },
      "source": [
        "%matplotlib inline"
      ]
    },
    {
      "source": [
        "\nGeometric example\n------------------\n\nA small example script showing the usage of the 'geographic' coordinates type\nfor ordinary kriging on a sphere.\n\n"
      ],
      "cell_type": "markdown",
      "metadata": {}
    },
    {
      "execution_count": null,
      "outputs": [],
      "cell_type": "code",
      "metadata": {
        "collapsed": false
      },
      "source": [
        "from pykrige.ok import OrdinaryKriging\nimport numpy as np\n\n# Make this example reproducible:\nnp.random.seed(89239413)\n\n# Generate random data following a uniform spatial distribution\n# of nodes and a uniform distribution of values in the interval\n# [2.0, 5.5]:\nN = 7\nlon = 360.0*np.random.random(N)\nlat = 180.0/np.pi*np.arcsin(2*np.random.random(N)-1)\nz   = 3.5*np.random.rand(N) + 2.0\n\n# Generate a regular grid with 60\u00b0 longitude and 30\u00b0 latitude steps:\ngrid_lon = np.linspace(0.0, 360.0, 7)\ngrid_lat = np.linspace(-90.0, 90.0, 7)\n\n# Create ordinary kriging object:\nOK = OrdinaryKriging(lon, lat, z, variogram_model='linear', verbose=False,\n                     enable_plotting=False, coordinates_type='geographic')\n\n# Execute on grid:\nz1, ss1 = OK.execute('grid', grid_lon, grid_lat)\n\n# Create ordinary kriging object ignoring curvature:\nOK = OrdinaryKriging(lon, lat, z, variogram_model='linear', verbose=False,\n                     enable_plotting=False)\n\n# Execute on grid:\nz2, ss2 = OK.execute('grid', grid_lon, grid_lat)\n\n# Print data at equator (last longitude index will show periodicity):\nprint(\"Original data:\")\nprint(\"Longitude:\",lon.astype(int))\nprint(\"Latitude: \",lat.astype(int))\nprint(\"z:        \",np.array_str(z, precision=2))\nprint(\"\\nKrige at 60\u00b0 latitude:\\n======================\")\nprint(\"Longitude:\",grid_lon)\nprint(\"Value:    \",np.array_str(z1[5,:], precision=2))\nprint(\"Sigma\u00b2:   \",np.array_str(ss1[5,:], precision=2))\nprint(\"\\nIgnoring curvature:\\n=====================\")\nprint(\"Value:    \",np.array_str(z2[5,:], precision=2))\nprint(\"Sigma\u00b2:   \",np.array_str(ss2[5,:], precision=2))\n\n##====================================OUTPUT==================================\n# >>> Original data:\n# >>> Longitude: [122 166  92 138  86 122 136]\n# >>> Latitude:  [-46 -36 -25 -73 -25  50 -29]\n# >>> z:         [ 2.75  3.36  2.24  3.07  3.37  5.25  2.82]\n# >>> \n# >>> Krige at 60\u00b0 latitude:\n# >>> ======================\n# >>> Longitude: [   0.   60.  120.  180.  240.  300.  360.]\n# >>> Value:     [ 5.32  5.14  5.3   5.18  5.35  5.61  5.32]\n# >>> Sigma\u00b2:    [ 2.19  1.31  0.41  1.22  2.1   2.46  2.19]\n# >>> \n# >>> Ignoring curvature:\n# >>> =====================\n# >>> Value:     [ 4.55  4.72  5.25  4.82  4.61  4.53  4.48]\n# >>> Sigma\u00b2:    [ 3.77  1.99  0.39  1.84  3.52  5.43  7.5 ]\n#\n# We can see that the data point at longitude 122, latitude 50 correctly\n# dominates the kriged results, since it is the closest node in spherical\n# distance metric, as longitude differences scale with cos(latitude).\n# When kriging using longitude / latitude linearly, the value for grid points\n# with longitude values further away as longitude is now incorrectly\n# weighted equally as latitude."
      ]
    }
  ],
  "nbformat": 4,
  "metadata": {
    "language_info": {
      "codemirror_mode": {
        "version": 3,
        "name": "ipython"
      },
      "file_extension": ".py",
      "name": "python",
      "nbconvert_exporter": "python",
      "mimetype": "text/x-python",
      "version": "3.5.2",
      "pygments_lexer": "ipython3"
    },
    "kernelspec": {
      "language": "python",
      "name": "python3",
      "display_name": "Python 3"
    }
  }
}