Demo: GeoPandas and CRS

Demo: GeoPandas and CRS#

UW Geospatial Data Analysis
CEE467/CEWA567
David Shean

GeoPandas Background#

https://geopandas.org/en/stable/getting_started/introduction.html
https://geopandas.org/data_structures.html

Key modules and packages#

Multiple levels of open-source

GeoPandas - high-level vector processing (most of what we’ll do), with dependencies:
- shapely (Python interface to GEOS) https://shapely.readthedocs.io/en/latest/manual.html
  - GEOS https://libgeos.org/
  - Handles geometry, spatial operations
- fiona (interface to GDAL/OGR) https://fiona.readthedocs.io/en/latest/README.html
  - GDAL/OGR https://gdal.org/
  - File input and output
- pyproj (interface to PROJ)
  - PROJ https://proj.org/
  - Cartographic projections and coordinate transformations library
- numpy
- pandas

import pandas as pd
import geopandas as gpd
import matplotlib.pyplot as plt

Load csv as Pandas DataFrame#

glas_fn = '../02_NumPy_Pandas_Matplotlib/data/GLAH14_tllz_conus_lulcfilt_demfilt.csv'

glas_df = pd.read_csv(glas_fn)

type(glas_df)

pandas.core.frame.DataFrame

glas_df.head()

	decyear	ordinal	lat	lon	glas_z	dem_z	dem_z_std	lulc
0	2003.139571	731266.943345	44.157897	-105.356562	1398.51	1400.52	0.33	31
1	2003.139571	731266.943346	44.150175	-105.358116	1387.11	1384.64	0.43	31
2	2003.139571	731266.943347	44.148632	-105.358427	1392.83	1383.49	0.28	31
3	2003.139571	731266.943347	44.147087	-105.358738	1384.24	1382.85	0.84	31
4	2003.139571	731266.943347	44.145542	-105.359048	1369.21	1380.24	1.73	31

Convert to GeoDataFrame#

#gpd.GeoDataFrame?

gpd.GeoDataFrame(glas_df)

	decyear	ordinal	lat	lon	glas_z	dem_z	dem_z_std	lulc
0	2003.139571	731266.943345	44.157897	-105.356562	1398.51	1400.52	0.33	31
1	2003.139571	731266.943346	44.150175	-105.358116	1387.11	1384.64	0.43	31
2	2003.139571	731266.943347	44.148632	-105.358427	1392.83	1383.49	0.28	31
3	2003.139571	731266.943347	44.147087	-105.358738	1384.24	1382.85	0.84	31
4	2003.139571	731266.943347	44.145542	-105.359048	1369.21	1380.24	1.73	31
...	...	...	...	...	...	...	...	...
65231	2009.775995	733691.238340	37.896222	-117.044399	1556.16	1556.43	0.00	31
65232	2009.775995	733691.238340	37.897769	-117.044675	1556.02	1556.43	0.00	31
65233	2009.775995	733691.238340	37.899319	-117.044952	1556.19	1556.44	0.00	31
65234	2009.775995	733691.238340	37.900869	-117.045230	1556.18	1556.44	0.00	31
65235	2009.775995	733691.238341	37.902420	-117.045508	1556.32	1556.44	0.00	31

65236 rows × 8 columns

Looks the same, let’s add geometry column!#

#gpd.points_from_xy?

mygeometry_array = gpd.points_from_xy(glas_df['lon'], glas_df['lat'])

mygeometry_array

<GeometryArray>
[<POINT (-105.357 44.158)>,  <POINT (-105.358 44.15)>,
 <POINT (-105.358 44.149)>, <POINT (-105.359 44.147)>,
 <POINT (-105.359 44.146)>, <POINT (-105.359 44.144)>,
 <POINT (-105.363 44.127)>, <POINT (-105.374 44.074)>,
 <POINT (-105.374 44.073)>, <POINT (-105.374 44.071)>,
 ...
 <POINT (-117.043 37.888)>,  <POINT (-117.043 37.89)>,
 <POINT (-117.044 37.892)>, <POINT (-117.044 37.893)>,
 <POINT (-117.044 37.895)>, <POINT (-117.044 37.896)>,
 <POINT (-117.045 37.898)>, <POINT (-117.045 37.899)>,
 <POINT (-117.045 37.901)>, <POINT (-117.046 37.902)>]
Length: 65236, dtype: geometry

A quick look at Point objects#

type(mygeometry_array[0])

shapely.geometry.point.Point

mygeometry_array[0]

../../_images/2e85d01b7a7ea09b598c0ab1f33e2886ed9f774aefeb26031bb9b933e97bbcfd.svg

p = mygeometry_array[0]

print(p)

POINT (-105.356562 44.157897)

p.x

-105.356562

p.y

44.157897

p.coords

<shapely.coords.CoordinateSequence at 0x7cdfecc8e550>

#list(p.coords)
p.coords[:]

[(-105.356562, 44.157897)]

Assign geometry column to GeoDataFrame#

glas_gdf = gpd.GeoDataFrame(glas_df, geometry=mygeometry_array)

glas_gdf.head()

	decyear	ordinal	lat	lon	glas_z	dem_z	dem_z_std	lulc	geometry
0	2003.139571	731266.943345	44.157897	-105.356562	1398.51	1400.52	0.33	31	POINT (-105.35656 44.1579)
1	2003.139571	731266.943346	44.150175	-105.358116	1387.11	1384.64	0.43	31	POINT (-105.35812 44.15018)
2	2003.139571	731266.943347	44.148632	-105.358427	1392.83	1383.49	0.28	31	POINT (-105.35843 44.14863)
3	2003.139571	731266.943347	44.147087	-105.358738	1384.24	1382.85	0.84	31	POINT (-105.35874 44.14709)
4	2003.139571	731266.943347	44.145542	-105.359048	1369.21	1380.24	1.73	31	POINT (-105.35905 44.14554)

glas_gdf.iloc[0]

decyear                        2003.139571
ordinal                      731266.943345
lat                              44.157897
lon                            -105.356562
glas_z                             1398.51
dem_z                              1400.52
dem_z_std                             0.33
lulc                                    31
geometry     POINT (-105.356562 44.157897)
Name: 0, dtype: object

glas_gdf.iloc[0].geometry

type(glas_gdf['geometry'])

geopandas.geoseries.GeoSeries

type(glas_gdf['lon'])

pandas.core.series.Series

Set CRS#

Right now, we just have x and y values for each point, but no idea what those numbers mean

glas_gdf.crs

glas_gdf.crs = 'EPSG:4326'

glas_gdf.crs

<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 ensemble
- Ellipsoid: WGS 84
- Prime Meridian: Greenwich

One-line DataFrame to GeoDataFrame constructor#

The above demonstration was interactive, step-by-step, but can do this in one shot

#glas_gdf = gpd.GeoDataFrame(glas_df, crs='EPSG:4326', geometry=gpd.points_from_xy(glas_df['lon'], glas_df['lat']))

Explore some GeoDataFrame attributes and methods#

glas_gdf.plot();

../../_images/7a3b191878a25754a3ec8e474b8aebe9dbe384efaee5c9b12e41fa40a81b653c.png

glas_gdf.plot(column='glas_z', legend=True)

<Axes: >

../../_images/4bbbed54d3b47fefb903b2ebb57369784e9f67393ff0c51fdc8e96c6e7f4d467.png

glas_gdf.total_bounds

array([-124.482406,   34.999455, -104.052336,   48.999727])

glas_gdf.geometry

       POINT (-105.35656 44.1579)
      POINT (-105.35812 44.15018)
      POINT (-105.35843 44.14863)
      POINT (-105.35874 44.14709)
      POINT (-105.35905 44.14554)
                    ...             
   POINT (-117.0444 37.89622)
  POINT (-117.04468 37.89777)
  POINT (-117.04495 37.89932)
  POINT (-117.04523 37.90087)
  POINT (-117.04551 37.90242)
Name: geometry, Length: 65236, dtype: geometry

#glas_gdf.explore()

Reproject the GeoDataFrame#

#glas_gdf.to_crs?

glas_gdf_proj = glas_gdf.to_crs('EPSG:3857')

glas_gdf_proj.head()

	decyear	ordinal	lat	lon	glas_z	dem_z	dem_z_std	lulc	geometry
0	2003.139571	731266.943345	44.157897	-105.356562	1398.51	1400.52	0.33	31	POINT (-11728238.834 5489909.71)
1	2003.139571	731266.943346	44.150175	-105.358116	1387.11	1384.64	0.43	31	POINT (-11728411.824 5488711.598)
2	2003.139571	731266.943347	44.148632	-105.358427	1392.83	1383.49	0.28	31	POINT (-11728446.444 5488472.212)
3	2003.139571	731266.943347	44.147087	-105.358738	1384.24	1382.85	0.84	31	POINT (-11728481.065 5488232.521)
4	2003.139571	731266.943347	44.145542	-105.359048	1369.21	1380.24	1.73	31	POINT (-11728515.574 5487992.837)

glas_gdf_proj.crs

<Projected CRS: EPSG:3857>
Name: WGS 84 / Pseudo-Mercator
Axis Info [cartesian]:
- X[east]: Easting (metre)
- Y[north]: Northing (metre)
Area of Use:
- name: World between 85.06°S and 85.06°N.
- bounds: (-180.0, -85.06, 180.0, 85.06)
Coordinate Operation:
- name: Popular Visualisation Pseudo-Mercator
- method: Popular Visualisation Pseudo Mercator
Datum: World Geodetic System 1984 ensemble
- Ellipsoid: WGS 84
- Prime Meridian: Greenwich

glas_gdf.plot();

glas_gdf_proj.plot();

../../_images/3b150726775e988cc42c39c0f88e09a72936791757bf0f788b04197aa989c03e.png

CRS and Projections#

There are many excellent references out there about coordinate systems and map projections. I’m not going to try to reproduce here. If you’re relatively new to all of this, please revisit resources in the reading assignment. Can also review:

https://www.axismaps.com/guide/map-projections
http://downloads2.esri.com/support/documentation/ao_/710Understanding_Map_Projections.pdf
https://courses.washington.edu/gis250/lessons/projection/

I particularly like this poster from the USGS: https://pubs.er.usgs.gov/publication/70047422

Coordinate Reference System (CRS)#

CRS is a “coordinate system with a datum”
AKA Spatial Reference System (SRS): https://en.wikipedia.org/wiki/Spatial_reference_system
Two types:
- Geodetic (AKA Geographic)
- Projected
https://datacarpentry.org/organization-geospatial/03-crs/

How to specify a CRS#

Many ways! https://geopandas.org/en/stable/docs/user_guide/projections.html

EPSG Codes
Proj strings
WKT and WKT2
- https://en.wikipedia.org/wiki/Well-known_text_representation_of_coordinate_reference_systems

!gdalsrsinfo EPSG:4326

PROJ.4 : +proj=longlat +datum=WGS84 +no_defs

OGC WKT2:2019 :
GEOGCRS["WGS 84",
    ENSEMBLE["World Geodetic System 1984 ensemble",
        MEMBER["World Geodetic System 1984 (Transit)"],
        MEMBER["World Geodetic System 1984 (G730)"],
        MEMBER["World Geodetic System 1984 (G873)"],
        MEMBER["World Geodetic System 1984 (G1150)"],
        MEMBER["World Geodetic System 1984 (G1674)"],
        MEMBER["World Geodetic System 1984 (G1762)"],
        MEMBER["World Geodetic System 1984 (G2139)"],
        MEMBER["World Geodetic System 1984 (G2296)"],
        ELLIPSOID["WGS 84",6378137,298.257223563,
            LENGTHUNIT["metre",1]],
        ENSEMBLEACCURACY[2.0]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433]],
    CS[ellipsoidal,2],
        AXIS["geodetic latitude (Lat)",north,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433]],
        AXIS["geodetic longitude (Lon)",east,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433]],
    USAGE[
        SCOPE["Horizontal component of 3D system."],
        AREA["World."],
        BBOX[-90,-180,90,180]],
    ID["EPSG",4326]]

!gdalsrsinfo EPSG:32610

PROJ.4 : +proj=utm +zone=10 +datum=WGS84 +units=m +no_defs

OGC WKT2:2019 :
PROJCRS["WGS 84 / UTM zone 10N",
    BASEGEOGCRS["WGS 84",
        ENSEMBLE["World Geodetic System 1984 ensemble",
            MEMBER["World Geodetic System 1984 (Transit)"],
            MEMBER["World Geodetic System 1984 (G730)"],
            MEMBER["World Geodetic System 1984 (G873)"],
            MEMBER["World Geodetic System 1984 (G1150)"],
            MEMBER["World Geodetic System 1984 (G1674)"],
            MEMBER["World Geodetic System 1984 (G1762)"],
            MEMBER["World Geodetic System 1984 (G2139)"],
            MEMBER["World Geodetic System 1984 (G2296)"],
            ELLIPSOID["WGS 84",6378137,298.257223563,
                LENGTHUNIT["metre",1]],
            ENSEMBLEACCURACY[2.0]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",4326]],
    CONVERSION["UTM zone 10N",
        METHOD["Transverse Mercator",
            ID["EPSG",9807]],
        PARAMETER["Latitude of natural origin",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",-123,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["Scale factor at natural origin",0.9996,
            SCALEUNIT["unity",1],
            ID["EPSG",8805]],
        PARAMETER["False easting",500000,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["(E)",east,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["(N)",north,
            ORDER[2],
            LENGTHUNIT["metre",1]],
    USAGE[
        SCOPE["Navigation and medium accuracy spatial referencing."],
        AREA["Between 126°W and 120°W, northern hemisphere between equator and 84°N, onshore and offshore. Canada - British Columbia (BC); Northwest Territories (NWT); Nunavut; Yukon. United States (USA) - Alaska (AK)."],
        BBOX[0,-126,84,-120]],
    ID["EPSG",32610]]

Map Projection check in#

https://en.wikipedia.org/wiki/Map_projection
Tradeoffs - no perfect projection, all have some distortion of distance, area, direction, shape
Infinite number of ways to represent the Earth’s 3D surface on a 2D Plane
- Plane (azimuthal), Cone, Cylinder
- Tangent, Secant (intersects sphere at two locations)
- Different parameters for different definitions
  - Center longitude, center latitude
  - Lines of “true scale”

Proj strings#

mycrs = glas_gdf.crs
mycrs

<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 ensemble
- Ellipsoid: WGS 84
- Prime Meridian: Greenwich

mycrs.to_epsg()

mycrs.to_proj4()
#Note warning

/opt/conda/lib/python3.11/site-packages/pyproj/crs/crs.py:1295: UserWarning: You will likely lose important projection information when converting to a PROJ string from another format. See: https://proj.org/faq.html#what-is-the-best-format-for-describing-coordinate-reference-systems
  proj = self._crs.to_proj4(version=version)

'+proj=longlat +datum=WGS84 +no_defs +type=crs'

glas_gdf_proj.crs.to_epsg()

glas_gdf_proj.crs.to_proj4()

'+proj=merc +a=6378137 +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null +wktext +no_defs +type=crs'

proj_str = '+proj=sinu'

glas_gdf.to_crs(proj_str)

	decyear	ordinal	lat	lon	glas_z	dem_z	dem_z_std	lulc	geometry
0	2003.139571	731266.943345	44.157897	-105.356562	1398.51	1400.52	0.33	31	POINT (-8427806.282 4891366.903)
1	2003.139571	731266.943346	44.150175	-105.358116	1387.11	1384.64	0.43	31	POINT (-8429029.663 4890508.871)
2	2003.139571	731266.943347	44.148632	-105.358427	1392.83	1383.49	0.28	31	POINT (-8429274.141 4890337.42)
3	2003.139571	731266.943347	44.147087	-105.358738	1384.24	1382.85	0.84	31	POINT (-8429518.899 4890165.747)
4	2003.139571	731266.943347	44.145542	-105.359048	1369.21	1380.24	1.73	31	POINT (-8429763.572 4889994.074)
...	...	...	...	...	...	...	...	...	...
65231	2009.775995	733691.238340	37.896222	-117.044399	1556.16	1556.43	0.00	31	POINT (-10294767.886 4195979.129)
65232	2009.775995	733691.238340	37.897769	-117.044675	1556.02	1556.43	0.00	31	POINT (-10294576.705 4196150.837)
65233	2009.775995	733691.238340	37.899319	-117.044952	1556.19	1556.44	0.00	31	POINT (-10294385.184 4196322.879)
65234	2009.775995	733691.238340	37.900869	-117.045230	1556.18	1556.44	0.00	31	POINT (-10294193.744 4196494.92)
65235	2009.775995	733691.238341	37.902420	-117.045508	1556.32	1556.44	0.00	31	POINT (-10294002.155 4196667.073)

65236 rows × 9 columns

Combining Vector Data#

#Grab the bundled world polygons
world = gpd.read_file("https://naciscdn.org/naturalearth/110m/cultural/ne_110m_admin_0_countries.zip")

world

	featurecla	scalerank	LABELRANK	SOVEREIGNT	SOV_A3	ADM0_DIF	LEVEL	TYPE	TLC	ADMIN	...	FCLASS_TR	FCLASS_ID	FCLASS_PL	FCLASS_GR	FCLASS_IT	FCLASS_NL	FCLASS_SE	FCLASS_BD	FCLASS_UA	geometry
0	Admin-0 country	1	6	Fiji	FJI	0	2	Sovereign country	1	Fiji	...	None	None	None	None	None	None	None	None	None	MULTIPOLYGON (((180 -16.06713, 180 -16.55522, ...
1	Admin-0 country	1	3	United Republic of Tanzania	TZA	0	2	Sovereign country	1	United Republic of Tanzania	...	None	None	None	None	None	None	None	None	None	POLYGON ((33.90371 -0.95, 34.07262 -1.05982, 3...
2	Admin-0 country	1	7	Western Sahara	SAH	0	2	Indeterminate	1	Western Sahara	...	Unrecognized	Unrecognized	Unrecognized	None	None	Unrecognized	None	None	None	POLYGON ((-8.66559 27.65643, -8.66512 27.58948...
3	Admin-0 country	1	2	Canada	CAN	0	2	Sovereign country	1	Canada	...	None	None	None	None	None	None	None	None	None	MULTIPOLYGON (((-122.84 49, -122.97421 49.0025...
4	Admin-0 country	1	2	United States of America	US1	1	2	Country	1	United States of America	...	None	None	None	None	None	None	None	None	None	MULTIPOLYGON (((-122.84 49, -120 49, -117.0312...
...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...
172	Admin-0 country	1	5	Republic of Serbia	SRB	0	2	Sovereign country	1	Republic of Serbia	...	None	None	None	None	None	None	None	None	None	POLYGON ((18.82982 45.90887, 18.82984 45.90888...
173	Admin-0 country	1	6	Montenegro	MNE	0	2	Sovereign country	1	Montenegro	...	None	None	None	None	None	None	None	None	None	POLYGON ((20.0707 42.58863, 19.80161 42.50009,...
174	Admin-0 country	1	6	Kosovo	KOS	0	2	Disputed	1	Kosovo	...	Admin-0 country	Unrecognized	Admin-0 country	Unrecognized	Admin-0 country	Admin-0 country	Admin-0 country	Admin-0 country	Unrecognized	POLYGON ((20.59025 41.85541, 20.52295 42.21787...
175	Admin-0 country	1	5	Trinidad and Tobago	TTO	0	2	Sovereign country	1	Trinidad and Tobago	...	None	None	None	None	None	None	None	None	None	POLYGON ((-61.68 10.76, -61.105 10.89, -60.895...
176	Admin-0 country	1	3	South Sudan	SDS	0	2	Sovereign country	1	South Sudan	...	None	None	None	None	None	None	None	None	None	POLYGON ((30.83385 3.50917, 29.9535 4.1737, 29...

177 rows × 169 columns

world.crs

<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 ensemble
- Ellipsoid: WGS 84
- Prime Meridian: Greenwich

world.geometry

    MULTIPOLYGON (((180 -16.06713, 180 -16.55522, ...
    POLYGON ((33.90371 -0.95, 34.07262 -1.05982, 3...
    POLYGON ((-8.66559 27.65643, -8.66512 27.58948...
    MULTIPOLYGON (((-122.84 49, -122.97421 49.0025...
    MULTIPOLYGON (((-122.84 49, -120 49, -117.0312...
                             ...                        
  POLYGON ((18.82982 45.90887, 18.82984 45.90888...
  POLYGON ((20.0707 42.58863, 19.80161 42.50009,...
  POLYGON ((20.59025 41.85541, 20.52295 42.21787...
  POLYGON ((-61.68 10.76, -61.105 10.89, -60.895...
  POLYGON ((30.83385 3.50917, 29.9535 4.1737, 29...
Name: geometry, Length: 177, dtype: geometry

world.plot();

../../_images/d80c4d94396adef325eb999a56a15d5d8f1d975cfcca534d0bef6e298c450151.png

Select a country by name#

world

	featurecla	scalerank	LABELRANK	SOVEREIGNT	SOV_A3	ADM0_DIF	LEVEL	TYPE	TLC	ADMIN	...	FCLASS_TR	FCLASS_ID	FCLASS_PL	FCLASS_GR	FCLASS_IT	FCLASS_NL	FCLASS_SE	FCLASS_BD	FCLASS_UA	geometry
0	Admin-0 country	1	6	Fiji	FJI	0	2	Sovereign country	1	Fiji	...	None	None	None	None	None	None	None	None	None	MULTIPOLYGON (((180 -16.06713, 180 -16.55522, ...
1	Admin-0 country	1	3	United Republic of Tanzania	TZA	0	2	Sovereign country	1	United Republic of Tanzania	...	None	None	None	None	None	None	None	None	None	POLYGON ((33.90371 -0.95, 34.07262 -1.05982, 3...
2	Admin-0 country	1	7	Western Sahara	SAH	0	2	Indeterminate	1	Western Sahara	...	Unrecognized	Unrecognized	Unrecognized	None	None	Unrecognized	None	None	None	POLYGON ((-8.66559 27.65643, -8.66512 27.58948...
3	Admin-0 country	1	2	Canada	CAN	0	2	Sovereign country	1	Canada	...	None	None	None	None	None	None	None	None	None	MULTIPOLYGON (((-122.84 49, -122.97421 49.0025...
4	Admin-0 country	1	2	United States of America	US1	1	2	Country	1	United States of America	...	None	None	None	None	None	None	None	None	None	MULTIPOLYGON (((-122.84 49, -120 49, -117.0312...
...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...
172	Admin-0 country	1	5	Republic of Serbia	SRB	0	2	Sovereign country	1	Republic of Serbia	...	None	None	None	None	None	None	None	None	None	POLYGON ((18.82982 45.90887, 18.82984 45.90888...
173	Admin-0 country	1	6	Montenegro	MNE	0	2	Sovereign country	1	Montenegro	...	None	None	None	None	None	None	None	None	None	POLYGON ((20.0707 42.58863, 19.80161 42.50009,...
174	Admin-0 country	1	6	Kosovo	KOS	0	2	Disputed	1	Kosovo	...	Admin-0 country	Unrecognized	Admin-0 country	Unrecognized	Admin-0 country	Admin-0 country	Admin-0 country	Admin-0 country	Unrecognized	POLYGON ((20.59025 41.85541, 20.52295 42.21787...
175	Admin-0 country	1	5	Trinidad and Tobago	TTO	0	2	Sovereign country	1	Trinidad and Tobago	...	None	None	None	None	None	None	None	None	None	POLYGON ((-61.68 10.76, -61.105 10.89, -60.895...
176	Admin-0 country	1	3	South Sudan	SDS	0	2	Sovereign country	1	South Sudan	...	None	None	None	None	None	None	None	None	None	POLYGON ((30.83385 3.50917, 29.9535 4.1737, 29...

177 rows × 169 columns

idx = (world['SOVEREIGNT'] == 'United States of America') & (world['TYPE'] == 'Country')

idx

    False
    False
    False
    False
     True
       ...  
  False
  False
  False
  False
  False
Length: 177, dtype: bool

us = world[idx]
us

	featurecla	scalerank	LABELRANK	SOVEREIGNT	SOV_A3	ADM0_DIF	LEVEL	TYPE	TLC	ADMIN	...	FCLASS_TR	FCLASS_ID	FCLASS_PL	FCLASS_GR	FCLASS_IT	FCLASS_NL	FCLASS_SE	FCLASS_BD	FCLASS_UA	geometry
4	Admin-0 country	1	2	United States of America	US1	1	2	Country	1	United States of America	...	None	None	None	None	None	None	None	None	None	MULTIPOLYGON (((-122.84 49, -120 49, -117.0312...

1 rows × 169 columns

us.plot()

<Axes: >

../../_images/e36f05a18530bd5ea974f350b56bbd47e01b15964f76f2a31edb02cc6c2193fb.png

us.crs

<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 ensemble
- Ellipsoid: WGS 84
- Prime Meridian: Greenwich

us.geometry

4    MULTIPOLYGON (((-122.84 49, -120 49, -117.0312...
Name: geometry, dtype: geometry

Compute area#

us.area
#Note warning!

/tmp/ipykernel_591/4249169899.py:1: UserWarning: Geometry is in a geographic CRS. Results from 'area' are likely incorrect. Use 'GeoSeries.to_crs()' to re-project geometries to a projected CRS before this operation.

  us.area

4    1122.281921
dtype: float64

Note: All calculations in GEOS/Shapely are simple euclidian geometry calculations in a 2D cartesian coordinate system!

Use an equal area projection!#

https://proj.org/operations/projections/cea.html

us_cea = us.to_crs('+proj=cea')
us_cea

	featurecla	scalerank	LABELRANK	SOVEREIGNT	SOV_A3	ADM0_DIF	LEVEL	TYPE	TLC	ADMIN	...	FCLASS_TR	FCLASS_ID	FCLASS_PL	FCLASS_GR	FCLASS_IT	FCLASS_NL	FCLASS_SE	FCLASS_BD	FCLASS_UA	geometry
4	Admin-0 country	1	2	United States of America	US1	1	2	Country	1	United States of America	...	None	None	None	None	None	None	None	None	None	MULTIPOLYGON (((-13674486.249 4793613.071, -13...

1 rows × 169 columns

us_cea.plot()

<Axes: >

../../_images/38afbc7762307100dacbdfa5856d528d0a12190f4503f34b147ba7444abfe0a0.png

us_cea.area

4    9.509851e+12
dtype: float64

us_cea.area/1E6

4    9.509851e+06
dtype: float64

us_cea = us.to_crs('+proj=tmerc +lon_0=-100')
us_cea.plot()

<Axes: >

../../_images/3c3196f1cf0b8a32a8b0b423f51b327f41ce42dacb4e0036db49bea2ecd5d6e3.png

#Back to Web Mercator
#us_cea.explore()

Combining on same plot#

#%matplotlib widget

f, ax = plt.subplots()
us.plot(ax=ax)
glas_gdf_proj.plot(ax=ax, color='r');

#Where are the US polygons?  Let's take a look...

../../_images/dd25d8c6d9ef860bdfcb832f84b471141d90837f9192c1ff753f833a0036fe70.png

Use a common projection!#

us_proj = us.to_crs(glas_gdf_proj.crs)
us_proj

	featurecla	scalerank	LABELRANK	SOVEREIGNT	SOV_A3	ADM0_DIF	LEVEL	TYPE	TLC	ADMIN	...	FCLASS_TR	FCLASS_ID	FCLASS_PL	FCLASS_GR	FCLASS_IT	FCLASS_NL	FCLASS_SE	FCLASS_BD	FCLASS_UA	geometry
4	Admin-0 country	1	2	United States of America	US1	1	2	Country	1	United States of America	...	None	None	None	None	None	None	None	None	None	MULTIPOLYGON (((-13674486.249 6274861.394, -13...

1 rows × 169 columns

us_proj.crs

<Projected CRS: EPSG:3857>
Name: WGS 84 / Pseudo-Mercator
Axis Info [cartesian]:
- X[east]: Easting (metre)
- Y[north]: Northing (metre)
Area of Use:
- name: World between 85.06°S and 85.06°N.
- bounds: (-180.0, -85.06, 180.0, 85.06)
Coordinate Operation:
- name: Popular Visualisation Pseudo-Mercator
- method: Popular Visualisation Pseudo Mercator
Datum: World Geodetic System 1984 ensemble
- Ellipsoid: WGS 84
- Prime Meridian: Greenwich

f, ax = plt.subplots()
us_proj.plot(ax=ax)
glas_gdf_proj.plot(ax=ax, color='r');

# Much better!

../../_images/859872d580306e4cb620a7e01271214815367b9fdd2b62eb984b77739e99e228.png

Demo: GeoPandas and CRS

Contents

Demo: GeoPandas and CRS#

GeoPandas Background#

Key modules and packages#

Load csv as Pandas DataFrame#

Convert to GeoDataFrame#

Looks the same, let’s add geometry column!#

A quick look at Point objects#

Assign geometry column to GeoDataFrame#

Set CRS#

One-line DataFrame to GeoDataFrame constructor#

Explore some GeoDataFrame attributes and methods#

Reproject the GeoDataFrame#

CRS and Projections#

Coordinate Reference System (CRS)#

How to specify a CRS#

Map Projection check in#

Proj strings#

Combining Vector Data#

Select a country by name#

Compute area#

Use an equal area projection!#

Combining on same plot#

Use a common projection!#