Use geo data with Shapely

- [Instructor] Sometimes, we like to do more than just draw points on a map. For example, you might like to know if your location which is a point, is inside a city limit. Cities are represented as polygons, and the question maps to, is this point inside that polygon? There are several libraries in Python to work with what's called GIS or Geographical Information System, and geometry in general. We're going to have a look at one library called Shapely. Shapely is not installed by default, so we'll install it.

Go to the Anaconda Navigator UI. Click on environment. Make sure that the Track Environment is selected, and also that the drop down is showing Not Installed. Then search for Shapely. Once you've found it, click on Apply, and OK. Now let's get back to our notebook. The two main data structure we're going to look at are point and polygon. So from shapely dot geometry import point and polygon.

Let's create a point, point of 1, 2, and then if you look at the point, the notebook actually renders it as an image. And you can look at the point x coordinate and see what's the value. Now let's create the polygon. So poly equal polygon, and to a polygon, we give a list of points. So, 0, 0 and then 0, 10 and then 10, 10 and then 10, 0.

And now we can ask the polygon, what is the area? And we can ask it, where's the center? We can also check if a point is inside the polygon. We'll say poly dot intersects with our point. And it's True, the point is inside the polygon. Let's try another point. Poly dot intersects and we'll give it a point of 10, 20. And this time, we'll get False. The plan is to create a column of points in our data frame.

You find an area and color points that are in that area in a different color. Let's start with creating a column of points. So our mdf point will be mdf, and we'll pick just two columns, the longitude and the latitude, and we're going to apply point with axis equal 1, meaning we're going to work only on the rows. And let's take a look at the head of the data frame. And now we see that we have a column of point objects inside our data frame.

Now, let's define a polygon which is the northeast section of our trek. So first, let's calculate some points. Mean longitude and max longitude are the mdf lng dot mean and mdf lng dot lat. And the same, mean latitude and max latitude is the mdf dot lat dot mean and mdf dot lat, sorry here it should be max, dot max.

And now, let's define our polygon. The poly will be a polygon with mean longitude and mean latitude, and then mean longitude and max latitude, and then max longitude and max latitude, and the last one is the max longitude and the mean latitude. Let's plot this polygon on the map.

This is a bit tricky since folium expects a series of longitude-latitude values. If you ask the polygon for its exterior, we'll get two arrays. So, the polygon dot exterior dot xy. And we see we get two arrays. We're going to do some shape shifting here. This is quite common and lucky for us, numpy makes it easy. First we're going to mesh the two arrays together with numpy dot stack. Numpy dot stack takes two arrays and places one on top of the other. So let's import numpy, import numpy as np.

And then np dot stack poly dot exterior dot xy. And this is almost what we want. We just like to rotate the array by 90 degrees to get the pair of coordinates per row. This is called transpose and is done with the dot T attribute. So let's add a dot T here, and now we get an array of points. Let's copy the cell above, first. So we're going to copy the code above and we're also going to copy the code that adds the point to the data frame.

And now, we are going to add the polygon. So, m the map, dot add child and were going to add a folium dot Polyline with the computation that we did here, and we say that the color that we want it is yellow. And now we see our area colored in yellow. Now let's change a bit, add marker, so that the points inside the polygon will be a different color. What we're going to do is say that the color is yellow if our polygon intersects with the point in the row, otherwise, it's going to be green.

And we're going to add here, that now the fill color of the marker is going to be the color. And let's run it. And we see that we got green dots outside and yellow dots inside. I think it's pretty cool that with very few lines of code we can run maps and display GIS information. If you work a lot with GIS information, I suggest you take look at the package called GeoPandas which adds GIS data types to Pandas data frames.