Figure 6: GeoPandas / spatial data - choropleth map - area with rivers and cities of the states of the USA #6
There are no files selected for viewing
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,35 @@ | ||
| { | ||
| "cells": [ | ||
| { | ||
| "cell_type": "code", | ||
| "execution_count": 1, | ||
| "id": "760146f6-d4c9-4e69-b539-f42e057c7945", | ||
| "metadata": {}, | ||
| "outputs": [], | ||
| "source": [ | ||
| "# Copy finale version of script from normal Python file" | ||
| ] | ||
| } | ||
| ], | ||
| "metadata": { | ||
| "kernelspec": { | ||
| "display_name": "Python 3 (ipykernel)", | ||
| "language": "python", | ||
| "name": "python3" | ||
| }, | ||
| "language_info": { | ||
| "codemirror_mode": { | ||
| "name": "ipython", | ||
| "version": 3 | ||
| }, | ||
| "file_extension": ".py", | ||
| "mimetype": "text/x-python", | ||
| "name": "python", | ||
| "nbconvert_exporter": "python", | ||
| "pygments_lexer": "ipython3", | ||
| "version": "3.12.3" | ||
| } | ||
| }, | ||
| "nbformat": 4, | ||
| "nbformat_minor": 5 | ||
| } |
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,53 @@ | ||
| import geopandas as gpd | ||
| import pygmt | ||
| # import numpy as np | ||
|
|
||
| provider = "https://naciscdn.org/naturalearth" | ||
| states = gpd.read_file(f"{provider}/50m/cultural/ne_50m_admin_1_states_provinces.zip") | ||
| rivers = gpd.read_file(f"{provider}/50m/physical/ne_50m_rivers_lake_centerlines.zip") | ||
| cities = gpd.read_file(f"{provider}/110m/cultural/ne_110m_populated_places_simple.zip") | ||
|
|
||
| states = states[states["admin"] == "United States of America"] | ||
| rivers = rivers[rivers.intersects(states.union_all())] | ||
| cities = cities[cities["adm0name"] == "United States of America"] | ||
|
|
||
| states["area_sqkm"] = states.to_crs(epsg=6933).area / 10 ** 9 | ||
| # states["area_sqkm_log"] = np.log10(states["area_sqkm"]) | ||
|
|
||
| fig = pygmt.Figure() | ||
| fig.basemap(projection="L-96/35/33/41/12c", region=[-126, -66, 25, 49], frame="+n") | ||
|
|
||
| pygmt.makecpt(cmap="hawaii", series=[0, states["area_sqkm"].max()], reverse=True) | ||
| fig.plot(data=states, cmap=True, pen="0.2p,gray50", fill="+z", aspatial="Z=area_sqkm") | ||
| fig.colorbar(frame="xaf+lArea (1000 km@+2@+)", position="jRB+o1.9c/0.2c+w3c/0.15c+ml") | ||
|
|
||
| fig.plot(data=rivers, pen="0.5p,dodgerblue4") | ||
|
|
||
| fig.plot(data=cities, style="s0.17c", fill="darkorange", pen="0.5p") | ||
| fig.text( | ||
| x=cities.geometry.x, | ||
| y=cities.geometry.y, | ||
| text=cities["name"], | ||
| offset="0.35c/0.2c", | ||
| justify="BC", | ||
| font="4.5p,Helvetica-Bold", | ||
| fill="white@30", | ||
| pen="0.2p,darkorange", | ||
| clearance="0.05c+tO", | ||
| ) | ||
|
|
||
| # Add Alaska and Hawaii separately in the lower left corner | ||
|
Member
There was a problem hiding this comment. Currently, we set the region for Alaska and Hawaii manually. Ideally, we should be able to get the region from the data directly. The tricky thing is that Hawaii has many small islands, and the # Add Alaska and Hawaii separately in the lower left corner
for name, xshift in zip(["Alaska", "Hawaii"], ["0.9c", "2.5c"]):
substates = states[(states["name"] == name)]
substates = substates.explode()
substates = substates[substates.to_crs(epsg=6933).area > 1.0e8].dissolve()
region = pygmt.info(substates, spacing=1)
with fig.shift_origin(xshift=xshift):
fig.plot(
data=substates,
region=region,
projection="M2c",
cmap=True,
pen="0.2p,gray50",
fill="+z",
aspatial="Z=area_sqkm",
) |
||
| for xshift, region in zip(["0.9c", "2.3c"], [[172, 230, 51, 72], [-168, -154, 18, 29]]): | ||
| with fig.shift_origin(xshift=xshift): | ||
| fig.plot( | ||
| data=states, | ||
| region=region, | ||
| projection="M2.5c", | ||
| cmap=True, | ||
| pen="0.2p,gray50", | ||
| fill="+z", | ||
| aspatial="Z=area_sqkm" | ||
| ) | ||
|
Comment on lines
+39
to
+50
Member
Author
There was a problem hiding this comment. Depending on how we describe this figure in the manuscript, maybe it makes sense to move the code for insets up to the code part for plotting the polygons of the main map (currently line 18).
Member
There was a problem hiding this comment. That would change the region and projection settings and we have to set region and projection again when plotting rivers and cities.
Member
Author
There was a problem hiding this comment. Yes, that's true. We do not use |
||
|
|
||
| fig.show() | ||
| fig.savefig(fname="Fig6_PyGMT_geopandas.png") | ||
Add this suggestion to a batch that can be applied as a single commit.
This suggestion is invalid because no changes were made to the code.
Suggestions cannot be applied while the pull request is closed.
Suggestions cannot be applied while viewing a subset of changes.
Only one suggestion per line can be applied in a batch.
Add this suggestion to a batch that can be applied as a single commit.
Applying suggestions on deleted lines is not supported.
You must change the existing code in this line in order to create a valid suggestion.
Outdated suggestions cannot be applied.
This suggestion has been applied or marked resolved.
Suggestions cannot be applied from pending reviews.
Suggestions cannot be applied on multi-line comments.
Suggestions cannot be applied while the pull request is queued to merge.
Suggestion cannot be applied right now. Please check back later.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
The original dataset has the
area_sqkmcolumn with zero values. That's why we compute the areas and store them in this column. But, since the area is in units of 1000 km^2, I feel the column name is more confusing. Maybe we just save it to a new column "area" instead?