<jupyter_start><jupyter_text>Introduction to Geometric ShapesThis jupyter notebook is a companion to [this](https://medium.com/@neilbardhan/location-location-location-46bba15d9b21?sk=44d4752f88733111f70c069456611fdc) post on Medium, which serves as an introduction to the basic building blocks of geospatial data. I strongly recommend going through these works in tandem as the material covered in both will help inform your knowledge of geospatial data.For this notebook, I will be using the [shapely](https://shapely.readthedocs.io/en/latest/) library for creating the geometric objects. It can be installed by running the command `pip install shapely` in your system terminal. I will also be using [mplleaflet](https://github.com/jwass/mplleaflet) to create interactive maps to display the geometric objects on a basemap to help us get a more intuitive understanding of these objects and their geometries. `mplleaflet` can be installed by running the command `pip install mplleaflet`.  Importing the libraries<jupyter_code>import os
import mplleaflet
import matplotlib.pyplot as plt
from IPython.display import Image
from IPython.display import JSON
from shapely.ops import polygonize
from shapely.geometry import (Point, LineString, Polygon, LinearRing)
os.chdir('..')<jupyter_output><empty_output><jupyter_text>Points A point is the representation of a real-world location using it's Latitude and Longitude coordinates. `Shapely` has a function called `Point` that takes as input the longitude and the latitude of a Position and returns a `POINT` object which is a tuple containing those coordinates.<jupyter_code>lat = 51.520
long = -0.093
barbican = Point([long, lat])
somewhere = Point([lat, long]) ## inherently incorrect, but serves as an illustration
print(barbican)
print(somewhere)<jupyter_output>POINT (-0.093 51.52)
POINT (51.52 -0.093)<jupyter_text>This brings up an important question, *Why is the input* `[long, lat]` *in contrast to the conventional `[lat, long]`?* That's because `Shapely`'s `Point` function resolves the earth's geometric coordinates (long/lat) to the Cartesian coordinate system, where longitude is analogous to values on the $X$ axis and latitude is analogous to values on the $Y$ axis. Therefore, a `[long, lat]` pair can be considered to be similar to a `[x, y]` pair. If we were to reverse this, we would end up in an incorrect position on the earth's surface. Let's map it out!  Plotting Points<jupyter_code>x = [barbican.x, somewhere.x]
y = [barbican.y, somewhere.y]
names = ['Barbican', 'Somewhere']

map_fig = plt.scatter(x, y)

for i, txt in enumerate(names):
    map_fig._axes.annotate(txt, (x[i] + 1, y[i] - 1))<jupyter_output><empty_output><jupyter_text>While we can use `matplotlib` for plotting our points, such plots provide no contextual information from a geographic perspective. We already know where these points lie in the coordinate space, what we would like to know is where they lie on the earth's surface.To achieve this, we almost always need the help of **basemaps**, previously designed layers containing the visual geospatial information that we are used to seeing on Google Maps or OpenStreetMaps. We can build our own basemaps, however `mplleaflet`, running on **[leaflet.js](https://leafletjs.com/examples/quick-start/)**, serves as a great ready to go tool for our mapping needs.<jupyter_code>mplleaflet.display(fig=map_fig.figure)<jupyter_output>C:\anaconda\lib\site-packages\IPython\core\display.py:694: UserWarning: Consider using IPython.display.IFrame instead
  warnings.warn("Consider using IPython.display.IFrame instead")<jupyter_text>Now this is information we can relate to! The 2 points here are - the correct point and the incorrect point. Can you guess which is the incorrect one?Accidentally interchanging the lat and the long provide us with a point off the eastern coast of Africa as opposed to a cultural center in the City of London. Here is what the correct map should look like with only the correct point plotted.<jupyter_code>map_fig = plt.scatter(barbican.x, barbican.y)
mplleaflet.display(fig=map_fig.figure)<jupyter_output><empty_output><jupyter_text>Lines A Line in it's simplest form connects two Points with a straight line. `Shapely` provides us with the `LineString` function to create `Line` geometric objects which takes as input at least two points. We have already considered The Barbican Center `[-0.093, 51.520]`, let's consider the Royal Observatory at Greenwich `[0.0, 51.476]`. A Line would simply connect these 2 Points.<jupyter_code>royal_observatory = Point([0.0, 51.476])

l1 = LineString([barbican, royal_observatory])
print(l1)<jupyter_output>LINESTRING (-0.093 51.52, 0 51.476)<jupyter_text>Plotting Lines Let's plot this line in `matplotlib` and overlay it on a basemap using `mplleaflet`. This is a non-directional, straight line connecting 2 points, The Barbican and The Royal Observatory in London.<jupyter_code>x = [barbican.x, royal_observatory.x]
y = [barbican.y, royal_observatory.y]

line_fig = plt.plot(x, y, color = 'blue')[0]
mplleaflet.display(fig=line_fig.figure, tiles = 'cartodb_positron')<jupyter_output><empty_output><jupyter_text>Polygons To help us better understand Polygons, we need to understand the concept of **Linear Rings**. A linear ring is a closed LineString with **four or more** Positions, where the first and last Positions are the same. A linear ring is the boundary of a surface or the boundary of a hole in a surface.Polygons are arrays of Linear Ring coordinates. To illustrate Polygons and Linear Rings, in addition to The Barbican `[-0.093, 51.520]` and the Royal Observatory `[0.0, 51.476]`, consider the Olympic Stadium `[-0.0165, 51.538]`. Let's assume that on our London day trip that we make a stop at the Olympic Stadium *before* heading to the Royal Observatory.<jupyter_code>## Create Point object for the Olympic Stadium
stadium = Point([-0.0165, 51.538])
## We have 3 points/Positions/Vertices for our Polygon
print(barbican)
print(royal_observatory)
print(stadium)<jupyter_output>POINT (-0.093 51.52)
POINT (0 51.476)
POINT (-0.0165 51.538)<jupyter_text>We feed in our Point coordinates to the `Polygon` function from `shapely` to get a `Polygon` object. Note that `shapely` is able to intuit that the first and last positions must be equal. Therefore, we needn't provide the last Point as an input.It is important to remember, that we must traverse our Points in clockwork order to get the correct Polygon.<jupyter_code>point_list = [barbican, stadium, royal_observatory]
poly = Polygon([[p.x, p.y] for p in point_list])
print(poly)<jupyter_output>POLYGON ((-0.093 51.52, -0.0165 51.538, 0 51.476, -0.093 51.52))<jupyter_text>Plotting Polygons We have created our Polygon, given the Points. Now let's plot our Polygon in `matplotlib` and `mplleaflet`.<jupyter_code>coords = poly.exterior.coords
pts = list(coords)
x,y = zip(*pts)
poly_fig = plt.plot(x,y)[0]
mplleaflet.display(fig=poly_fig.figure, tiles = 'cartodb_positron')<jupyter_output><empty_output><jupyter_text>Holes in Polygons Polygons are not just closed areas, but can have internal linear rings or holes nested in exterior linear rings. The image below shows what that would look like on a map.![polygon_holes](inner_ring.png) Now would be a good time to introduce **geoJSON** files. So far, we have been working with hardcoded values. We live in the age of Big Data and rarely come across values in the tens and hundreds. geoJSON files are a lightweight, JSON standard file-format used to store geometric data.Working with geoJSON files in Python is very easy. We will be using the [`geopandas`](https://geopandas.org/) library to work with geoJSON files, `geopandas` is an extension of the `pandas` library and makes working with geoJSON a breeze.  geoJSON and geopandas  Importing the librariesWe will be using `geopandas` for reading in `.geojson` files. We will also be using a web-tool from [geojson.io](https://geojson.io) called `geojsonio` to display our geometric object in the browser.<jupyter_code>import geojsonio
import geopandas as gpd<jupyter_output><empty_output><jupyter_text>Importing the dataset<jupyter_code>DATA_PATH = os.getcwd() + "//data//"
fname = "inner_ring.geojson"
poly_ring = gpd.read_file(DATA_PATH + fname)
poly_ring.head()
print(type(poly_ring))<jupyter_output><class 'geopandas.geodataframe.GeoDataFrame'><jupyter_text>`geopandas` reads in geospatial data from a multitude of formats and returns `geopandas` DataFrames. This abstraction of the geospatial data helps us perform a large variety of operations on our data by leveraging the convenience and efficiency of data frames.  Visualizing the data We can visualize geoJSON geometries in 2 ways -   1. Using `geojsonio` - This is the most convenient way as we simply pass our geometry to a 3rd party application to do all the plotting.   2. Using `geopandas` and `mplleaflet` - We use our tried and tested method of first plotting out our geometry in `matplotlib` and then overlaying the figure on a basemap.  Using `geojsonio`<jupyter_code>poly_ring_json = poly_ring.to_json()
geojsonio.display(poly_ring_json)<jupyter_output><empty_output><jupyter_text>Running the above cell will open a new tab in your browser with the geoJSON data being displayed on a map in the left hand side panel and the raw data in the right hand side panel. [geojson.io](https://www.geojson.io) is a powerful web tool, you can add additional layers, edit the geoJSON file to add more geometries, save and share your work.  Using `geopandas` and `mplleaflet` `geopandas` has a very handy `plot` function that plots the geometry of our `geopandas` DataFrame. The `plot` function returns a `matplotlib` object that we can overlay on top of our `mplleaflet` basemap.<jupyter_code>poly_fig = poly_ring.plot(alpha = 0.4, color = 'red')
mplleaflet.display(fig = poly_fig.figure, tiles = 'cartodb_positron')<jupyter_output>C:\anaconda\lib\site-packages\IPython\core\display.py:694: UserWarning: Consider using IPython.display.IFrame instead
  warnings.warn("Consider using IPython.display.IFrame instead")<jupyter_text>Multi Geometries All the geometries we have seen have been single geometries; `Points`, `Lines` and `Polygons`. Because of the complexities of real-world shapes and figures, there are certain limitations of using single shapes as often the physical world contains entities that aren't just a single contiguous shape. For example, countries like the United Kingdom and the United States have multiple disconnected parts that need multiple shapes to correctly represent one geopolitical entity. GeoJSON has versions of each of the three basic types with Multi stuck on the front: `MultiPolygon`, `MultiLineString`, `MultiPoint`. Together they give us something of a solution for multiple shapes.The way Multi Features are created is the same across all the types: **everything moves down a step of nesting**. The coordinates of the Barbican are represented as `[-0.093, 51.52]`, so a `MultiPoint` of the Barbican, the Olympic Stadium and the Royal Observatory will look like `[[-0.093, 51.520], [-0.0165, 51.538], [0.0, 51.476]]`.  Multi Points<jupyter_code>## Reading in a file containing multiple points
file_name = 'multi_points.geojson'
multi_points = gpd.read_file(DATA_PATH + file_name)
multi_points.head()<jupyter_output><empty_output><jupyter_text>A `MultiPoint` geometry contains multiple points in a single geometry and can be visualized on a basemap as shown below.<jupyter_code>mlp_fig = multi_points.plot(alpha = 1, color = 'blue')
mplleaflet.display(fig = mlp_fig.figure, tiles = 'cartodb_positron')<jupyter_output>C:\anaconda\lib\site-packages\IPython\core\display.py:694: UserWarning: Consider using IPython.display.IFrame instead
  warnings.warn("Consider using IPython.display.IFrame instead")<jupyter_text>Multi Lines<jupyter_code>## Reading in a file containing multiple lines
file_name = 'multi_line.geojson'
multi_lines = gpd.read_file(DATA_PATH + file_name)
multi_lines.head()<jupyter_output><empty_output><jupyter_text>The logic of Multiple Geometries extends to MultiLines and MultiPolygons as well. Here, we have a single geometry containing 2 different `LineStrings` as visualized below.<jupyter_code>mll_fig = multi_lines.plot(alpha = 1, color = 'blue')
mplleaflet.display(fig = mll_fig.figure, tiles = 'cartodb_positron')<jupyter_output>C:\anaconda\lib\site-packages\IPython\core\display.py:694: UserWarning: Consider using IPython.display.IFrame instead
  warnings.warn("Consider using IPython.display.IFrame instead")<jupyter_text>Multi Polygons<jupyter_code>## Reading in a file containing multiple polygons
file_name = 'hawaii.geojson'
multi_poly = gpd.read_file(DATA_PATH + file_name)
multi_poly.head()<jupyter_output><empty_output><jupyter_text>A `MultiPolygon` geometry contains multiple polygon geometries within it to accurately reflect different shapes that comprise one geopolitical entity. For example, the US state of Hawaii contains multiple islands, each is it's own polygon, but is part of a `MultiPolygon` geometry.<jupyter_code>hawaii_fig = multi_poly.plot(alpha = 0.4, color = 'red')
mplleaflet.display(fig = hawaii_fig.figure, tiles = 'cartodb_positron')<jupyter_output>C:\anaconda\lib\site-packages\IPython\core\display.py:694: UserWarning: Consider using IPython.display.IFrame instead
  warnings.warn("Consider using IPython.display.IFrame instead")<jupyter_text>FeatureCollection A Feature is the combination of the geometries and properties of a single entity in one JSON object. A FeatureCollection expands upon this by aggregating multiple features into one JSON object.<jupyter_code>## Reading in a file containing multiple polygons
file_name = 'Emergency_Services.geojson'
fc = gpd.read_file(DATA_PATH + file_name)
fc.head()<jupyter_output><empty_output><jupyter_text>A FeatureCollection is simply an array of Features. A FeatureCollection also must have Feature objects under the key `features`. The array needs to contain Feature objects only. When parsed by `geopandas` into a DataFrame, a FeatureCollection will yield multiple rows in the data frame, one for each feature. The dataframe above is a dataset containing the geometries of 4 hospitals in the London Borough of Lambeth. The source geoJSON looks as follows -<jupyter_code>print(fc.to_json(indent=2))
fc_fig = fc.plot(alpha = 1, color = 'blue')
mplleaflet.display(fig = fc_fig.figure, tiles = 'cartodb_positron')<jupyter_output>C:\anaconda\lib\site-packages\IPython\core\display.py:694: UserWarning: Consider using IPython.display.IFrame instead
  warnings.warn("Consider using IPython.display.IFrame instead")
{'fraction_non_alphanumeric': 0.08259889321582291, 'fraction_numerical': 0.01789984286397486, 'mean_word_length': 5.730114942528735, 'pattern_counts': {'":': 0, '<': 79, '<?xml version=': 0, '>': 79, 'https://': 7, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '22', 'dir': 'jupyter-structured-clean-dedup', 'id': '206064', 'license': 'MIT', 'n_tokens_mistral': 4616, 'n_tokens_neox': 4290, 'n_words': 1721, 'path': 'intro-gis/src/geometric_shapes.ipynb', 'repo_name': 'NeilBardhan/gis', 'size': '14637'}