From 640eda7e5d6449d6afc85ab3fe1a217a52de2458 Mon Sep 17 00:00:00 2001 From: Louis Moresi Date: Wed, 24 Apr 2024 13:33:59 +1000 Subject: [PATCH] Add notebook versions of Numpy --- .../NumpyAndScipy-IPYNB/0-NumpyAndScipy.ipynb | 39 + .../1-IntroductionToNumpy.ipynb | 805 ++++++++++++++++++ .../2-Application-TheGameOfLife.ipynb | 175 ++++ .../3-Discussion-TheGameOfLife.ipynb | 412 +++++++++ .../4-IntroductionToScipy.ipynb | 85 ++ .../5-ScipyInterpolate.ipynb | 324 +++++++ .../6-ScipySpatialAndMeshing.ipynb | 314 +++++++ .../NumpyAndScipy-IPYNB/7-ScipyOptimize.ipynb | 269 ++++++ .../8-Numpy+ScipyMiscellany.ipynb | 160 ++++ 9 files changed, 2583 insertions(+) create mode 100644 Notebooks/Themes/NumpyAndScipy-IPYNB/0-NumpyAndScipy.ipynb create mode 100644 Notebooks/Themes/NumpyAndScipy-IPYNB/1-IntroductionToNumpy.ipynb create mode 100644 Notebooks/Themes/NumpyAndScipy-IPYNB/2-Application-TheGameOfLife.ipynb create mode 100644 Notebooks/Themes/NumpyAndScipy-IPYNB/3-Discussion-TheGameOfLife.ipynb create mode 100644 Notebooks/Themes/NumpyAndScipy-IPYNB/4-IntroductionToScipy.ipynb create mode 100644 Notebooks/Themes/NumpyAndScipy-IPYNB/5-ScipyInterpolate.ipynb create mode 100644 Notebooks/Themes/NumpyAndScipy-IPYNB/6-ScipySpatialAndMeshing.ipynb create mode 100644 Notebooks/Themes/NumpyAndScipy-IPYNB/7-ScipyOptimize.ipynb create mode 100644 Notebooks/Themes/NumpyAndScipy-IPYNB/8-Numpy+ScipyMiscellany.ipynb diff --git a/Notebooks/Themes/NumpyAndScipy-IPYNB/0-NumpyAndScipy.ipynb b/Notebooks/Themes/NumpyAndScipy-IPYNB/0-NumpyAndScipy.ipynb new file mode 100644 index 0000000..5b190e0 --- /dev/null +++ b/Notebooks/Themes/NumpyAndScipy-IPYNB/0-NumpyAndScipy.ipynb @@ -0,0 +1,39 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "e8a4e5c8", + "metadata": {}, + "source": [ + "# Numpy and Scipy\n", + "\n", + "In this section we introduce the `numpy` and `scipy` packages that you will certainly use and will need to know for this course. \n", + "\n", + "`numpy` focuses on the efficient manipulation of structured data that in many languages would already be included as some sort of *array* data type. `numpy` provides routines to do many operations on arrays of data such as sorting, finding the largest or smallest element, transposing data and so on. It also provides many mathematical operations including the expected (element-by-element) functions (trigonometry, exponentials etc). \n", + "\n", + "`numpy` also provides an extensive overloading of the standard operators for manipulating lists of data (etc) that make array operations quite natural and compact in python. It is very important to understand what `numpy` is doing because it is almost impossible to avoid close encounters with these operators in everyday python. \n", + "\n", + "`numpy` has been very carefully optimised and, used correctly, will give you access to very efficient (mostly that means *fast*) code without having to leave the python environment. This means you should always look to see if there is a `numpy` function to do what you need *and never write a python loop if you can help it*. \n", + "\n", + " - [Introduction To Numpy](1-IntroductionToNumpy.md) - This notebook goes into more detail about why we need `numpy`, what the operator overloading looks like, and explores the speed benefits of using `numpy`. The other benefits are readability of your code (as long as you understand the operator overloading) and reliability because you are using pre-built code, not writing your own loops. \n", + " - [The game of life 1](2-Application-TheGameOfLife.md) - This notebook sets up an exercise: the game of life, to show how `numpy` works in a real problem. We will see how to approach problems of structured data so that we can most of the work efficiently with `numpy`.\n", + " - [The game of life 2](3-Discussion-TheGameOfLife.md) - Here we have a chance to reflect on the choices we made in the previous notebook and think about how things could be done better (or differently, anyway). \n", + "\n", + " These notebooks are not intended to be comprehensive because you will have lots of opportunity to practice looking at `numpy` code when we look at satellite images and basemaps in the [`cartopy` mapping](Notebooks/Themes/Mapping/0-Maps_with_Cartopy.md) examples, and the [`stripy` spherical meshing](Notebooks/Themes/SphericalMeshing/0-Stripy.md) examples. \n", + "\n", + " `scipy` is a collection of useful algorithms loosely bundled together. If you have ever come across the book \"Numerical Recipes\" then you will know what `scipy` is trying to do: there is usually some reasonable implementation in `scipy` for whatever problem you are faced with: interpolation, optimisation, meshing, image manipulation. It is neither comprehensive nor definitive but, like a Swiss-Army Knife, more useful in a pinch than most other tools.\n", + "\n", + " - [Introduction](4-IntroductionToScipy.md) - Really just a list of modules that `scipy` provides. Each of the modules is a collection of related functionality that has a mostly-consistent user interface." + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/Notebooks/Themes/NumpyAndScipy-IPYNB/1-IntroductionToNumpy.ipynb b/Notebooks/Themes/NumpyAndScipy-IPYNB/1-IntroductionToNumpy.ipynb new file mode 100644 index 0000000..6a34c2c --- /dev/null +++ b/Notebooks/Themes/NumpyAndScipy-IPYNB/1-IntroductionToNumpy.ipynb @@ -0,0 +1,805 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "613ddf55", + "metadata": {}, + "source": [ + "# Numpy data structures\n", + "\n", + "When we looked at python data structures, it was obvious that the only way to deal with arrays of values (matrices / vectors etc) would be via lists and lists of lists.\n", + "\n", + "This is slow and inefficient in both execution and writing code. \n", + "\n", + "`numpy` attempts to fix this." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "d1ac3f36", + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "\n", + "## This is a list of everything in the module\n", + "np.__all__" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "81a40ffe", + "metadata": {}, + "outputs": [], + "source": [ + "an_array = np.array([0,1,2,3,4,5,6], dtype=np.float)\n", + "\n", + "print (an_array)\n", + "print( type(an_array))\n", + "help(an_array)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "bd3d9103", + "metadata": {}, + "outputs": [], + "source": [ + "print(an_array.dtype)\n", + "print(an_array.shape)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "bc33c25e", + "metadata": {}, + "outputs": [], + "source": [ + "A = np.zeros((4,4))\n", + "print (A)\n", + "\n", + "print (A.shape)\n", + "print (A.diagonal())\n", + "\n", + "A[0,0] = 2.0\n", + "print (A)\n", + "\n", + "np.fill_diagonal(A, 1.0)\n", + "print (A)\n", + "\n", + "B = A.diagonal()\n", + "B[0] = 2.0" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "bc423f17", + "metadata": {}, + "outputs": [], + "source": [ + "for i in range(0,A.shape[0]):\n", + " A[i,i] = 1.0\n", + "print (A)\n", + "print()\n", + "\n", + "A[:,2] = 2.0\n", + "print (A)\n", + "print()\n", + "\n", + "A[2,:] = 4.0\n", + "print (A)\n", + "print()\n", + "\n", + "print (A.T)\n", + "print()\n", + "\n", + "A[...] = 0.0\n", + "print (A)\n", + "print()\n", + "\n", + "for i in range(0,A.shape[0]):\n", + " A[i,:] = float(i)\n", + " \n", + "print (A)\n", + "print()\n", + "\n", + "for i in range(0,A.shape[0]):\n", + " A[i,:] = i\n", + " \n", + "print (A)\n", + "print()\n", + "\n", + "print (A[::2,::2] )\n", + "\n", + "\n", + "print (A[::-1,::-1])" + ] + }, + { + "cell_type": "markdown", + "id": "a1cbfbba", + "metadata": {}, + "source": [ + "## Speed" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "2b993bec", + "metadata": {}, + "outputs": [], + "source": [ + "%%timeit\n", + "\n", + "B = np.zeros((1000,1000))\n", + "for i in range(0,1000):\n", + " for j in range(0,1000):\n", + " B[i,j] = 2.0\n", + " \n", + " " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "eaedbd19", + "metadata": {}, + "outputs": [], + "source": [ + "%%timeit\n", + "\n", + "B = np.zeros((1000,1000))\n", + "B[:,:] = 2.0" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "c484e3d7", + "metadata": {}, + "outputs": [], + "source": [ + "%%timeit\n", + "\n", + "B = np.zeros((1000,1000))\n", + "B[...] = 2.0" + ] + }, + { + "cell_type": "markdown", + "id": "914e423f", + "metadata": {}, + "source": [ + "## Views of the data (are free)\n", + "\n", + "It costs very little to look at data in a different way (e.g. to view a 2D array as a 1D vector).\n", + "\n", + "Making a copy is a different story" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "89c133e1", + "metadata": {}, + "outputs": [], + "source": [ + "print( A.reshape((2,8)))\n", + "print()\n", + "\n", + "print( A.reshape((-1)))\n", + "print( A.ravel())\n", + "print()\n", + "\n", + "print( A.reshape((1,-1)))\n", + "print()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "9915b009", + "metadata": {}, + "outputs": [], + "source": [ + "%%timeit \n", + "\n", + "A.reshape((1,-1))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "f668ac45", + "metadata": {}, + "outputs": [], + "source": [ + "%%timeit\n", + "\n", + "elements = A.shape[0]*A.shape[1]\n", + "B = np.empty(elements)\n", + "B[...] = A[:,:].reshape(elements)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "7d9d0134", + "metadata": {}, + "outputs": [], + "source": [ + "%%timeit\n", + "\n", + "elements = A.shape[0]*A.shape[1]\n", + "B = np.empty(elements)\n", + "\n", + "for i in range(0,A.shape[0]):\n", + " for j in range(0,A.shape[1]):\n", + " B[i+j*A.shape[1]] = A[i,j]" + ] + }, + { + "cell_type": "markdown", + "id": "d4fe6b1a", + "metadata": {}, + "source": [ + "__Exercise:__ Try this again for a 10000x10000 array" + ] + }, + { + "cell_type": "markdown", + "id": "6a76b384", + "metadata": {}, + "source": [ + "## Indexing / broadcasting\n", + "\n", + "In numpy, we can index an array by explicitly specifying elements, by specifying slices, by supplying lists of indices (or arrays), we can also supply a boolean array of the same shape as the original array which will select / return an array of all those entries where `True` applies.\n", + "\n", + "Although some of these might seem difficult to use, they are often the result of other numpy operations. For example `np.where` converts a truth array to a list of indices." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "9b7e5394", + "metadata": {}, + "outputs": [], + "source": [ + "AA = np.zeros((100,100))\n", + "\n", + "AA[10,11] = 1.0\n", + "AA[99,1] = 2.0\n", + "\n", + "cond = np.where(AA >= 1.0)\n", + "print (cond)\n", + "print (AA[cond])\n", + "print (AA[ AA >= 1])\n", + "print(AA>=1.0)" + ] + }, + { + "cell_type": "markdown", + "id": "e2f52c77", + "metadata": {}, + "source": [ + "---\n", + "\n", + "Broadcasting is a way of looping on arrays which have \"compatible\" but unequal sizes.\n", + "\n", + "For example, the element-wise multiplication of 2 arrays\n", + "\n", + "```python\n", + "a = np.array([1.0, 2.0, 3.0])\n", + "b = np.array([2.0, 2.0, 2.0]) \n", + "print a * b\n", + "```\n", + "has an equivalent:\n", + "\n", + "```python\n", + "a = np.array([1.0, 2.0, 3.0])\n", + "b = np.array([2.0]) \n", + "print a * b\n", + "```\n", + "\n", + "or \n", + "\n", + "```python\n", + "a = np.array([1.0, 2.0, 3.0])\n", + "b = 2.0 \n", + "print a * b\n", + "```\n", + "\n", + "in which the \"appropriate\" interpretation of `b` is made in each case to achieve the result." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "5035c2f0", + "metadata": {}, + "outputs": [], + "source": [ + "a = np.array([1.0, 2.0, 3.0])\n", + "b = np.array([2.0, 2.0, 2.0]) \n", + "\n", + "print( a * b)\n", + "\n", + "b = np.array([2.0]) \n", + "\n", + "print (a * b)\n", + "print (a * 2.0)" + ] + }, + { + "cell_type": "markdown", + "id": "5dc71cc4", + "metadata": {}, + "source": [ + "Arrays are compatible as long as each of their dimensions (shape) is either equal to the other or 1. \n", + "\n", + "Thus, above, the multplication works when `a.shape` is `(1,3)` and `b.shape` is either `(1,3)` or `(1,1)`\n", + "\n", + "(Actually, these are `(3,)` and `(1,)` in the examples above ..." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "bfbe6143", + "metadata": {}, + "outputs": [], + "source": [ + "print (a.shape)\n", + "print (b.shape)\n", + "print ((a*b).shape)\n", + "print ((a+b).shape)\n", + "\n", + "\n", + "aa = a.reshape(1,3)\n", + "bb = b.reshape(1,1)\n", + "\n", + "print (aa.shape)\n", + "print (bb.shape)\n", + "print ((aa*bb).shape)\n", + "print ((aa+bb).shape)" + ] + }, + { + "cell_type": "markdown", + "id": "6b1e625e", + "metadata": {}, + "source": [ + "In multiple dimensions, the rule applies but, perhaps, is less immediately intuitive:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "af174f2c", + "metadata": {}, + "outputs": [], + "source": [ + "a = np.array([[ 0.0, 0.0, 0.0],\n", + " [10.0,10.0,10.0],\n", + " [20.0,20.0,20.0],\n", + " [30.0,30.0,30.0]])\n", + "b = np.array([[1.0,2.0,3.0]])\n", + "print (a + b)\n", + "print ()\n", + "print (a.shape)\n", + "print (b.shape)\n", + "print ((a+b).shape)" + ] + }, + { + "cell_type": "markdown", + "id": "94d2f521", + "metadata": {}, + "source": [ + "Note that this also works for" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "251bc0da", + "metadata": {}, + "outputs": [], + "source": [ + "a = np.array([[ 0.0, 0.0, 0.0],\n", + " [10.0,10.0,10.0],\n", + " [20.0,20.0,20.0],\n", + " [30.0,30.0,30.0]])\n", + "\n", + "b = np.array([1.0,2.0,3.0])\n", + "print (a + b)\n", + "print ()\n", + "print (a.shape)\n", + "print (b.shape)\n", + "print ((a+b).shape)" + ] + }, + { + "cell_type": "markdown", + "id": "daf16cb4", + "metadata": {}, + "source": [ + "But not for" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "068b8736", + "metadata": {}, + "outputs": [], + "source": [ + "a = np.array([[ 0.0, 0.0, 0.0],\n", + " [10.0,10.0,10.0],\n", + " [20.0,20.0,20.0],\n", + " [30.0,30.0,30.0]])\n", + "\n", + "b = np.array([[1.0],[2.0],[3.0]])\n", + "\n", + "print (a.shape)\n", + "print (b.shape)\n", + "print ((a+b).shape)" + ] + }, + { + "cell_type": "markdown", + "id": "98ccf085", + "metadata": {}, + "source": [ + "## Vector Operations" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "e8dab0e2", + "metadata": {}, + "outputs": [], + "source": [ + "X = np.arange(0.0, 2.0*np.pi, 0.0001)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "63949ceb", + "metadata": {}, + "outputs": [], + "source": [ + "print (X)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ed3316df", + "metadata": {}, + "outputs": [], + "source": [ + "import math\n", + "\n", + "math.sin(X)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "bc799bfe", + "metadata": {}, + "outputs": [], + "source": [ + "np.sin(X)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "8d2dddb0", + "metadata": {}, + "outputs": [], + "source": [ + "S = np.sin(X)\n", + "C = np.cos(X)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "138b667a", + "metadata": {}, + "outputs": [], + "source": [ + "S2 = S**2 + C**2\n", + "print (S2)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "9e501edb", + "metadata": {}, + "outputs": [], + "source": [ + "print (S2 - 1.0)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "7dfa5200", + "metadata": {}, + "outputs": [], + "source": [ + "test = np.isclose(S2,1.0)\n", + "print (test)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "2afaca1a", + "metadata": {}, + "outputs": [], + "source": [ + "print (np.where(test == False))\n", + "print (np.where(S2 == 0.0))" + ] + }, + { + "cell_type": "markdown", + "id": "44e76389", + "metadata": {}, + "source": [ + "---\n", + "\n", + "__Exercise__: find out how long it takes to compute the sin, sqrt, power of a 1000000 length vector (array). How does this speed compare to using the normal `math` functions element by element in the array ? What happens if X is actually a complex array ?\n", + "\n", + "__Hints:__ you might find it useful to know about:\n", + " - `np.linspace` v `np.arange`\n", + " - `np.empty_like` or `np.zeros_like`\n", + " - the python `enumerate` function\n", + " - how to write a table in markdown\n", + " \n", + "| description | time | notes |\n", + "|-----------------|--------|-------|\n", + "| `np.sin` | ? | |\n", + "| `math.sin` | ? | |\n", + "| | ? | - |" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "0098c885", + "metadata": { + "solution": "hidden", + "solution_first": true + }, + "outputs": [], + "source": [ + "X = np.linspace(0.0, 2.0*np.pi, 10000000)\n", + "print (X.shape)\n", + "\n", + "# ... " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "d7cdbf0a", + "metadata": { + "solution": "hidden" + }, + "outputs": [], + "source": [ + "%%timeit\n", + "S = np.sin(X)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "cc863008", + "metadata": { + "solution": "hidden" + }, + "outputs": [], + "source": [ + "%%timeit\n", + "\n", + "S = np.empty_like(X)\n", + "for i, x in enumerate(X):\n", + " S[i] = math.sin(x)\n", + " " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "fec47c66", + "metadata": {}, + "outputs": [], + "source": [ + "XX = np.ones((100,100,100))\n", + "XX.shape" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b7240192", + "metadata": {}, + "outputs": [], + "source": [ + "%%timeit\n", + "np.sin(XX)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "01a5565b", + "metadata": {}, + "outputs": [], + "source": [ + "%%timeit\n", + "\n", + "for i in range(0,100):\n", + " for j in range(0,100):\n", + " for k in range(0,100):\n", + " math.sin(XX[i,j,k])\n", + " " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4a3dca53", + "metadata": { + "solution": "hidden" + }, + "outputs": [], + "source": [ + "X = np.linspace(0.0, 2.0*np.pi, 10000000)\n", + "Xj = X + 1.0j\n", + "print (Xj.shape, Xj.dtype)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "d29c851a", + "metadata": { + "solution": "hidden" + }, + "outputs": [], + "source": [ + "%%timeit\n", + "\n", + "Sj = np.sin(Xj)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "0c3b3dcf", + "metadata": { + "solution": "hidden" + }, + "outputs": [], + "source": [ + "import cmath" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "5d5aea18", + "metadata": { + "solution": "hidden" + }, + "outputs": [], + "source": [ + "%%timeit\n", + "\n", + "Sj = np.empty_like(Xj)\n", + "for i, x in enumerate(Xj):\n", + " Sj[i] = cmath.sin(x)\n", + " " + ] + }, + { + "cell_type": "markdown", + "id": "9afcc350", + "metadata": { + "solution": "hidden", + "solution_first": true + }, + "source": [ + "__Exercise__: look through the functions below from numpy, choose 3 of them and work out how to use them on arrays of data. Write a few lines to explain what you find. \n", + "\n", + " - `np.max` v. `np.argmax`\n", + " - `np.where`\n", + " - `np.logical_and`\n", + " - `np.fill_diagonal`\n", + " - `np.count_nonzero`\n", + " - `np.isinf` and `np.isnan`\n", + " \n", + " \n", + "Here is an example: \n", + "\n", + "`np.concatenate` takes a number of arrays and glues them together. For 1D arrays this is simple:\n", + "\n", + "```python\n", + "\n", + "A = np.array([1.0,1.0,1.0,1.0])\n", + "B = np.array([2.0,2.0,2.0,2.0])\n", + "C = np.array([3.0,3.0,3.0,3.0])\n", + "\n", + "R = np.concatenate((A,B,C))\n", + "\n", + "# array([ 1., 1., 1., 1., 2., 2., 2., 2., 3., 3., 3., 3.])\n", + "```\n", + "\n", + "an equivalent statement is `np.hstack((A,B,C))` but note the difference with `np.vstack((A,B,C))`\n", + "\n", + "With higher dimensional arrays, the gluing takes place along one `axis`:\n", + "\n", + "```python\n", + "A = np.array(([1.0,1.0,1.0,1.0],[2.0,2.0,2.0,2.0]))\n", + "B = np.array(([3.0,3.0,3.0,3.0],[4.0,4.0,4.0,4.0]))\n", + "C = np.array(([5.0,5.0,5.0,5.0],[6.0,6.0,6.0,6.0]))\n", + "\n", + "R = np.concatenate((A,B,C))\n", + "print R\n", + "print\n", + "R = np.concatenate((A,B,C), axis=1)\n", + "print R\n", + "\n", + "```" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "cf313947", + "metadata": { + "solution": "hidden" + }, + "outputs": [], + "source": [ + "# Test the results here \n", + "\n", + "A = np.array(([1.0,1.0,1.0,1.0],[2.0,2.0,2.0,2.0]))\n", + "B = np.array(([3.0,3.0,3.0,3.0],[4.0,4.0,4.0,4.0]))\n", + "C = np.array(([5.0,5.0,5.0,5.0],[6.0,6.0,6.0,6.0]))\n", + "\n", + "R = np.concatenate((A,B,C))\n", + "print R\n", + "print\n", + "R = np.concatenate((A,B,C), axis=1)\n", + "print R" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/Notebooks/Themes/NumpyAndScipy-IPYNB/2-Application-TheGameOfLife.ipynb b/Notebooks/Themes/NumpyAndScipy-IPYNB/2-Application-TheGameOfLife.ipynb new file mode 100644 index 0000000..dd3a35f --- /dev/null +++ b/Notebooks/Themes/NumpyAndScipy-IPYNB/2-Application-TheGameOfLife.ipynb @@ -0,0 +1,175 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "08889dfa", + "metadata": {}, + "source": [ + "# The game of life \n", + "\n", + "The universe of the Game of Life ([John Conway](https://en.wikipedia.org/wiki/John_Horton_Conway)) is an infinite two-dimensional orthogonal grid of square cells, each of which is in one of two possible states, live or dead. Every cell interacts with its eight neighbours, which are the cells that are directly horizontally, vertically, or diagonally adjacent. At each step in time, the following transitions occur:\n", + "\n", + "* Any live cell with fewer than two live neighbours dies, as if by needs caused by underpopulation.\n", + "* Any live cell with more than three live neighbours dies, as if by overcrowding.\n", + "* Any live cell with two or three live neighbours lives, unchanged, to the next generation.\n", + "* Any dead cell with exactly three live neighbours becomes a live cell.\n", + "\n", + "See [this page](http://web.stanford.edu/~cdebs/GameOfLife/) for some more information. And, note, thanks to Dan Sandiford for setting up this example." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "574c208e", + "metadata": {}, + "outputs": [], + "source": [ + "%pylab inline\n", + "import numpy as np\n", + "import matplotlib.pyplot as plt" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "8e71af57", + "metadata": {}, + "outputs": [], + "source": [ + "start = np.array([[1,0,0,0,0,0],\n", + " [0,0,0,1,0,0],\n", + " [0,1,0,1,0,0],\n", + " [0,0,1,1,0,0],\n", + " [0,0,0,0,0,0],\n", + " [0,0,0,0,0,1]])" + ] + }, + { + "cell_type": "markdown", + "id": "06b684c9", + "metadata": {}, + "source": [ + "We will talk more about plotting later, but for now we can use this without digging deeper:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "715cd367", + "metadata": {}, + "outputs": [], + "source": [ + "plt.imshow(start, interpolation='nearest', cmap=\"gray\") " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "6ad00cbb", + "metadata": {}, + "outputs": [], + "source": [ + "print start[4:8,4:8] # neighbours of start[5,5]\n", + "print start[1:4,1:4] # neighbours of start[2,2]\n", + "#print start[?:?] # neighbours of start[1,1]\n", + "#print start[?:?] # neighbours of start[0,0]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "19c24617", + "metadata": {}, + "outputs": [], + "source": [ + "live_neighbours = np.empty(start.shape)\n", + "for index, value in np.ndenumerate(start):\n", + " #Need to add 2, becase the slicing works like 'up to but not including'\n", + " x0 = max(0,(index[0]-1))\n", + " x1 = max(0,(index[0]+2))\n", + " y0 = max(0,(index[1]-1))\n", + " y1 = max(0,(index[1]+2))\n", + " subarray = start[x0:x1, y0:y1]\n", + " live_neighbours[index] = subarray.sum() - value # need to subtract actual value at that cell..." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "d59a7f17", + "metadata": {}, + "outputs": [], + "source": [ + "live_neighbours" + ] + }, + { + "cell_type": "markdown", + "id": "9802b3b7", + "metadata": {}, + "source": [ + "---\n", + "\n", + "__Exercise:__ Your task is to write a function that \"runs\" the game of life. This should be possible by filling out the two function templates below. \n", + "\n", + " - conditions: boundaries are always dead" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "d2f888e9", + "metadata": {}, + "outputs": [], + "source": [ + "def get_neighbours(start):\n", + "\"\"\"\n", + "This function gets the number of live neighbours in the binary array start\n", + "start : np.ndarray\n", + "\"\"\" " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "dc6b660f", + "metadata": {}, + "outputs": [], + "source": [ + "def game_of_life(start, n):\n", + " \"\"\"\n", + " this function runs the game of life for n steps...\n", + " start : np.ndarray (0s and 1s)\n", + " n: int number of steps \n", + " \"\"\"\n", + " assert (1>= start.min() >= 0) and (1>= start.max() >= 0), \"array must be ones and zeros\"\n", + " \n", + " current = np.copy(start)\n", + " \n", + " while n:\n", + " neighbours = get_neighbours(current)\n", + " \n", + " for index, value in np.ndenumerate(current):\n", + " print(index, value)\n", + " # Apply the rules to current\n", + " if ...\n", + " \n", + " else ...\n", + " \n", + " n -= 1 \n", + " \n", + " \n", + " return current" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python3", + "language": "python", + "name": "python3" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/Notebooks/Themes/NumpyAndScipy-IPYNB/3-Discussion-TheGameOfLife.ipynb b/Notebooks/Themes/NumpyAndScipy-IPYNB/3-Discussion-TheGameOfLife.ipynb new file mode 100644 index 0000000..d463ecf --- /dev/null +++ b/Notebooks/Themes/NumpyAndScipy-IPYNB/3-Discussion-TheGameOfLife.ipynb @@ -0,0 +1,412 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": null, + "id": "59ba1a16", + "metadata": { + "hide_input": true + }, + "outputs": [], + "source": [ + "%matplotlib inline\n", + "import numpy as np\n", + "import matplotlib.pyplot as plt" + ] + }, + { + "cell_type": "markdown", + "id": "2a50f790", + "metadata": {}, + "source": [ + "# The game of life (Advanced)\n", + "\n", + "This is an example of a __cellular automaton__ - a collection of rules that act on the state of a system and determine the updated state. These often have interesting, _emergent_ behaviours that are not obviously inherent in the rules.\n", + "\n", + "There is a similarity with Molecular Dynamics, Discrete Element Modeling, Agent-Based modeling and so on where simple rules acting locally produce bulk properties of the system. \n", + "\n", + "The game of life is simple enough in rules but can produce interesting emergent behaviour. It is a good template to explore some of the issues we always face in building models.\n", + "\n", + "Let's look a building a more challenging version of the problem - larger, faster, better analysis." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a951c6c5", + "metadata": {}, + "outputs": [], + "source": [ + "# The rules of the game\n", + "\n", + "size = 100\n", + "lonely = 2\n", + "crowded = 4\n", + "breeding = 3" + ] + }, + { + "cell_type": "markdown", + "id": "95ed937b", + "metadata": {}, + "source": [ + "We don't want to keep setting this thing up by hand. Here is a way to make an array of zeros and ones at random" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "e9c769fd", + "metadata": {}, + "outputs": [], + "source": [ + "state = np.round(np.random.random_sample((size,size) )).astype(int)\n", + "\n", + "print state\n", + "print \n", + "print np.count_nonzero(state)\n", + "print" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "9691c7d0", + "metadata": {}, + "outputs": [], + "source": [ + "plt.imshow(state, interpolation='nearest', cmap=\"gray\") " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "54fb765f", + "metadata": {}, + "outputs": [], + "source": [ + "def live_neighbours(state):\n", + " ln = np.empty_like(state)\n", + " for index, value in np.ndenumerate(state):\n", + " #Need to add 2, becase the slicing works like 'up to but not including'\n", + " x0 = max(0,(index[0]-1))\n", + " x1 = max(0,(index[0]+2))\n", + " y0 = max(0,(index[1]-1))\n", + " y1 = max(0,(index[1]+2))\n", + " subarray = state[x0:x1, y0:y1]\n", + " ln[index] = subarray.sum() - value # need to subtract actual value at that cell...\n", + " \n", + " return ln" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4c70f8a0", + "metadata": {}, + "outputs": [], + "source": [ + "ln = live_neighbours(state)\n", + "\n", + "print ln\n", + "\n", + "plt.imshow(ln, interpolation='nearest', cmap=\"jet\") \n", + "plt.colorbar()" + ] + }, + { + "cell_type": "markdown", + "id": "50a2de0b", + "metadata": {}, + "source": [ + "Let's use numpy to locate *quickly* the regions where the death conditions apply \n", + "\n", + "```python\n", + "death1 = np.where(ln >= crowded)\n", + "```\n", + "\n", + "Note the form that this returns - the coordinates where the condition applies" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "48d284f1", + "metadata": {}, + "outputs": [], + "source": [ + "death1 = np.where(ln >= crowded)\n", + "death2 = np.where(ln < lonely)\n", + "\n", + "print death1" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "19cdc891", + "metadata": {}, + "outputs": [], + "source": [ + "dying_cells = np.zeros_like(state)\n", + "dying_cells[death1] = 1\n", + "dying_cells[death2] = 1\n", + "\n", + "## \n", + "\n", + "plt.imshow(dying_cells, interpolation='nearest', cmap=\"gray_r\") \n", + "plt.colorbar()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "9b91b933", + "metadata": {}, + "outputs": [], + "source": [ + "birth = np.where(ln == breeding)\n", + "\n", + "baby_cells = np.zeros_like(state)\n", + "baby_cells[birth] = 1\n", + "\n", + "plt.imshow(baby_cells, interpolation='nearest', cmap=\"gray_r\") \n", + "plt.colorbar()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "6f3e99b5", + "metadata": {}, + "outputs": [], + "source": [ + "## These should be complementary\n", + "\n", + "check = dying_cells + baby_cells\n", + "print check.max()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "e823e632", + "metadata": { + "solution": "hidden" + }, + "outputs": [], + "source": [ + "# Logical operations on the entire array to update state\n", + "\n", + "new_state = state.copy()\n", + "new_state[death1] = 0\n", + "new_state[death2] = 0\n", + "new_state[birth] = 1\n", + "\n", + "plt.imshow(new_state, interpolation='nearest', cmap=\"gray_r\",) \n", + "plt.colorbar()\n", + "plt.savefig(\"Frame-1.png\", dpi=300, format=\"png\" )" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "64d482b4", + "metadata": {}, + "outputs": [], + "source": [ + "## Bundle this up to make a sequence of images\n", + "\n", + "state_list = [state.copy()]\n", + "\n", + "for i in range(0, 250):\n", + " \n", + " if i%50 == 0:\n", + " print i\n", + "\n", + " ln = live_neighbours(state)\n", + "\n", + " state[ln >= crowded] = 0\n", + " state[ln < lonely] = 0\n", + " state[ln == breeding] = 1\n", + "\n", + " state_list.append(state.copy())\n", + " " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "80a9f737", + "metadata": {}, + "outputs": [], + "source": [ + "state_list" + ] + }, + { + "cell_type": "markdown", + "id": "5efce189", + "metadata": { + "solution": "hidden", + "solution_first": true + }, + "source": [ + "### Exercise\n", + "\n", + "\n", + "__Class Exercise:__ Now let's try to make a series of frames as images for a number of frames of the game of life. We could also try changing the parameters of the game to see what patterns develop. \n", + "\n", + " - What resolution can we reach ? \n", + " - Can we work out ways to make this run faster ?\n", + " - What if there is more than one happy state where new cells grow ?" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "eb276e84", + "metadata": { + "solution": "hidden" + }, + "outputs": [], + "source": [ + "def update_state(state):\n", + " new_state = state.copy()\n", + "\n", + " ln = live_neighbours(state)\n", + " \n", + " death1 = np.where(ln >= crowded)\n", + " death2 = np.where(ln < lonely)\n", + " birth = np.where(ln == breeding)\n", + "\n", + " new_state[death1] = 0\n", + " new_state[death2] = 0\n", + " new_state[birth] = 1\n", + " \n", + " return new_state\n", + " " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "d6ac76af", + "metadata": { + "solution": "hidden" + }, + "outputs": [], + "source": [ + "# Initialise\n", + "\n", + "state = np.round(np.random.random_sample((size,size) )).astype(int)\n", + "steps = 100\n", + "\n", + "# Loops and save image\n", + "\n", + "for step in range(0, steps):\n", + " \n", + " new_state = update_state(state)\n", + " print \"Step {} - living cells - {}\".format(step, np.count_nonzero(new_state))\n", + " \n", + " plt.imshow(new_state, interpolation='nearest', cmap=\"gray_r\",) \n", + " plt.savefig(\"Frame-{:04d}.png\".format(step), dpi=300, format=\"png\" )\n", + " \n", + " # We could keep all of these for later or discard\n", + " state = new_state" + ] + }, + { + "cell_type": "markdown", + "id": "1650bf55", + "metadata": {}, + "source": [ + "### Speedup\n", + "\n", + "What is taking the time in this computation ?\n", + "\n", + " - Boundaries always reset to zero, do away with submatrix tests ?\n", + " - Fastest way to reset boundaries ?\n", + " - Logic tests or multiply or add / cap the max" + ] + }, + { + "cell_type": "markdown", + "id": "97656061", + "metadata": {}, + "source": [ + "## Animation" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "3afff391", + "metadata": {}, + "outputs": [], + "source": [ + "from matplotlib import animation, rc\n", + "from IPython.display import HTML\n", + "\n", + "fig, ax = plt.subplots()\n", + "fig.set_size_inches((7,7))\n", + "animage = ax.imshow(state_list[0], interpolation='bilinear', cmap=plt.cm.bone)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a5425cb7", + "metadata": {}, + "outputs": [], + "source": [ + "def init():\n", + " animage.set_data(state_list[0])\n", + " return (animage,)\n", + "\n", + "def animate(i):\n", + " animage.set_data(state_list[i])\n", + " return (animage,)\n", + "\n", + "anim = animation.FuncAnimation(fig, animate, init_func=init,\n", + " frames=len(state_list), interval=50, \n", + " blit=True)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "d1406994", + "metadata": {}, + "outputs": [], + "source": [ + "HTML(anim.to_jshtml())" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "5d23f255", + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "63210200", + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/Notebooks/Themes/NumpyAndScipy-IPYNB/4-IntroductionToScipy.ipynb b/Notebooks/Themes/NumpyAndScipy-IPYNB/4-IntroductionToScipy.ipynb new file mode 100644 index 0000000..002dd16 --- /dev/null +++ b/Notebooks/Themes/NumpyAndScipy-IPYNB/4-IntroductionToScipy.ipynb @@ -0,0 +1,85 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "fa992992", + "metadata": {}, + "source": [ + "# scipy \n", + "\n", + "\n", + "`scipy` is a collection of tools that are commonly encountered in scientific environments. The implementations are generally robust and reliable. The documentation varies in quality but is generally helpful. \n", + "\n", + "It is impossible to cover all the aspects of the entire package because it includes modules to cover an enormous range of functionality:\n", + "\n", + " - `cluster`: clustering algorithms (data analysis)\n", + " - `constants`: physical / mathematical constants\n", + " - `fftpack`: fourier transforms\n", + " - `integrate`: numerical integration routines\n", + " - `interpolate`: various data interpolation routines\n", + " - `linalg`: linear algebra routines (matrices etc)\n", + " - `ndimage`: N dimensional image manipulation and resizing (also array manipulations)\n", + " - `odr`: Orthogonal distance regression\n", + " - `optimize`: Optimization / root finding\n", + " - `signal`: Signal processing\n", + " - `sparse`: Sparse matrix packages\n", + " - `spatial`: Spatial data / triangulation / Voronoi diagrams \n", + " - `special`: Special (mathematical) functions \n", + " - `stats`: Statistics routines\n", + " \n", + " \n", + "Here we will just skim some of the available sub-modules and go through some of the examples provided." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4a80fb6f", + "metadata": {}, + "outputs": [], + "source": [ + "## importing scipy does very little (mostly the numpy namespace)\n", + "import numpy as np\n", + "import scipy\n", + "\n", + "## to use the modules import each one separately\n", + "from scipy import constants as constants\n", + "\n", + "print (constants.electron_mass)\n", + "print ()\n", + "print (constants.find(\"electron mass\"))\n", + "print ()\n", + "print (constants.value(\"electron mass\"),\" \", constants.unit(\"electron mass\"))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "f604a8e1", + "metadata": {}, + "outputs": [], + "source": [ + "## Note, this is actually the numpy array as well\n", + "\n", + "help(scipy.ndarray)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "24ae1e23", + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/Notebooks/Themes/NumpyAndScipy-IPYNB/5-ScipyInterpolate.ipynb b/Notebooks/Themes/NumpyAndScipy-IPYNB/5-ScipyInterpolate.ipynb new file mode 100644 index 0000000..d0e29e6 --- /dev/null +++ b/Notebooks/Themes/NumpyAndScipy-IPYNB/5-ScipyInterpolate.ipynb @@ -0,0 +1,324 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "77bc70c8", + "metadata": {}, + "source": [ + "# `scipy.interpolate`\n", + "\n", + "This module provides general interpolation capability for data in 1, 2, and higher dimensions. This list of features is from the documentation:\n", + "\n", + " - A class representing an interpolant (interp1d) in 1-D, offering several interpolation methods.\n", + "\n", + " - Convenience function griddata offering a simple interface to interpolation in N dimensions (N = 1, 2, 3, 4, ...). Object-oriented interface for the underlying routines is also available.\n", + " \n", + " - Functions for 1- and 2-dimensional (smoothed) cubic-spline interpolation, based on the FORTRAN library FITPACK. There are both procedural and object-oriented interfaces for the FITPACK library.\n", + " \n", + " - Interpolation using Radial Basis Functions." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "0064684b", + "metadata": {}, + "outputs": [], + "source": [ + "%matplotlib inline\n", + "import matplotlib.pyplot as plt\n", + "import numpy as np" + ] + }, + { + "cell_type": "markdown", + "id": "2c92627c", + "metadata": {}, + "source": [ + "## 1D data" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "33ba8b7b", + "metadata": {}, + "outputs": [], + "source": [ + "from scipy.interpolate import interp1d" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "7022a712", + "metadata": {}, + "outputs": [], + "source": [ + "x = np.linspace(0, 10, num=11, endpoint=True)\n", + "y = np.cos(-x**2/9.0)\n", + "f = interp1d(x, y, kind='linear') # default if kind=None\n", + "f2 = interp1d(x, y, kind='cubic')\n", + "f3 = interp1d(x, y, kind='nearest')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "31289fdc", + "metadata": {}, + "outputs": [], + "source": [ + "xnew = np.linspace(0, 10, num=41, endpoint=True)\n", + "plt.plot(x, y, 'o', xnew, f(xnew), '-', xnew, f2(xnew), '--', xnew, f3(xnew), '.-')\n", + "plt.legend(['data', 'linear', 'cubic', 'nearest'], loc='best')\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "id": "54318a6d", + "metadata": {}, + "source": [ + "## nD data\n", + "\n", + "There are fewer approaches to n-dimensional data, the evaluation for arbitrary dimensions is always for points on an n dimensional grid." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "712947b8", + "metadata": {}, + "outputs": [], + "source": [ + "from scipy.interpolate import griddata\n", + "\n", + "def func(x, y):\n", + " return x*(1-x)*np.cos(4*np.pi*x) * np.sin(4*np.pi*y**2)**2" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "785ac641", + "metadata": { + "solution": "hidden", + "solution_first": true + }, + "outputs": [], + "source": [ + "# A regular grid array of x,y coordinates\n", + "\n", + "grid_x, grid_y = np.mgrid[0:1:100j, 0:1:200j] # see np.info(np.mgrid) for an explanation of the 200j !!" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "9d931d1c", + "metadata": { + "solution": "hidden" + }, + "outputs": [], + "source": [ + "np.info(np.mgrid)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "97bef2a3", + "metadata": {}, + "outputs": [], + "source": [ + "# A random sampling within the same area\n", + "\n", + "points = np.random.rand(1000, 2)\n", + "values = func(points[:,0], points[:,1])\n", + "\n", + "# Resample from the values at these points onto the regular mesh\n", + "\n", + "grid_z0 = griddata(points, values, (grid_x, grid_y), method='nearest')\n", + "grid_z1 = griddata(points, values, (grid_x, grid_y), method='linear')\n", + "grid_z2 = griddata(points, values, (grid_x, grid_y), method='cubic')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "7011fa59", + "metadata": {}, + "outputs": [], + "source": [ + "plt.subplot(221)\n", + "plt.imshow(func(grid_x, grid_y).T, extent=(0,1,0,1), origin='lower', cmap='jet')\n", + "plt.plot(points[:,0], points[:,1], 'k.', ms=1)\n", + "plt.title('Original')\n", + "plt.subplot(222)\n", + "plt.imshow(grid_z0.T, extent=(0,1,0,1), origin='lower', cmap='jet')\n", + "plt.title('Nearest')\n", + "plt.subplot(223)\n", + "plt.imshow(grid_z1.T, extent=(0,1,0,1), origin='lower', cmap='jet')\n", + "plt.title('Linear')\n", + "plt.subplot(224)\n", + "plt.imshow(grid_z2.T, extent=(0,1,0,1), origin='lower', cmap='jet')\n", + "plt.title('Cubic')\n", + "plt.gcf().set_size_inches(6, 6)\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "id": "78734ed8", + "metadata": {}, + "source": [ + "## Splines \n", + "\n", + "Which have the added benefit of giving smooth derivative information" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "2eb82317", + "metadata": {}, + "outputs": [], + "source": [ + "from scipy.interpolate import splrep, splev" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "c1386e7d", + "metadata": {}, + "outputs": [], + "source": [ + "x = np.arange(0, 2*np.pi+np.pi/4, 2*np.pi/8)\n", + "y = np.sin(x)\n", + "tck = splrep(x, y, s=0)\n", + "xnew = np.arange(0, 2*np.pi, np.pi/50)\n", + "ynew = splev(xnew, tck, der=0)\n", + "yder = splev(xnew, tck, der=1)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "782a6f32", + "metadata": {}, + "outputs": [], + "source": [ + "plt.figure()\n", + "plt.plot(x, y, 'x', xnew, ynew, xnew, np.sin(xnew), x, y, 'b')\n", + "plt.legend(['Linear', 'Cubic Spline', 'True'])\n", + "plt.axis([-0.05, 6.33, -1.05, 1.05])\n", + "plt.title('Cubic-spline interpolation')\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "9a26e2ed", + "metadata": {}, + "outputs": [], + "source": [ + "plt.figure()\n", + "plt.plot(xnew, yder, xnew, np.cos(xnew),'--')\n", + "plt.legend(['Cubic Spline', 'True'])\n", + "plt.axis([-0.05, 6.33, -1.05, 1.05])\n", + "plt.title('Derivative estimation from spline')\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "id": "74f69e58", + "metadata": {}, + "source": [ + "__2D splines__ are also available" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "5f771f58", + "metadata": {}, + "outputs": [], + "source": [ + "from scipy.interpolate import bisplrep, bisplev\n", + "\n", + "# Gridded function (at low resolution ... doesn't need to be gridded data here)\n", + "\n", + "x, y = np.mgrid[-1:1:20j, -1:1:20j]\n", + "z = (x+y) * np.exp(-6.0*(x*x+y*y))\n", + "\n", + "plt.figure()\n", + "plt.pcolor(x, y, z, cmap='jet')\n", + "plt.colorbar()\n", + "plt.title(\"Sparsely sampled function.\")\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "fa19082c", + "metadata": {}, + "outputs": [], + "source": [ + "xnew, ynew = np.mgrid[-1:1:70j, -1:1:70j]\n", + "\n", + "## Create the spline-representation object tck\n", + "\n", + "tck = bisplrep(x, y, z, s=0)\n", + "znew = bisplev(xnew[:,0], ynew[0,:], tck)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "58f8e353", + "metadata": {}, + "outputs": [], + "source": [ + "plt.figure()\n", + "plt.pcolor(xnew, ynew, znew, cmap='jet')\n", + "plt.colorbar()\n", + "plt.title(\"Interpolated function.\")\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "id": "66015890", + "metadata": {}, + "source": [ + "## See also\n", + "\n", + " - Radial basis function interpolation for scattered data in n dimensions (slow for large numbers of points): `from scipy.interpolate import Rbf`\n", + " - `scipy.ndimage` for fast interpolation operations on image-like arrays\n", + " - B-splines on regular arrays are found in the `scipy.signal` module" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "23516808", + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/Notebooks/Themes/NumpyAndScipy-IPYNB/6-ScipySpatialAndMeshing.ipynb b/Notebooks/Themes/NumpyAndScipy-IPYNB/6-ScipySpatialAndMeshing.ipynb new file mode 100644 index 0000000..2cd4daa --- /dev/null +++ b/Notebooks/Themes/NumpyAndScipy-IPYNB/6-ScipySpatialAndMeshing.ipynb @@ -0,0 +1,314 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "0126597d", + "metadata": {}, + "source": [ + "# `scipy.spatial` \n", + "\n", + "scipy.spatial can compute triangulations, Voronoi diagrams, and convex hulls of a set of points, by leveraging the `Qhull` library.\n", + "\n", + "Moreover, it contains `KDTree` implementations for nearest-neighbor point queries, and utilities for distance computations in various metrics.\n", + "\n", + "## Triangulations (qhull)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a6cf0541", + "metadata": {}, + "outputs": [], + "source": [ + "%matplotlib inline\n", + "\n", + "import numpy as np\n", + "from scipy.spatial import Delaunay, ConvexHull, Voronoi\n", + "import matplotlib.pyplot as plt\n", + "\n", + "\n", + "points = np.random.rand(30, 2) # 30 random points in 2-D\n", + "\n", + "tri = Delaunay(points)\n", + "hull = ConvexHull(points)\n", + "voronoi = Voronoi(points)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "c5a5db6e", + "metadata": {}, + "outputs": [], + "source": [ + "print (\"Neighbour triangles\\n\",tri.neighbors[0:5])\n", + "print (\"Simplices\\n\", tri.simplices[0:5])\n", + "print (\"Points\\n\", points[tri.simplices[0:5]])" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "d7788996", + "metadata": {}, + "outputs": [], + "source": [ + "from scipy.spatial import delaunay_plot_2d\n", + "delaunay_plot_2d(tri)\n", + "pass" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "14003e95", + "metadata": {}, + "outputs": [], + "source": [ + "from scipy.spatial import convex_hull_plot_2d\n", + "\n", + "convex_hull_plot_2d(hull)\n", + "pass" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "e5268af7", + "metadata": {}, + "outputs": [], + "source": [ + "from scipy.spatial import voronoi_plot_2d\n", + "\n", + "voronoi_plot_2d(voronoi)\n", + "pass" + ] + }, + { + "cell_type": "markdown", + "id": "7239dbc1", + "metadata": {}, + "source": [ + "## KDtree\n", + "\n", + "Allows very fast point to point searches." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "31af4941", + "metadata": {}, + "outputs": [], + "source": [ + "from scipy.spatial import KDTree, cKDTree" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "655bdbbd", + "metadata": {}, + "outputs": [], + "source": [ + "tree = cKDTree(points)\n", + "\n", + "print (tree.data)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "453cfc3a", + "metadata": {}, + "outputs": [], + "source": [ + "%%timeit\n", + "\n", + "tree.query((0.5,0.5))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "e17cce19", + "metadata": {}, + "outputs": [], + "source": [ + "test_points = np.random.rand(1000, 2) # 1000 random points in 2-D" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a92643f3", + "metadata": {}, + "outputs": [], + "source": [ + "%%timeit\n", + "\n", + "tree.query(test_points) " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "6b53e867", + "metadata": {}, + "outputs": [], + "source": [ + "more_points = np.random.rand(10000, 2) # 1000 random points in 2-D\n", + "\n", + "big_tree = KDTree(more_points)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "172f7215", + "metadata": {}, + "outputs": [], + "source": [ + "%%timeit\n", + "\n", + "KDTree(more_points)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "6c6a2e5f", + "metadata": {}, + "outputs": [], + "source": [ + "%%timeit\n", + "\n", + "big_tree.query(test_points) " + ] + }, + { + "cell_type": "markdown", + "id": "80a0ea5b", + "metadata": {}, + "source": [ + "## Compare this to the brute-force version\n", + "\n", + "\n", + "At what point does it make sense to use kdTree and not brute-force distance tests ?\n", + "\n", + "The brute force method takes a fixed time per sample point and a fixed cost associated with the N-neighbour distance computation (but this can be vectorised efficiently)." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "8d983118", + "metadata": {}, + "outputs": [], + "source": [ + "# Brute force version\n", + "\n", + "def brute_force_distance(pts, spt):\n", + "\n", + " d = pts - spt\n", + " d2 = d**2\n", + " distances2 = np.einsum('ij->i',d2)\n", + " \n", + " nearest = np.argsort(distances2)[0]\n", + " \n", + " return np.sqrt(distances2[nearest]), nearest\n", + "\n", + "# print np.einsum('ij->i',distances2)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a35254c0", + "metadata": {}, + "outputs": [], + "source": [ + "print (brute_force_distance(more_points, (0.0,0.0)))\n", + "print (big_tree.query((0.0,0.0)))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "1a5c34e5", + "metadata": {}, + "outputs": [], + "source": [ + "%%timeit\n", + "\n", + "brute_force_distance(points, (0.5,0.5))\n", + "brute_force_distance(points, (0.0,0.0))\n", + "brute_force_distance(points, (0.25,0.25))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "8c639bbb", + "metadata": {}, + "outputs": [], + "source": [ + "%%timeit\n", + "\n", + "tree.query(np.array([(0.5,0.5), (0.0,0.0), (0.25,0.25)]))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "14b05b48", + "metadata": {}, + "outputs": [], + "source": [ + "%%timeit\n", + "\n", + "brute_force_distance(more_points, (0.5,0.5))\n", + "# brute_force_distance(more_points, (0.0,0.0))\n", + "# brute_force_distance(more_points, (0.25,0.25))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "0b684264", + "metadata": {}, + "outputs": [], + "source": [ + "%%timeit\n", + "\n", + "big_tree.query(np.array([(0.5,0.5), (0.0,0.0), (0.25,0.25)]))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "756eda49", + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "2e5a88ce", + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/Notebooks/Themes/NumpyAndScipy-IPYNB/7-ScipyOptimize.ipynb b/Notebooks/Themes/NumpyAndScipy-IPYNB/7-ScipyOptimize.ipynb new file mode 100644 index 0000000..b0e7709 --- /dev/null +++ b/Notebooks/Themes/NumpyAndScipy-IPYNB/7-ScipyOptimize.ipynb @@ -0,0 +1,269 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "0178c079", + "metadata": {}, + "source": [ + "# `scipy.optimize`\n", + "\n", + "\n", + "See the following page for a detailed list of routines: [scipy.optimize reference docs](https://docs.scipy.org/doc/scipy/reference/optimize.html#module-scipy.optimize)\n", + "\n", + "The module contains:\n", + "\n", + " - Unconstrained and constrained minimization of multivariate scalar functions (minimize) using a variety of algorithms (e.g. BFGS, Nelder-Mead simplex, Newton Conjugate Gradient, COBYLA or SLSQP)\n", + " - Global (brute-force) optimization routines (e.g. basinhopping, differential_evolution)\n", + " - Least-squares minimization (least_squares) and curve fitting (curve_fit) algorithms\n", + " - Scalar univariate functions minimizers (minimize_scalar) and root finders (newton)\n", + " - Multivariate equation system solvers (root) using a variety of algorithms (e.g. hybrid Powell, Levenberg-Marquardt or large-scale methods such as Newton-Krylov).\n", + " \n", + "Here are a couple of useful examples for simple problems." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "19bb7201", + "metadata": {}, + "outputs": [], + "source": [ + "%matplotlib inline\n", + "import matplotlib.pyplot as plt\n", + "\n", + "import numpy as np" + ] + }, + { + "cell_type": "markdown", + "id": "7f55929f", + "metadata": {}, + "source": [ + "## Curve fitting" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "60590844", + "metadata": {}, + "outputs": [], + "source": [ + "from scipy.optimize import curve_fit" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "842c3ec3", + "metadata": {}, + "outputs": [], + "source": [ + "def func(x, a, b, c):\n", + " \"\"\"\n", + " Fitting function for the data\n", + " \"\"\"\n", + " return a * np.exp(-b * x) + c" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "522c4f34", + "metadata": {}, + "outputs": [], + "source": [ + "xdata = np.linspace(0, 4, 500)\n", + "y = func(xdata, 2.5, 1.3, 0.5)\n", + "y_noise = 0.2 * np.random.normal(size=xdata.size)\n", + "ydata = y + y_noise\n", + "plt.plot(xdata, ydata, 'b-', label='data')\n", + "pass" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "252d6397", + "metadata": {}, + "outputs": [], + "source": [ + "# Fit for the parameters a, b, c of the function func\n", + "popt, pcov = curve_fit(func, xdata, ydata)\n", + "print (\"func(x, a={}, b={}, c={})\".format(popt[0], popt[1], popt[2]))\n", + "\n", + "# Now func (x, *popt) will give the y values of the fitted function\n", + "\n", + "plt.plot(xdata, ydata, 'b-', label='data')\n", + "plt.plot(xdata, func(xdata, *popt), 'r-', label='fit', linewidth=5)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "f9d67540", + "metadata": {}, + "outputs": [], + "source": [ + "# Constraints on a, b, c\n", + "poptc, pcovc = curve_fit(func, xdata, ydata, bounds=([0.,0.,1.], [3.0, 3.0, 1.5]))\n", + "print (\"func(x, a={}, b={}, c={})\".format(poptc[0], poptc[1], poptc[2]))\n", + "\n", + "plt.plot(xdata, ydata, 'b-', label='data')\n", + "plt.plot(xdata, func(xdata, *popt), 'r-', label='fit', linewidth=3)\n", + "plt.plot(xdata, func(xdata, *poptc), 'g-', label='fit-with-bounds', linewidth=10)" + ] + }, + { + "cell_type": "markdown", + "id": "78ab7160", + "metadata": {}, + "source": [ + "## Least squares\n", + "\n", + "We can also use the general least squares machinery for a more adjustable curve fitting method" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "d4587f39", + "metadata": {}, + "outputs": [], + "source": [ + "from scipy.optimize import least_squares" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "76da4a46", + "metadata": {}, + "outputs": [], + "source": [ + "def gen_data(t, a, b, c, noise=0, n_outliers=0, random_state=0):\n", + " y = a + b * np.exp(t * c)\n", + "\n", + " rnd = np.random.RandomState(random_state)\n", + " error = noise * rnd.randn(t.size)\n", + " outliers = rnd.randint(0, t.size, n_outliers)\n", + " error[outliers] *= 10\n", + "\n", + " return y + error\n", + "\n", + "a = 0.5\n", + "b = 2.0\n", + "c = -1\n", + "t_min = 0\n", + "t_max = 10\n", + "n_points = 15\n", + "\n", + "t_train = np.linspace(t_min, t_max, n_points)\n", + "y_train = gen_data(t_train, a, b, c, noise=0.1, n_outliers=3)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "d0587cf0", + "metadata": {}, + "outputs": [], + "source": [ + "# Fitting function of three parameters in 2 variables\n", + "def fun(x, t, y):\n", + " return x[0] + x[1] * np.exp(x[2] * t) - y\n", + "\n", + "x0 = np.array([1.0, 1.0, 0.0])" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "d5983cc2", + "metadata": {}, + "outputs": [], + "source": [ + "# Solve the least squares problem\n", + "\n", + "res_lsq = least_squares(fun, x0, args=(t_train, y_train))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a719e423", + "metadata": {}, + "outputs": [], + "source": [ + "# And the lsq with two different robust loss functions. \n", + "# The parameter f_scale is set to 0.1, meaning that inlier \n", + "# residuals should not significantly exceed 0.1 (the noise level used).\n", + "\n", + "res_soft_l1 = least_squares(fun, x0, loss='soft_l1', f_scale=0.1,\n", + " args=(t_train, y_train))\n", + "\n", + "res_log = least_squares(fun, x0, loss='cauchy', f_scale=0.1,\n", + " args=(t_train, y_train))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ed4d9a08", + "metadata": {}, + "outputs": [], + "source": [ + "# Compute dense arrays of these fits to plot\n", + "\n", + "t_test = np.linspace(t_min, t_max, n_points * 10)\n", + "y_true = gen_data(t_test, a, b, c)\n", + "y_lsq = gen_data(t_test, *res_lsq.x)\n", + "y_soft_l1 = gen_data(t_test, *res_soft_l1.x)\n", + "y_log = gen_data(t_test, *res_log.x)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "d70f56f7", + "metadata": {}, + "outputs": [], + "source": [ + "plt.plot(t_train, y_train, 'o')\n", + "plt.plot(t_test, y_true, 'k', linewidth=2, label='true')\n", + "plt.plot(t_test, y_lsq, label='linear loss')\n", + "plt.plot(t_test, y_soft_l1, label='soft_l1 loss')\n", + "plt.plot(t_test, y_log, label='cauchy loss')\n", + "plt.xlabel(\"t\")\n", + "plt.ylabel(\"y\")\n", + "plt.legend()\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "56ad7cff", + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "71255e4c", + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/Notebooks/Themes/NumpyAndScipy-IPYNB/8-Numpy+ScipyMiscellany.ipynb b/Notebooks/Themes/NumpyAndScipy-IPYNB/8-Numpy+ScipyMiscellany.ipynb new file mode 100644 index 0000000..96eeefc --- /dev/null +++ b/Notebooks/Themes/NumpyAndScipy-IPYNB/8-Numpy+ScipyMiscellany.ipynb @@ -0,0 +1,160 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "32d8763b", + "metadata": { + "deletable": true, + "editable": true + }, + "source": [ + "# Data storage\n", + "\n", + "Python provides file read write and object serialisation / reconstruction (python pickle module). `numpy` provides methods for storing and retrieving structured arrays quickly and efficiently (including data compression). `scipy` provides some helper functions for common file formats such as `netcdf` and `matlab` etc etc.\n", + "\n", + "Sometimes data are hard won - and reformatting them into easily retrieved files can be a lifesaver." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "27d4405e", + "metadata": { + "deletable": true, + "editable": true + }, + "outputs": [], + "source": [ + "import numpy as np\n", + "from scipy.io import netcdf" + ] + }, + { + "cell_type": "markdown", + "id": "6c40fc05", + "metadata": { + "deletable": true, + "editable": true + }, + "source": [ + "## text-based data files\n", + "\n", + "Completely portable and human readable, text-based formats are common outputs from simple scripted programs, web searches, program logs etc. Reading and writing them is trivial, and it is easy to append information to a file. The only problem is that the conversion can be extremely slow. \n", + "\n", + "This example is taken from our mapping exercise and shows the benefit of converting data to binary formats.\n", + "\n", + "---\n", + "\n", + "```python\n", + "\n", + "# Seafloor age data and global image - data from Earthbyters\n", + "\n", + "# The data come as ascii lon / lat / age tuples with NaN for no data. \n", + "# This can be loaded with ...\n", + "\n", + "age = np.loadtxt(\"global_age_data.3.6.xyz\")\n", + "age_data = age.reshape(1801,3601,3) # I looked at the data and figured out what numbers to use\n", + "age_img = age_data[:,:,2]\n", + "\n", + "# But this is super slow, so I have just stored the Age data on the grid (1801 x 3601) which we can reconstruct easily\n", + "\n", + "datasize = (1801, 3601, 3)\n", + "age_data = np.empty(datasize)\n", + "\n", + "ages = np.load(\"global_age_data.3.6.z.npz\")[\"ageData\"]\n", + "\n", + "lats = np.linspace(90, -90, datasize[0])\n", + "lons = np.linspace(-180.0,180.0, datasize[1])\n", + "\n", + "arrlons,arrlats = np.meshgrid(lons, lats)\n", + "\n", + "age_data[...,0] = arrlons[...]\n", + "age_data[...,1] = arrlats[...]\n", + "age_data[...,2] = ages[...]\n", + "```\n", + "\n", + "\n", + "---\n", + "\n", + "The timing comparison is astonishing\n", + "\n", + "On my laptop the numpy binary file is about a million times faster to read. I cut out the lat/lon values from this file to save some space, but this would add, at most, a factor of 3 to the npz timing." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "d2506d09", + "metadata": { + "deletable": true, + "editable": true + }, + "outputs": [], + "source": [ + "%%timeit\n", + "\n", + "gadtxt = np.loadtxt(\"../../Data/Resources/global_age_data.3.6.xyz\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "0bfec973", + "metadata": { + "deletable": true, + "editable": true + }, + "outputs": [], + "source": [ + "%%timeit\n", + "\n", + "gadnpy = np.load(\"../../Data/Resources/global_age_data.3.6.z.npz\")" + ] + }, + { + "cell_type": "markdown", + "id": "f912982d", + "metadata": { + "deletable": true, + "editable": true + }, + "source": [ + "## netcdf\n", + "\n", + "Many earth-science datasets are stored in the `netcdf` format. `scipy` provides a reader for this format" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "f34a489b", + "metadata": { + "deletable": true, + "editable": true + }, + "outputs": [], + "source": [ + "nf = netcdf.netcdf_file(filename=\"../../Data/Resources/velocity_AU.nc\")\n", + "\n", + "from pprint import pprint # pretty printer for python objects\n", + "\n", + "pprint( nf.dimensions )\n", + "pprint( nf.variables )\n", + "\n", + "print (nf.variables[\"lat\"].data.shape)\n", + "print (nf.variables[\"lon\"].data.shape)\n", + "print (nf.variables['ve'].data.shape)\n", + "print (nf.variables['vn'].data.shape)" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +}