W8. Working with Spatial Data.ipynb

{"nbformat":4,"nbformat_minor":0,"metadata":{"kernelspec":{"name":"python3","display_name":"Python 3"},"colab":{"name":"W8. Working with Spatial Data.ipynb","provenance":[],"collapsed_sections":[]}},"cells":[{"cell_type":"markdown","metadata":{"id":"znT1hoQwk8U_"},"source":["## Quantitative Methods 2:  Data Science and Visualisation\n","\n","## Workshop 8: Working with Spatiotemporal Data\n","In this workshop, we will work with data that information about space and time, and show different ways of presenting this data, with the goal of producing fully-fledged maps.\n","\n","### Aims:\n","\n","- Plot and summarise spatial data\n","- Create simple point maps\n","- Understand the basics of projection"]},{"cell_type":"markdown","metadata":{"id":"Gb-Q2SoEk8VA"},"source":["## Downloading the Data\n","Let's grab the data we will need this week from our course website and save it into our data folder. If you've not already created a data folder then do so using the following command. \n","\n","Don't worry if it generates an error, that means you've already got a data folder."]},{"cell_type":"code","metadata":{"id":"fHsvcU2wk8VA"},"source":["!mkdir data"],"execution_count":null,"outputs":[]},{"cell_type":"code","metadata":{"id":"eCLj9srLk8VE"},"source":["!mkdir data/wk8\n","!curl https://s3.eu-west-2.amazonaws.com/qm2/wk8/tweet_data.csv -o ./data/wk8/tweet_data.csv"],"execution_count":null,"outputs":[]},{"cell_type":"markdown","metadata":{"id":"SSy1Czvpk8VG"},"source":["`------------------------------`"]},{"cell_type":"code","metadata":{"id":"FruzOm_ck8VH"},"source":["import pandas\n","import matplotlib.pyplot as plt\n","import numpy as np\n","import pylab\n","\n","%matplotlib inline\n","plt.style.use('ggplot')\n","pylab.rcParams['figure.figsize'] = (10, 8)"],"execution_count":null,"outputs":[]},{"cell_type":"markdown","metadata":{"id":"Q90Q4ptBk8VJ"},"source":["## Point and Areal Data\n","\n","We're going to look at some *point data*, data which has a spatial location but not an extent - this can be contrasted with *areal data*, where data is reported or represented as covering, or relating to, a specific region geography. An example of this second category would be the released census geography - which, as we saw earlier in the term, is reported on a bespoke areal unit called the Output Area.\n","\n","Today we will look at point data. For the purpose, we'll be looking at some data from twitter - data which has detailed spatial position as well as time and date."]},{"cell_type":"markdown","metadata":{"id":"6mlXy13bk8VK"},"source":["## The birds sing a pretty song\n","\n","Let's start by loading in the twitter data and running head() to take a look at what the dataset contains. The dataset has information about tweeters but not the content of the tweets:"]},{"cell_type":"code","metadata":{"id":"x4BrQQIbk8VL"},"source":["data_path = \"./data/wk8/tweet_data.csv\"\n","\n","tweets = pandas.read_csv(data_path, parse_dates=[1], infer_datetime_format=True, encoding = 'latin1')\n","tweets.head()"],"execution_count":null,"outputs":[]},{"cell_type":"markdown","metadata":{"id":"Bu8G-Noqk8VQ"},"source":["## Exercise\n","What does each data column represent?"]},{"cell_type":"markdown","metadata":{"id":"AnXS0B5Gk8VR"},"source":["## Working with datetime data\n","Datetime data is a little trickier to work with - it has structure which allows the extraction of hours, minutes, seconds, and so on. For example, we can just take the time part: "]},{"cell_type":"code","metadata":{"id":"f-NRm_lCk8VR"},"source":["tweets['dateT'].dt.time.head()"],"execution_count":null,"outputs":[]},{"cell_type":"markdown","metadata":{"id":"iqcYLIPGk8VV"},"source":["## Exercise\n","What temporal extent does the data cover? How do we need to structure our approach?"]},{"cell_type":"markdown","metadata":{"id":"wzY725LNk8VW"},"source":["## Exercise\n","Create a new column in the dataframe which stores the \"minute\" component of the timestamp, and use it to create a histogram of the data over the course of an hour in five minute intervals. Make sure that your graph includes a title and labelled axes."]},{"cell_type":"markdown","metadata":{"id":"OZsXOrAAk8VX"},"source":["You can also execuate this code in one line if you don't mind seeing a lot of dots - I've not included the parameters to give this some polish -  we will leave that as an exercise."]},{"cell_type":"code","metadata":{"id":"SfXPisUOk8VX"},"source":["tweets['dateT'].dt.minute.hist()"],"execution_count":null,"outputs":[]},{"cell_type":"markdown","metadata":{"id":"G7qgj0qtk8Vc"},"source":["## My First Map\n","The beauty of this data is that the data points have x and y values, so plotting them as a scatter graph will give us out first approximation (with caveats) of a map of the data:"]},{"cell_type":"code","metadata":{"id":"Haz-K2Fvk8Vd"},"source":["tweets.plot(\n","    kind='scatter',\n","    x='Lon',\n","    y='Lat',\n","    title=\"Location of Tweets\")\n","plt.xlabel(\"Longitude [degrees]\")\n","plt.ylabel(\"Latitude [degrees]\")"],"execution_count":null,"outputs":[]},{"cell_type":"markdown","metadata":{"id":"dUgAPch0k8Vh"},"source":["Notice that we've laid out the code so it's easier to see the multiple arguments in plot() - this is just the same as:"]},{"cell_type":"code","metadata":{"id":"hEwu_E9Wk8Vi"},"source":["tweets.plot(kind='scatter',x='Lon',y='Lat', title=\"Location of Tweets\")\n","plt.xlabel(\"Longitude [degrees]\")\n","plt.ylabel(\"Latitude [degrees]\")"],"execution_count":null,"outputs":[]},{"cell_type":"markdown","metadata":{"id":"ObYF6G-mk8Vm"},"source":["## Exercise\n","Change the style of the above map using the optional arguments:\n","\n","- *alpha = *: to set the opacity - 1 being opaque and 0 being transparent. Set the transparecny so that you can see busy areas *and* individual points\n","- *color=*: to set the colo*u*r to 'red'\n","- *s=*: so set the point size to 50"]},{"cell_type":"markdown","metadata":{"id":"XGUll77Ak8Vn"},"source":["So we have something that looks a bit like a heat map, and even looks a bit Gaussian. Let's see what a histogram of this data looks like:"]},{"cell_type":"code","metadata":{"id":"7FfBw7cck8Vo"},"source":["tweets['Lat'].hist(bins = 50)\n","plt.xlabel(\"Latitude [degrees]\")\n","plt.ylabel(\"Number of Tweets\")"],"execution_count":null,"outputs":[]},{"cell_type":"code","metadata":{"id":"XKPbwX-ok8Vr"},"source":["ax = tweets['Lon'].hist(bins=50)\n","plt.xlabel(\"Longitude [degrees]\")\n","ax.set_ylabel(\"Number of Tweets\")"],"execution_count":null,"outputs":[]},{"cell_type":"markdown","metadata":{"id":"wI3Usvrdk8Vu"},"source":["## The Naivest Projection\n","\n","The above tweet plot implicitly uses the *equirectangular* projection, which maps longitude onto the x axis, and latitude onto the y axis. What is the problem with this ? \n","\n","Projection is hugely complex and mathematically fiddly - luckily, we'll be working with packages which mostly do the heavy lifting for us. It's still worth thinking about projection a bit, as the process of taking points on a sphere and translating that to a flat surface is never a perfect one.\n","\n","If  we look at the picture below, then clearly the distance between $p$ and $q$ is the length of the curve on the sphere rather than the straight line between them. The closer $p$ and $q$ are the more the distance is like a straight line, and we can use a *linear mapping* - i.e. the x coordinate is a linear function of lon, and the y axis is a linear function of latitude."]},{"cell_type":"code","metadata":{"id":"EqJ0RjpTk8Vu"},"source":["from IPython.display import Image\n","Image(\"https://s3.eu-west-2.amazonaws.com/qm2/wk8/great-circle-distance.png\")"],"execution_count":null,"outputs":[]},{"cell_type":"markdown","metadata":{"id":"lKuEBAyik8Vw"},"source":["## Extension: Not-so-great circles\n","\n","Calculating great circle distances is the \"real\" way of figuring out the distance between two points on a sphere is fairly complex. Thankfully, there's a small angle approximation. \n","\n","If the two points (1 and 2) have latitudes $\\phi_1$ and $\\phi_2$ and longitudes $\\lambda_1$ and $\\lambda_2$, then let $\\Delta\\phi = \\phi_1 - \\phi_2$ and $\\Delta\\lambda = \\lambda_1 - \\lambda_2$ , where $(\\phi_1,\\lambda_1) ,(\\phi_2,\\lambda_2)$ are two points given in (latitude, longitude).  \n","\n","If $\\Delta\\phi$ and $\\Delta\\lambda$ are small enough, you can calculate the distance $D$ with : \n","\n","$$D = R \\sqrt{(\\Delta\\phi)^2+(cos(\\bar{\\phi})\\Delta\\lambda)^2}$$\n","\n","\n","where $\\bar{\\phi}$ is the mean latitude of the two points, $\\frac{1}{2}(\\phi_1+\\phi_2)$, and R is the radius of the earth. \n","\n","How small is small enough if we want to use this approximation? Well, it depends on how much error you want to incur. But generally if the angles are much less than one radian, you'll incur small errors. Radians, you say? Yes, everything in the above equations assumes angles are expressed in radians. 1 radian is about 57 degrees, but there are more precise definitions, and python has a utility function for converting between the two.\n","\n","For reference, the errors accumulated over the size of London are tens of metres."]},{"cell_type":"markdown","metadata":{"id":"HGbLZSCik8Vx"},"source":["## British Values\n","\n","The British National Grid provides projected values in metres, so we can get by without doing projections \"on the fly\" just yet. If we plot these values, it will look pretty similar, for the reasons outlined above - over the few km of London, most projection methods are quite close to the linear mapping we've done."]},{"cell_type":"code","metadata":{"id":"PybmDEITk8Vx"},"source":["ax = tweets.plot(\n","    kind='scatter',\n","    x='OSGB_Lon',\n","    y='OSGB_Lat',\n","    title=\"Location of Tweets\")\n","ax.set_xlabel(\"projected Longitude\")\n","ax.set_ylabel(\"projected Latitude\")"],"execution_count":null,"outputs":[]},{"cell_type":"markdown","metadata":{"id":"fnnTg55Sk8V3"},"source":["## Exercise: Describing the Data\n","\n","1) Find the data centroid (lat, lon)\n","\n","2) Calculate the x, y and total extent of the data, in km (or miles). (Use the projected [OSGB] data for that.)\n","\n","Hint: use commands which capture the maximum, minimum and mean of the data - describe() is a useful one here.\n"]},{"cell_type":"markdown","metadata":{"id":"fMu-1Lywk8V4"},"source":["## Exercise\n","\n","How might we go about calculating the geographical extent which contained 95% of tweets? Assuming the distribution is Gaussian in both variables, estimate a) the latitude the limits which contain 95% of tweets, b) the longitude limits which contain 95% of tweets. Then c) and add lines showing these limits to the tweet graph and d) save the figure as an image using plt.savefig(*filename*).\n","\n","What proportion of tweets is held within this box?\n","\n","How do you think the following elements influence the above result?\n","\n","- The 2D nature of the data\n","\n","- Asymmetry of the Gaussian (i.e. if $\\sigma_x \\neq \\sigma_y$)\n","\n","- Whether the data is Gaussian!\n","\n","- What other approaches could you take with this data?"]},{"cell_type":"code","metadata":{"id":"ubROHXCbk8V5"},"source":[""],"execution_count":null,"outputs":[]},{"cell_type":"code","metadata":{"id":"rie9_47yk8V9"},"source":["Image(\"https://s3.eu-west-2.amazonaws.com/qm2/wk8/tweets.png\")"],"execution_count":null,"outputs":[]},{"cell_type":"markdown","metadata":{"id":"1oE_v9nCk8WB"},"source":["## Working in 2D\n","\n","It's clear that we can learn from 1D about how we can approach 2D data. But there are limitations, and treating 2D as two sets of 1D data doesn't work for everything. We need to find ways to carry out histograms and other aggregations in 2D - and the first of these is hexbinning."]},{"cell_type":"markdown","metadata":{"id":"kcO_9JuTk8WC"},"source":["## Hexbinning\n","\n","We can also use a [hexbin clustering](http://pandas-docs.github.io/pandas-docs-travis/visualization.html#hexagonal-bin-plot) method, which is similar to binning in a histogram, the more points we have in hexbin the warmer the color. Here, we count the number data points in each hexagons, in the same way that we count the number of data in each bin for 1D data\n","\n","Luckily, the code is very easy to execute, and requires only small changes:"]},{"cell_type":"code","metadata":{"id":"6a2lM8i3k8WD"},"source":["ax = tweets.plot(\n","    kind='hexbin',\n","    x='OSGB_Lon', y='OSGB_Lat',\n","    gridsize=50,\n","    title=\"Tweet Density (Hex Bin)\",\n","    cmap='coolwarm',\n","    )\n","plt.xlabel(\"projected Longitude\")\n","plt.ylabel(\"projected Latitude\")"],"execution_count":null,"outputs":[]},{"cell_type":"markdown","metadata":{"id":"eMRBanwjk8WG"},"source":["Possible values are: Spectral, summer, coolwarm, Wistia_r, pink_r, Set1, Set2, Set3, brg_r, Dark2, prism, PuOr_r, afmhot_r, terrain_r, PuBuGn_r, RdPu, gist_ncar_r, gist_yarg_r, Dark2_r, YlGnBu, RdYlBu, hot_r, gist_rainbow_r, gist_stern, PuBu_r, cool_r, cool, gray, copper_r, Greens_r, GnBu, gist_ncar, spring_r, gist_rainbow, gist_heat_r, Wistia, OrRd_r, CMRmap, bone, gist_stern_r, RdYlGn, Pastel2_r, spring, terrain, YlOrRd_r, Set2_r, winter_r, PuBu, RdGy_r, spectral, rainbow, flag_r, jet_r, RdPu_r, gist_yarg, BuGn, Paired_r, hsv_r, bwr, cubehelix, Greens, PRGn, gist_heat, spectral_r, Paired, hsv, Oranges_r, prism_r, Pastel2, Pastel1_r, Pastel1, gray_r, jet, Spectral_r, gnuplot2_r, gist_earth, YlGnBu_r, copper, gist_earth_r, Set3_r, OrRd, gnuplot_r, ocean_r, brg, gnuplot2, PuRd_r, bone_r, BuPu, Oranges, RdYlGn_r, PiYG, CMRmap_r, YlGn, binary_r, gist_gray_r, Accent, BuPu_r, gist_gray, flag, bwr_r, RdBu_r, BrBG, Reds, Set1_r, summer_r, GnBu_r, BrBG_r, Reds_r, RdGy, PuRd, Accent_r, Blues, autumn_r, autumn, cubehelix_r, nipy_spectral_r, ocean, PRGn_r, Greys_r, pink, binary, winter, gnuplot, RdYlBu_r, hot, YlOrBr, coolwarm_r, rainbow_r, Purples_r, PiYG_r, YlGn_r, Blues_r, YlOrBr_r, seismic, Purples, seismic_r, RdBu, Greys, BuGn_r, YlOrRd, PuOr, PuBuGn, nipy_spectral, afmhot"]},{"cell_type":"markdown","metadata":{"id":"9fIZCKNak8WH"},"source":["## Extension: Splitting Mapping Data by Time\n","\n","In the next section, we use both space and time to show different geographical distributions at different times. We'll select on index, splitting the dataset in two."]},{"cell_type":"code","metadata":{"id":"GDrnlw3Mk8WH"},"source":["early = tweets[:750]\n","late = tweets[750:len(tweets)]"],"execution_count":null,"outputs":[]},{"cell_type":"code","metadata":{"id":"Z-tHtta_k8WK"},"source":["early.head()"],"execution_count":null,"outputs":[]},{"cell_type":"code","metadata":{"id":"MsLNYICnk8WM"},"source":["late.head()"],"execution_count":null,"outputs":[]},{"cell_type":"markdown","metadata":{"id":"qFXrA6xok8WP"},"source":["## Exercise\n","Plot both sets of tweets onto the same axes so they can be compared. Try and make your plot look like the image below."]},{"cell_type":"code","metadata":{"id":"mCd1DETGk8WP"},"source":["Image(\"https://s3.eu-west-2.amazonaws.com/qm2/wk8/two_times.png\")"],"execution_count":null,"outputs":[]},{"cell_type":"code","metadata":{"id":"bx8UUKrmk8WS"},"source":[""],"execution_count":null,"outputs":[]},{"cell_type":"markdown","metadata":{"id":"1YVR82sJk8WU"},"source":["We can visually inspect the spatial plots of the two time frames using hexbin plots; in this case there's not much to see..."]},{"cell_type":"code","metadata":{"id":"FoJZapDFk8WV"},"source":["ax = early.plot(\n","    kind='hexbin',\n","    x='OSGB_Lon', \n","    y='OSGB_Lat',\n","    gridsize=50,\n","    title=\"Tweet Density (Hex Bin)\",\n","    cmap='coolwarm',\n","    )\n","plt.xlabel(\"projected Longitude\")\n","plt.ylabel(\"projected Latitude\")"],"execution_count":null,"outputs":[]},{"cell_type":"code","metadata":{"id":"GlPmW4S_k8WY"},"source":["ax = late.plot(\n","    kind='hexbin',\n","    x='OSGB_Lon', \n","    y='OSGB_Lat',\n","    gridsize=50,\n","    title=\"Tweet Density (Hex Bin)\",\n","    cmap='coolwarm',\n","    )\n","plt.xlabel(\"projected Longitude\")\n","plt.ylabel(\"projected Latitude\")"],"execution_count":null,"outputs":[]},{"cell_type":"code","metadata":{"id":"d2uBZJNek8Wa"},"source":[""],"execution_count":null,"outputs":[]}]}