<jupyter_start><jupyter_text>Analyze Geospatial Data in Python: GeoPandas and ShapelySource: https://www.learndatasci.com/tutorials/geospatial-data-python-geopandas-shapely/ This notebook is the learning note from the article above.Geospatial data have a lot of value. Our Geospatial series will teach you how to extract this value as a data scientist.- This 1st article introduces you to the mindset and tools needed to deal with geospatial data. It also includes a reincarnation of what has become known as the first spatial data analysis ever conducted: [<NAME>'s investigation of the 1854 Broad Street cholera outbreak](https://en.wikipedia.org/wiki/1854_Broad_Street_cholera_outbreak).- The 2nd article will dive deeper into the geospatial python framework by showing you how to conduct your own spatial analysis.- The 3rd article will apply machine learning to geospatial data.  Intro to Geospatial DataGeospatial data describe any object or feature on Earth's surface. Common examples include:- Where should a brand locate its next store?- How does the weather impact regional sales?- What's the best route to take in a car?- Which area will be hit hardest by a hurricane?- How does ice cap melting relate to carbon emissions?- Which areas will be at the highest risk of fires?Answers to these questions are valuable, making spatial data skills a great addition to any data scientist's toolset.  The BasicsLet's start by learning to speak the language of geospatial data. At the end of this section, you will know about:1. Vector vs. raster data2. Geographic Reference Systems (CRS)3. The difference between Georeferencing and Geocoding.  VectorVector data represent geometries in the world. When you open a navigation map, you see vector data. The road network, the buildings, the restaurants, and ATMs are all vectors with their associated attributes.Note: Vectors are mathematical objects. Unlike rasters, you can zoom into vectors without losing resolution.There are three main types of vector data:PointsLines. Connecting points creates a line.Polygons. Connecting lines with an enclosed area generate a polygon.  VectorVector data represent geometries in the world. When you open a navigation map, you see vector data. The road network, the buildings, the restaurants, and ATMs are all vectors with their associated attributes.*Note*: Vectors are mathematical objects. Unlike rasters, you can zoom into vectors without losing resolution.There are three main types of vector data:- Points- Lines. Connecting points creates a line.- Polygons. Connecting lines with an enclosed area generate a polygon.We can use vectors to present features and properties on the Earth’s surface. You'll most often see vectors stored in shapefiles (.shp).Specific attributes that define properties will generally accompany vectors. For example, properties of a building (e.g., its name, address, price, date built) can accompany a polygon.  RasterRaster data is a grid of pixels. Each pixel within a raster has a value, such as color, height, temperature, wind velocity, or other measurements.Whereas the default view in Google maps contains vectors, the satellite view contains raster satellite images stitched together. Each pixel in the satellite image has a value/color associated with it. Each pixel in an elevation map represents a specific height. `Raster == Image with Pixels`.These are not your usual images. They contain RGB data that our eyes can see, and multispectral or even hyperspectral information from outside the visible electromagnetic spectrum. Instead of being limited to only 3 channels/colors (RGB), we can get images with many channels.Things that are invisible to the naked eye, absorbing only a small part of the electromagnetic spectrum, can be revealed in other electromagnetic frequencies.Raster VS Vector Table|                 Vector                 |            Raster            |   ||:--------------------------------------:|:----------------------------:|---|| Points, Lines, Polygons                | Pixels                       |   || Geometric Objects, Infinitely Scalable | Fixed Grid, Fixed Resolution |   || .svg, .shp                             | .jpg, .png, .tif             |   |![](./images/raster_vs_vector.jpg)  Coordinate Reference System (CRS)To identify exact locations on the surface of the Earth, we use a geographic coordinate system.Even though the Earth is a 3-dimensional sphere, we use a 2-dimensional coordinate system of longitude (vertical lines running north-south) and latitude (horizontal lines running east-west) to identify a position on the Earth's surface. Converting a 3D sphere (the globe) into a 2D coordinate system introduces some distortions. We will explore those distortions in the next section on Map Projections. *Note: No CRS is perfect*Any choice of CRS involves a tradeoff that distorts one or all of the following:shapescale/distanceareaVery Important!!! Most mistakes in geospatial analyses come from choosing the wrong CRS for the desired operation. If you do not want to spend days and nights debugging, read this section thoroughly!![](./images/longitude-and-latitude-simple.width-1200.jpg) **Note: No CRS is perfect**Any choice of CRS involves a tradeoff that distorts one or all of the following:- shape- scale/distance- area**Very Important!!!**Most mistakes in geospatial analyses come from choosing the wrong CRS for the desired operation. If you do not want to spend days and nights debugging, read this section thoroughly! Common CRS pitfalls:Mixing coordinate systems: When combining datasets, the spatial objects MUST have the same reference system. Be sure to convert everything to the same CRS. We show you how to perform this conversion below.Calculating areas: Use an equal-area CRS before measuring a shape's area.Calculating distances: Use an equidistant CRS when calculating distances between objects.  Map Projections A map projection flattens a globe's surface by transforming coordinates from the Earth's curved surface into a flat plane.Because the Earth is not flat (I hope we agree here), any projection of the Earth into a 2D plane is a mere approximation of reality.In reality, the Earth is a geoid, meaning an irregularly-shaped ball that is not quite a sphere. The most well-known projection is the Mercator projection. As shown in the above gif, the Mercator projection inflates objects that are far from the equator.These inflations lead to some surprising revelations of our ignorance, like how the USA, China, India, and Europe all fit inside Africa.  GeoreferencingGeoreferencing is the process of assigning coordinates to vectors or rasters to project them on a model of the Earth’s surface. It is what allows us to create layers of maps.With just a click within Google Maps, you can change seamlessly from satellite view to road network view. Georeferencing makes that switch possible. GeocodingGeocoding is the process of converting a human-readable address into a set of geographic coordinates.There several libraries that handle geocoding for you. In Python, `geopandas` has a `geocoding utility` that we'll cover in the following article.In this notebook, we'll learn about `geopandas` and `shapely`, two of the most useful libraries for geospatial analysis with Python.- Shapely - a library that allows manipulation and analysis of planar geometry objects.  - Geopandas - a library that allows you to process shapefiles representing tabular data (like pandas), where every row is associated with a geometry. It provides access to many spatial functions for applying geometries, plotting maps, and geocoding. Geopandas internally uses shapely for defining geometries.  ShapelyWhat can you put into geometry?The basic shapely objects are points, lines, and polygons, but you can also define multiple objects in the same object. Then you have multipoints, multilines and multipolygons. These are useful for objects defined by various geometries, such as countries with islands.Let's see how that looks:<jupyter_code>import shapely
from shapely.geometry import Point, LineString, Polygon, MultiPoint, MultiLineString, MultiPolygon<jupyter_output><empty_output><jupyter_text>Shapely defines a point by its x, y coordinates, like so:<jupyter_code>Point(0,0)<jupyter_output><empty_output><jupyter_text>We can calculate the distance between shapely objects, such as two points:<jupyter_code>a = Point(0, 0)
b = Point(1, 0)
a.distance(b)<jupyter_output><empty_output><jupyter_text>Multiple points can be placed into a single object:<jupyter_code>MultiPoint([(0,0), (0,1), (1,1), (1,0)])<jupyter_output><empty_output><jupyter_text>A sequence of points form a line object:<jupyter_code>line = LineString([(0,0),(1,2), (0,1)])
line<jupyter_output><empty_output><jupyter_text>The length and bounds of a line are available with the length and bounds attributes:<jupyter_code>print(f'Length of line {line.length}')
print(f'Bounds of line {line.bounds}')<jupyter_output>Length of line 3.6502815398728847
Bounds of line (0.0, 0.0, 1.0, 2.0)<jupyter_text>A polygon is also defined by a series of points:<jupyter_code>pol = Polygon([(0,0), (0,1), (1,1), (1,0)])
pol<jupyter_output><empty_output><jupyter_text>Polygons also have helpful attributes, such as area:<jupyter_code>pol.area<jupyter_output><empty_output><jupyter_text>There are other useful functions where geometries interact, such as checking if the polygon `pol` intersects with the `line` from above:<jupyter_code>pol.intersects(line)<jupyter_output><empty_output><jupyter_text>We can also calculate the intersection:<jupyter_code>pol.intersection(line)<jupyter_output><empty_output><jupyter_text>But what is this object?<jupyter_code>print(pol.intersection(line))<jupyter_output>GEOMETRYCOLLECTION (POINT (0 1), LINESTRING (0 0, 0.5 1))<jupyter_text>It's a `GeometryCollection`, which is a collection of different types of geometries.Pretty straightforward and intuitive so far! You can do so much more with the shapely library, so be sure to check the [docs](https://shapely.readthedocs.io/en/stable/index.html).  Geopandas BasicsAnother tool for working with geospatial data is geopandas. As we know, pandas DataFrames represent tabular datasets. Similarly, geopandas DataFrames represent tabular data with two extensions:- The geometry column defines a point, line, or polygon associated with the rest of the columns. This column is a collection of shapely objects. Whatever you can do with shapely objects, you can also do with the geometry object.- The CRS is the coordinate reference system of the geometry column that tells us where a point, line, or polygon lies on the Earth's surface. Geopandas maps a geometry onto the Earth's surface (e.g., WGS84). The easiest way to install geopandas on Windows is to use Anaconda with the following command: ```bashconda install -c conda-forge geopandas```<jupyter_code>import matplotlib
import geopandas as gpd<jupyter_output><empty_output><jupyter_text>Loading datasetLet's start by loading a dataset shipped with geopandas, called 'naturalearth_lowres'. This dataset includes the geometry of each country in the world, accompanied by some further details such as Population and GDP estimates.<jupyter_code>world_gdf = gpd.read_file(
    gpd.datasets.get_path('naturalearth_lowres')
)
world_gdf<jupyter_output><empty_output><jupyter_text>CRSThe dataframe also includes a CRS that maps the polygons defined in the geometry column to the Earth's surface.<jupyter_code>world_gdf.crs<jupyter_output><empty_output><jupyter_text>In our case, the CRS is EPSG:4326. That CRS uses Latitude and Longitude in degrees as coordinates.**Note - Components of a CRS**- Datum - The reference system, which in our case defines the starting point of measurement (Prime Meridian) and the model of the shape of the Earth (Ellipsoid). The most common Datum is WGS84, but it is not the only one.- Area of use - In our case, the are of use is the whole world, but there are many CRS that are optimized for a particular area of interest.- Axes and Units - Usually, longitude and latitude are measured in degrees. Units for x, y coordinates are often measured in meters. Let's see an application for which we have to change the CRS.Let's measure the population density of each country! We can measure the area of each geometry but bear in mind that we need first convert to an equal-area projection that uses meters as units.<jupyter_code>world_gdf = world_gdf.to_crs("+proj=eck4 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs")
world_gdf.crs<jupyter_output><empty_output><jupyter_text>We can now calculate each country's population density by dividing the population estimate by the area. Note: We can access the area of the geometries as we would regular columns. Although no column contains geometry areas, the area is an attribute of the geometry objects.?<jupyter_code>world_gdf['pop_density'] = world_gdf.pop_est / world_gdf.area * 10**6

world_gdf.sort_values(by='pop_density', ascending=False)<jupyter_output><empty_output><jupyter_text>Notice how the geometry objects now have values that are in totally different units than before.Just looking at the dataframe above, we can quickly identify the outliers. Bangladesh has a population density of around 1174 persons/km2. Antarctica has a near-zero population density, with only 810 people living in a vast space.It's always better to visualize maps, though. So let's visualize!  VisualizationWe can call `.plot()` on `world_gdf` just like a pandas dataframe:<jupyter_code>figsize = (20, 11)
world_gdf.plot('pop_density', legend=True, figsize=figsize)<jupyter_output><empty_output><jupyter_text>The above map doesn't look very helpful, so let's make it better by doing the following:1. Change to the Mercator projection since it's more familiar.2. Convert the colorbar to a logscale.We can pass different arguments to the plot function as you would directly on matplotlib.<jupyter_code>norm_color = matplotlib.colors.LogNorm(vmin=world_gdf.pop_density.min(), vmax=world_gdf.pop_density.max(), clip=True)

world_gdf.to_crs('epsg:4326').plot("pop_density", 
                                   figsize=figsize, 
                                   legend=False, norm=norm_color)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.0649505506067195, 'fraction_numerical': 0.008557199971943607, 'mean_word_length': 4.945788156797331, 'pattern_counts': {'":': 0, '<': 72, '<?xml version=': 0, '>': 72, 'https://': 3, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '18', 'dir': 'jupyter-structured-clean-dedup', 'id': '650786', 'license': 'MIT', 'n_tokens_mistral': 3901, 'n_tokens_neox': 3646, 'n_words': 1829, 'path': '10.Geospatial_Data/Analyze_Geospatial_Data.ipynb', 'repo_name': 'cuicaihao/Data_Science_Python', 'size': '14260'}