This article is the first out of three of our geospatial series.

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: John Snow's investigation of the 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.

Geospatial 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.

Let'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 data
2. Geographic Reference Systems (CRS)
3. The difference between Georeferencing and Geocoding.

Vector 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:

1. Points
2. Lines. Connecting points creates a line.
3. 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.

Raster 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

To identify exact locations on the surface of the Earth, we use a geographic coordinate system.

For example, try searching for 37.971441, 23.725665 on Google Maps. Those two numbers point to an exact place - the Parthenon in Athens, Greece. The two numbers are coordinates defined by the CRS.

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:

• 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.

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.

Georeferencing 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.

Geocoding is the process of converting a human-readable address into a set of geographic coordinates. For example, the Parthenon in Athens, Greece, is at latitude 37.988579 and longitude 23.7285437.

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 article, 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. pip install shapely.

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.

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:

Shapely defines a point by its x, y coordinates, like so:

We can calculate the distance between shapely objects, such as two points:

Multiple points can be placed into a single object:

A sequence of points form a line object:

The length and bounds of a line are available with the length and bounds attributes:

A polygon is also defined by a series of points:

Polygons also have helpful attributes, such as area:

There are other useful functions where geometries interact, such as checking if the polygon pol intersects with the line from above:

We can also calculate the intersection:

But what is this object?

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.

Another 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).

Let's see this in practice.

Let'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.

If you ignore the geometry column (a shapely object), this looks like a regular dataframe. All the columns are pretty much self-explanatory.

The dataframe also includes a CRS that maps the polygons defined in the geometry column to the Earth's surface.

In our case, the CRS is EPSG:4326. That CRS uses Latitude and Longitude in degrees as coordinates.

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 need to convert to an equal-area projection that uses meters as units.

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.

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 $ 1175 \space persons / km^2$. 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!

We can call .plot() on world_gdf just like a pandas dataframe:

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. Using the parameter to_crs('epsg:4326') can do that.
2. Convert the colorbar to a logscale, which can be achieved using matplotlib.colors.LogNorm(vmin=world.pop_density.min(), vmax=world.pop_density.max())

We can pass different arguments to the plot function as you would directly on matplotlib.

Up until now, we've gone over the basics of shapely and geopandas, but now it's time we move to a complete case study.

We'll use modern Python tools to redo John Snow's analysis identifying the source of the 1854 cholera outbreak on London's Broad Street. In contrast to his Game of Thrones counterpart, London's John Snow did now something: the source of cholera. He learned it doing the first-ever geospatial analysis!

For this example, we'll use the data from Robin's blog. Robin did the work to digitize Snow's original map and data.

Let's first retrieve the data and unzip it in our current directory:

Let's see what's inside:

Ignore the file extensions for a moment, and let's see what we have here.

Vector Data

• Cholera_Deaths : number of deaths at a given spatial coordinate
• Pumps : location of water pumps

We can ignore the other files for the vector data and only deal with the '.shp' files. Called shapefiles, .shp is a standard format for vector objects.

Raster Data

• OSMap_Grayscale : raster - georeferenced grayscale map of the area from OpenStreet Maps (OSM)
• OSMap : raster - georeferenced map of the area from OpenStreet Maps (OSM)
• SnowMap : raster - digitized and georeferenced John Snow's original map

We can ignore the other files for the raster data and only deal with the '.tif' files. '.tif' is the most common format for storing raster and image data.

We've already imported geopandas and matplotlib, so all we need for the rest of the analysis is to import contextily for plotting maps. Let's import those now:

Let's read in the Cholera_Death.shp and Pumps.shp files into geopandas:

The output looks exactly like a pandas dataframe. The only difference with geopandas' dataframes is the geometry column, which is our vector dataset's essence. In our case, it includes the point coordinates of the deaths as John Snow logged them.

Let's see what the CRS data looks like:

The other difference is that correctly defined shapefiles include metadata articulating their Coordinate Reference System (CRS). In this case, it is EPSG:27700.

Let's now briefly look at the pumps data:

Similarly, pumps_df holds the positions of the water pumps near Broad Street.

Here are the pumps CRS data:

We can now plot the deaths and pumps data on a map of London's Broad Street.

We'll start building the plot by first charting deaths:

With a reference to ax, we can then plot the pumps in their locations, marking them with a red X. Let's also make the figure larger.

We would now like to show a map of London's Broad Street underneath the data. This is where we can use contextily to read the CRS data:

Now let's look at the same data, but on John Snow's original map. We can do this by changing the source parameter to SnowMap.tif, like so:

John Snow understood that most cholera deaths were clustered around a specific water pump at the intersection of Broad Street and Lexington Street (red X near the middle of the map). He attributed the outbreak to an infected water supply at that pump.

It's interesting to see how little the area has changed since 1854.

In this article, we have had a small glimpse of what you can do with geospatial data:

• We covered the basic notions that you need to understand to work with geospatial data. You should know the difference between a vector vs. raster and between geocoding vs. georeferencing. You also learned about projections, CRSs, and that Africa is HUGE!
• We covered the basics of shapely and geopandas, allowing us to work with geospatial vectors.
• Lastly, we reincarnated the first geospatial analysis. We found the infected water pump that was the source of the 1854 cholera outbreak in London.

Follow us for the following articles where we:

1. dive deeper into geopandas, preparing and analyzing a geospatial dataset
2. do machine learning with geospatial data!

After this series, you'll be ready to carry out your own spatial analysis and identify patterns in our world!

A good place to find free spatial datasets is rtwilson's list of free spatial data sources.

You can download free satellite imagery from NASA's portal or Copernicus.

Consider enrolling in a course to learn more about how to handle spatial data.