Managing Projections of Spatial Data

 Tuseday, 13 October 2020, Malika, Dakar, Senegal

Managing Projections

In this practical tutorial, we will talk about:

1. Coordinate Reference Systems
2. Practicing coordinate reference system to loaded xlsx format dataframe

Introduction:

Coordinate Reference System is very important for spatial data, it's a also a specific feature for spatial data format. For example, the pandas dataframe doesn't need crs (abbreviation of coordinate reference system).

1. Coordinate Reference Systems (CRS)

The geometric shapes in GeoSeries or GeoDataFrame is a collection of coordinates somewhere in a space. The coordinates just tell the python how these coordinates are related on the earth. If your spatial data doesn't have a correct crs, then the visualization might be different than your expectation and can also result an ambiguity while storing the data. If your spatial data doesn't have the geometries shape that's not also a spatial data.\ Setting a projection for your data is very important, but mostly we don't need it in geopandas if we use geopandas.read_file(), because in such case, the data is mostly including automatically projection information. We can ensure that by using GeoSeries.crs or geopandasdataframe.crs.

a. Most common format of coordinate reference systems

The three following crs format are the main for spatial data crs format:

• EPSG
• WKT
• pyproj

EPSG \ EPSG or European Petroleum Survey Group is a 5 digits number that represents a particular crs. The common correct format of EPSG is:\ WGS84 Lattitude/Longitude 'EPSG 4326' is the most common crs.

We can check the crs of our data for example.\ Note: The following data is uploaded from geopandas datasets

In [6]:

import geopandas
world = geopandas.read_file(geopandas.datasets.get_path('naturalearth_lowres'))
world.crs

Out[6]:

<Geographic 2D CRS: EPSG:4326>
Name: WGS 84
Axis Info [ellipsoidal]:
- Lat[north]: Geodetic latitude (degree)
- Lon[east]: Geodetic longitude (degree)
Area of Use:
- name: World
- bounds: (-180.0, -90.0, 180.0, 90.0)
Datum: World Geodetic System 1984
- Ellipsoid: WGS 84
- Prime Meridian: Greenwich

The datum refers to the 0,0 reference for the coordinate system used in the projection.

What about a data that's not from geopandas datasets?\ let's try to read an uploaded data somewhere from my computer:

In [7]:

france = geopandas.read_file('input/gadm36_FRA.gpkg')
france.crs

Out[7]:

<Geographic 2D CRS: EPSG:4326>
Name: WGS 84
Axis Info [ellipsoidal]:
- Lat[north]: Geodetic latitude (degree)
- Lon[east]: Geodetic longitude (degree)
Area of Use:
- name: World
- bounds: (-180.0, -90.0, 180.0, 90.0)
Datum: World Geodetic System 1984
- Ellipsoid: WGS 84
- Prime Meridian: Greenwich

The crs of our local data have same crs of our geopandas dataset

WKT\ WKT or Well Known Text is a markup text format to repsent vector geometric object in map or spatial reference systems, mostly used by software like ArcMap or ESRI.\ We use WGS 84 reference coordinate system to locate a longitude/latitude coordinate. WKT can be Point, Polygon, LineString or MultiPolygon etc.\ Points(0 0)\ LineString(0 0, 1 1, 1 2)\ Polygon((0 0, 4 0, 4 4, 0 4, 0 0), (1 1, 2 1, 2 2, 1 2, 1 1))\ MultiPoint((0 0), (1 2))

In [62]:

from shapely.geometry import Point, Polygon
point = Point(-16.22639, 12.81028) 
print(point)

POINT (-16.22639 12.81028)

In [63]:

point.wkt

Out[63]:

'POINT (-16.22639 12.81028)'

In [64]:

point.coords

Out[64]:

<shapely.coords.CoordinateSequence at 0x114adaa8>

In [66]:

type(point)

Out[66]:

shapely.geometry.point.Point

In [67]:

type(point.wkt)

Out[67]:

str

That's it, wkt as defined above is a text markup format for vector geometry representation in map.

Pyproj

PROJ is "a generic coordinate transformation software that transforms geospatial coordinates from one coordinate reference system (CRS) to another. This includes cartographic projections as well as geodetic transformations".\ The .crs attribue of GeoDataFrame and GeoSeries stores the CRS information as a pyproj.crs.

In [69]:

type(world.crs)

Out[69]:

pyproj.crs.crs.CRS

In [70]:

type(gdf_ville_senegal.crs)

Out[70]:

NoneType

We see 'NoneType', because gdf_ville_senegal is nto having even a coordinate reference system, the crs in that data is still in wkt.

2.Practicing coordinate reference system to loaded xlsx format dataframe

Here, we go to our next section of this practical tutorial, we will see how to set a coordinate reference system (crs) for our loaded data

In [46]:

import pandas as pd
ville_senegal = pd.read_excel('input/villes_senegal.xlsx')
ville_senegal.head()

Out[46]:

	city	Region	Long	Lat
0	Bigno	Dakar	-16.22639	12.81028
1	Dakar	Dakar	-17.44406	14.69370
2	Daara	Louga	-15.47993	15.34844
3	Diawara	Matam	-12.54374	15.02196
4	Diofior	Fatick	-16.66667	14.18333

As i don't have geometry column here, my data is actually in pandas dataframe format

In [47]:

type(ville_senegal)

Out[47]:

pandas.core.frame.DataFrame

I will create the coordinates (i mean geometry column that contains the coordinate) from longitude and latitude of my data (ville_senegal).

In [48]:

gdf_ville_senegal = geopandas.GeoDataFrame(ville_senegal, geometry=geopandas.points_from_xy(ville_senegal['Long'], ville_senegal['Lat']))

In [49]:

gdf_ville_senegal.head()

Out[49]:

	city	Region	Long	Lat	geometry
0	Bigno	Dakar	-16.22639	12.81028	POINT (-16.22639 12.81028)
1	Dakar	Dakar	-17.44406	14.69370	POINT (-17.44406 14.69370)
2	Daara	Louga	-15.47993	15.34844	POINT (-15.47993 15.34844)
3	Diawara	Matam	-12.54374	15.02196	POINT (-12.54374 15.02196)
4	Diofior	Fatick	-16.66667	14.18333	POINT (-16.66667 14.18333)

Oooh great, let see the type of our data now

In [50]:

type(gdf_ville_senegal)

Out[50]:

geopandas.geodataframe.GeoDataFrame

Let see the crs of our data

In [52]:

gdf_ville_senegal.crs

No result or no output, but let's use print

In [32]:

print(gdf_ville_senegal.crs)

None

Specifying the crs for our data

In [74]:

gdf_ville_senegal.crs = 'EPSG:4326'

In [75]:

gdf_ville_senegal.crs

Out[75]:

<Geographic 2D CRS: EPSG:4326>
Name: WGS 84
Axis Info [ellipsoidal]:
- Lat[north]: Geodetic latitude (degree)
- Lon[east]: Geodetic longitude (degree)
Area of Use:
- name: World
- bounds: (-180.0, -90.0, 180.0, 90.0)
Datum: World Geodetic System 1984
- Ellipsoid: WGS 84
- Prime Meridian: Greenwich

That's it, we have correct reference system for our data now.

Note: We notify in our introduction that mainly when the data is read by geopandas.read_file(), mostly it contains the crs, but here, we uploaded an excel format with pandas and convert it to geopandas format then added the crs.

A small break \ Practicing Geospatial in Python, made my life again more stranger, I've studied History, civilizations and International Affairs at the University, but i'm addicted with informatics, the main reason is i did not had a chance to do that before, while now, from Coursea to open knowledge everywhere changed everything.\ Let's continue with our subject

We will plot world data with the non changed crs and see

In [76]:

# Checking the crs of world 
world.crs

Out[76]:

<Geographic 2D CRS: EPSG:4326>
Name: WGS 84
Axis Info [ellipsoidal]:
- Lat[north]: Geodetic latitude (degree)
- Lon[east]: Geodetic longitude (degree)
Area of Use:
- name: World
- bounds: (-180.0, -90.0, 180.0, 90.0)
Datum: World Geodetic System 1984
- Ellipsoid: WGS 84
- Prime Meridian: Greenwich

In [77]:

# Visualize
ax = world.plot()
ax.set_title('WGS84')

Out[77]:

Text(0.5, 1.0, 'WGS84')

In [89]:

set(world_without_antarctica.continent) # Here we removed 'Antarctica'

Out[89]:

{'Africa', 'Asia', 'Europe', 'North America', 'Oceania', 'South America'}

In [94]:

# Change the crs projection from EPSG: 4326 to EPSGE:3395 (Mercator). First we will keep the globe without Antarctica
world_without_antarctica = world[(world.name!='Antarctica') & (world.name != 'Fr. S. Antarctic Lands')]
# Reproject to Mercator after removing 'Antarctica'
world_without_antarctica = world_without_antarctica.to_crs('EPSG:3395')
ax = world_without_antarctica.plot()
ax.set_title('Mercator in World without Antarctica')

Out[94]:

Text(0.5, 1.0, 'Mercator in World without Antarctica')

In [92]:

world.plot()

Out[92]:

<AxesSubplot:>

In [95]:

world_without_antarctica.crs

Out[95]:

<Projected CRS: EPSG:3395>
Name: WGS 84 / World Mercator
Axis Info [cartesian]:
- E[east]: Easting (metre)
- N[north]: Northing (metre)
Area of Use:
- name: World - between 80°S and 84°N
- bounds: (-180.0, -80.0, 180.0, 84.0)
Coordinate Operation:
- name: World Mercator
- method: Mercator (variant A)
Datum: World Geodetic System 1984
- Ellipsoid: WGS 84
- Prime Meridian: Greenwich

Conclusion

Finally, we come to the end of our practical tutorial in Managing Projection for Spatial Data.\ We sow two important with managing projection that can be in any condition with spatial data.\ There are cases where the CRS is not recognized properly, in such cases, your data might include 'towgs84' which can give problems in recognizing the crs. Other case, might be the data have different axis order of the EPSG definition, for the solution you can refer to geopandas.org/projections.html.