How to draw Europe map using Cartopy

We can easily draw an Africa Map using Cartopy by setting the extents to [-13, 45, 30, 70]:

ax.set_extent([-13, 45, 30, 70])

Complete code example

The code above produces the image shown above:

import cartopy.crs as ccrs
import cartopy.feature as cf
from matplotlib import pyplot as plt

proj = ccrs.Miller()
ax = plt.axes(projection=proj)
ax.set_extent([-13, 45, 30, 70])
ax.stock_img()

ax.add_feature(cf.COASTLINE, lw=2)
ax.add_feature(cf.BORDERS)
# Make figure larger
plt.gcf().set_size_inches(20, 10)

# Save figure as SVG
plt.savefig("Europe.svg")

 

Posted by Uli Köhler in Cartopy, Geography, Python

How to make Cartopy coastline or border lines thicker (line width)

By standard, Cartopy draws every feature with the same line width:

ax.add_feature(cf.COASTLINE)
ax.add_feature(cf.BORDERS)

We can easily increasing the line width by adding e.g. lw=2 to the ax.add_feature() call:

ax.add_feature(cf.COASTLINE, lw=2)
ax.add_feature(cf.BORDERS)

Complete code example

This example produces the image with a wider coast line line width as shown above

import cartopy.crs as ccrs
import cartopy.feature as cf
from matplotlib import pyplot as plt

proj = ccrs.PlateCarree()
ax = plt.axes(projection=proj)
ax.set_extent([-23, 55, -35, 40])
ax.stock_img()

ax.add_feature(cf.COASTLINE, lw=2)
ax.add_feature(cf.BORDERS)
# Make figure larger
plt.gcf().set_size_inches(20, 10)

# Save figure as SVG
plt.savefig("Africa-Standard.svg")

 

Posted by Uli Köhler in Cartopy, Geography, Python

How to draw Africa map using Cartopy

We can easily draw an Africa Map using Cartopy by setting the extents to [-23, 55, -35, 40]:

ax.set_extent([-23, 55, -35, 40])

Complete code example

The code above produces the image shown above:

import cartopy.crs as ccrs
import cartopy.feature as cf
from matplotlib import pyplot as plt

proj = ccrs.PlateCarree()
ax = plt.axes(projection=proj)
ax.set_extent([-23, 55, -35, 40])
ax.stock_img()

ax.add_feature(cf.COASTLINE, lw=2)
ax.add_feature(cf.BORDERS)
# Make figure larger
plt.gcf().set_size_inches(20, 10)

# Save figure as SVG
plt.savefig("Africa.svg")

 

Posted by Uli Köhler in Cartopy, Geography, Python

How to draw straight line between two coordinates using Cartopy

In our previous posts Minimal Geodetic example using Cartopy and How to increase Geodetic resolution / accuracy / smoothness in Cartopy we have shown how to create a geodetic line on a map. While a geodetic is defined to be the shortest line on the Earth’s surface, it is not a straight line on a map projection.

In order to plot a geodetic we use:

plt.plot([lon1, lon2], [lat1, lat2], transform=ccrs.Geodetic())

In order to plot a straight line, we need to use the same projection as we used to create the map instead of transform=ccrs.Geodetic().

For example, if we created the map using

plt.axes(projection=ccrs.PlateCarree())

we need to plot the line using transform=ccrs.PlateCarree()

plt.plot([lon1, lon2], [lat1, lat2], transform=ccrs.PlateCarree())

In order to avoid errors, I strongly recommend using just one instance of the projection and assigning it to a common variable, for example:

proj = ccrs.PlateCarree()
ax = plt.axes(projection=proj)
plt.plot([-75, 77.23], [43, 28.61], transform=proj)

Note that for reasons currently unknown to me at the moment, this only works for some projections at the moment. Using ccrs.PlateCarree() and ccrs.Miller() works, but using ccrs.Mollweide() does not work!

Complete example

This code reproduces the image shown above:

import cartopy.crs as ccrs
import cartopy.feature as cf
from matplotlib import pyplot as plt

proj = ccrs.PlateCarree()
ax = plt.axes(projection=proj)
ax.stock_img()
ax.add_feature(cf.BORDERS)
# Add straight line between two points
# Format: plot([lon1, lon2], [lat1, lat2])
plt.plot([-75, 77.23], [43, 28.61], linestyle='--',
         color='blue', linewidth=8,
         transform=proj)
# Make figure larger
plt.gcf().set_size_inches(20, 10)

# Save figure as SVG
plt.savefig("Cartopy-Straight-Line-PlateCarree.svg")

 

Posted by Uli Köhler in Cartopy, Geography, Python

How to increase Geodetic resolution / accuracy / smoothness in Cartopy

In our previous post we have detailed how to draw a geodetic in Cartopy. However, as you can see in the resulting map, the geodetic line is broken up into several clearly visible segments:

In order to fix that, we need to subclass the original projection of the plot. In this example, we’re subclassing ccrs.Mollweide():

import cartopy.crs as ccrs

class HighResMollweide(ccrs.Mollweide):
    @property
    def threshold(self): return 100.0

Note that the default threshold for the Mollweide projection is 100000.0 – you can check for yourself using print(ccrs.Mollweide().threshold)

Now we can use that projection in plt.axes():

ax = plt.axes(projection=HighResMollweide())

Complete example code

This code reproduced the high-resolution geodetic image shown above:

import cartopy.crs as ccrs
import cartopy.feature as cf
from matplotlib import pyplot as plt

class HighResMollweide(ccrs.Mollweide):
    @property
    def threshold(self): return 100.0

ax = plt.axes(projection=HighResMollweide())
ax.stock_img()
ax.add_feature(cf.BORDERS)
# Add geodetic between two points
# Format: plot([lon1, lon2], [lat1, lat2])
plt.plot([-75, 77.23], [43, 28.61],
         color='blue', linewidth=2,
         transform=ccrs.Geodetic()
         )
# Make figure larger
plt.gcf().set_size_inches(20, 10)

# Save figure as SVG
plt.savefig("Cartopy-Geodetic-HiRes.svg")

Thanks to @ajdawson on StackOverflow for the original hint on how to solve this!

Posted by Uli Köhler in Cartopy, Geography, Python

Minimal Geodetic example using Cartopy

This minimal example shows you how to plot a geodetic line between two points. It is closely based on the cartopy with matplotlib intro.

import cartopy.crs as ccrs
import cartopy.feature as cf
from matplotlib import pyplot as plt
ax = plt.axes(projection = ccrs.Mollweide())
ax.stock_img()
ax.add_feature(cf.BORDERS)
# Add geodetic between two points
# Format: plot([lon1, lon2], [lat1, lat2])
plt.plot([-75, 77.23], [43, 28.61],
         color='blue', linewidth=2,
         transform=ccrs.Geodetic()
         )

 

Complete example code

This code reproduces the image shown above:

import cartopy.crs as ccrs
import cartopy.feature as cf
ax = plt.axes(projection = ccrs.Mollweide())
ax.stock_img()
ax.add_feature(cf.BORDERS)
# Add geodetic between two points
# Format: plot([lon1, lon2], [lat1, lat2])
plt.plot([-75, 77.23], [43, 28.61],
         color='blue', linewidth=2,
         transform=ccrs.Geodetic()
         )
# Make figure larger
plt.gcf().set_size_inches(20, 10)

# Save figure as SVG
plt.savefig("Cartopy-Geodetic.svg")

 

Posted by Uli Köhler in Cartopy, Geography, Python

How to add colored background to Cartopy map

Want to get from this black and white map

to this colored map

in just one line of code? Simply use cartopy’s stock_img():

ax.stock_img()

Complete example code

import cartopy.crs as ccrs
import cartopy.feature as cf
from matplotlib import pyplot as plt
ax = plt.axes(projection = ccrs.Mollweide())
ax.stock_img()
ax.add_feature(cf.COASTLINE)
ax.add_feature(cf.BORDERS)
# Make figure larger
plt.gcf().set_size_inches(20, 10)

# Save figure as SVG
plt.savefig("Cartopy-Colored.svg")

Posted by Uli Köhler in Cartopy, Geography, Python

How to draw country borders in Cartopy

Use

import cartopy.crs as ccrs
import cartopy.feature as cf
from matplotlib import pyplot as plt
ax = plt.axes(projection = ccrs.Mercator())
# This will add borders
ax.add_feature(cf.BORDERS)

The following code shows you a minimal example of how to plot country borders (and coastlines) using cartopy:

import cartopy.crs as ccrs
import cartopy.feature as cf
ax = plt.axes(projection = ccrs.Mercator())
ax.add_feature(cf.COASTLINE)
ax.add_feature(cf.BORDERS)
# Make figure larger
plt.gcf().set_size_inches(20, 10)

# Save figure as SVG
plt.savefig("Cartopy-Borders.svg")

Posted by Uli Köhler in Cartopy, Geography, Python

Cartopy minimal example with Coastlines

This example shows you a minimal example of how to plot coastlines using cartopy:

import cartopy.crs as ccrs
import cartopy.feature as cf
from matplotlib import pyplot as plt
ax = plt.axes(projection = ccrs.Mercator())
ax.add_feature(cf.COASTLINE)
# Make figure larger
plt.gcf().set_size_inches(20, 10)

# Save figure as SVG
plt.savefig("Cartopy-Coastlines.svg")

Posted by Uli Köhler in Cartopy, Geography, Python

How to draw Africa-focused map using matplotlib basemap

Note that matplotlib basemap is deprecated in favour of cartopy !

Also see:

This code allows you to draw an map that is focused on Africa using matplotlib basemap:

my_map = Basemap(projection='ortho', lat_0=10, lon_0=13, resolution='l')

Complete example code

This code reproduces the image shown above:

from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt
my_map = Basemap(projection='ortho', lat_0=10, lon_0=13, resolution='l')
my_map.drawcoastlines(linewidth=1)
my_map.drawcountries(linewidth=0.5)

# Make plot larger
plt.gcf().set_size_inches(20, 10)
# Save to file
plt.savefig("Africa.svg")

 

Posted by Uli Köhler in Geography, Python

How to draw Europe-focused map using matplotlib basemap

Also see:

This code allows you to draw an map that is focused on Europe using matplotlib basemap:

my_map = Basemap(projection='ortho', lat_0=55, lon_0=13, resolution='l')

Complete example code

This code reproduces the image shown above:

from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt
my_map = Basemap(projection='ortho', lat_0=55, lon_0=13, resolution='l')
my_map.drawcoastlines(linewidth=1)
my_map.drawcountries(linewidth=0.5)

# Make plot larger
plt.gcf().set_size_inches(20, 10)
# Save to file
plt.savefig("Europe.svg")

 

Posted by Uli Köhler in Geography, Python

How to draw South-America-focused map using matplotlib basemap

Note that matplotlib basemap is deprecated in favour of cartopy !

Also see:

This code allows you to draw an map that is focused on South America using matplotlib basemap:

my_map = Basemap(projection='ortho', lat_0=-13, lon_0=-55, resolution='l')

Complete example code

This code reproduces the image shown above:

from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt
my_map = Basemap(projection='ortho', lat_0=-13, lon_0=-55, resolution='l')
my_map.drawcoastlines(linewidth=1)
my_map.drawcountries(linewidth=0.5)

# Make plot larger
plt.gcf().set_size_inches(20, 10)
# Save to file
plt.savefig("South America.svg")

 

Posted by Uli Köhler in Geography, Python

How to draw North-America-focused map using matplotlib basemap

Note that matplotlib basemap is deprecated in favour of cartopy !

Also see:

This code allows you to draw an map that is focused on North America using matplotlib basemap:

from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt
my_map = Basemap(projection='ortho', lat_0=45 ,lon_0=-100, resolution='l')
my_map.drawcoastlines(linewidth=1)
my_map.drawcountries(linewidth=0.5)

# Make plot larger
plt.gcf().set_size_inches(20, 10)
# Save to file
plt.savefig("North America.svg")

 

Posted by Uli Köhler in Geography, Python

How to install matplotlib basemap on Ubuntu/Debian to fix ModuleNotFoundError: No module named ‘mpl_toolkits.basemap’

Note that matplotlib basemap is deprecated in favour of cartopy !

Problem:

When running your Python script, you see this error message:

---------------------------------------------------------------------------
ModuleNotFoundError                       Traceback (most recent call last)
<ipython-input-1-b1872316c43a> in <module>
----> 1 from mpl_toolkits.basemap import Basemap

ModuleNotFoundError: No module named 'mpl_toolkits.basemap'

But installing basemap using pip install basemap doesn’t work:

ERROR: Could not find a version that satisfies the requirement basemap (from versions: none)
ERROR: No matching distribution found for basemap

Solution:

First, we need to install the libgeos library and development headers:

sudo apt -y install libgeos-dev

Now we can install basemap using

pip3 install git+https://github.com/matplotlib/basemap/

This will usually take a couple of minutes.

Now run your script again, the ModuleNotFoundError should have disappeared now.

Posted by Uli Köhler in Python

How to fix Python ModuleNotFoundError: No module named ‘osmium’

Problem:

When you running your python script, you see an error message like

Traceback (most recent call last):
  File "run.py", line 2, in <module>
    import osmium as osm
ModuleNotFoundError: No module named 'osmium

Solution:

Install pyosmium using

pip install osmium

or

pip3 install osmium

 

 

Posted by Uli Köhler in OpenStreetMap, Python

Minimal example how to read .osm.pbf file using Python & osmium

This minimal example uses the osmium python binding to read an .osm.pbf file and count the number of nodes, ways and relations.

First, we need to download a suitable dataset. For this example, we’ll be using kenya-latest.osm.pbf which you can download from Geofabrik:

wget https://download.geofabrik.de/africa/kenya-latest.osm.pbf

Now we can run the script:

#!/usr/bin/env python3
import osmium as osm

class OSMHandler(osm.SimpleHandler):
    def __init__(self):
        osm.SimpleHandler.__init__(self)
        self.node_count = 0
        self.way_count = 0
        self.relation_count = 0

    def node(self, n):
        self.node_count += 1

    def way(self, w):
        self.way_count += 1

    def relation(self, r):
        self.relation_count += 1

osmhandler = OSMHandler()
osmhandler.apply_file("kenya-latest.osm.pbf")
print(f'Number of nodes: {osmhandler.node_count}')
print(f'Number of way: {osmhandler.way_count}')
print(f'Number of relations: {osmhandler.relation_count}')

Example output:

Number of nodes: 29442019
Number of way: 3188550
Number of relations: 2225

 

Posted by Uli Köhler in OpenStreetMap, Python

How to iterate column names including the index column in pandas

Just want to iterate over your DataFrame’s column names without the index column? See How to iterate column names in pandas!

In order to iterate over all column names including the index column name, use

for column_name in [df.index.name] + list(df.columns):
    print(column_name)

For example, we can print the name of all columns including the index column of the TechOverflow time series example dataset:

import pandas as pd

# Load pre-built time series example dataset
df = pd.read_csv("https://datasets.techoverflow.net/timeseries-example.csv", parse_dates=["Timestamp"])
df.set_index("Timestamp", inplace=True)

for column_name in [df.index.name] + list(df.columns):
    print(column_name)

which will print

Timestamp
Sine
Cosine

 

Posted by Uli Köhler in pandas, Python

How to get the name of the index column in pandas

In order to get the name of the index column for a pandas DataFrame, use

df.index.name

For example, we can print the name of the index column of the TechOverflow time series example dataset:

import pandas as pd

# Load pre-built time series example dataset
df = pd.read_csv("https://datasets.techoverflow.net/timeseries-example.csv", parse_dates=["Timestamp"])
df.set_index("Timestamp", inplace=True)

print(df.index.name)

which will print

Timestamp
Posted by Uli Köhler in pandas, Python

How to iterate column names in pandas

In order to iterate column names in pandas, use

for column in df.columns:
    print(columns)

For example, for the TechOverflow pandas time series example dataset, df.columnswill be

Index(['Sine', 'Cosine'], dtype='object')

so iterating over the columns using

# Load pre-built time series example dataset
df = pd.read_csv("https://datasets.techoverflow.net/timeseries-example.csv", parse_dates=["Timestamp"])
df.set_index("Timestamp", inplace=True)

# Iterate over columns
for column in df.columns:
    print(column)

will print

Sine
Cosine

Note that df.columns will not show the index column!

Posted by Uli Köhler in pandas, Python

How to replace pandas values by NaN by threshold

When processing pandas datasets, often you need to remove values above or below a given threshold from a dataset. One way to “remove” values from a dataset is to replace them by NaN (not a number) values which are typically treated as “missing” values.

For example: In order to replace values of the xcolumn by NaNwhere the x column is< 0.75 in a DataFrame df, use this snippet:

import numpy as np

df["x"][df["x"] < -0.75] = np.nan

For example, we can run this on the TechOverflow pandas time series example dataset. The original dataset has two columns: Sine and Cosine and looks like this:

After running

df["Sine"][df["Sine"] < -0.75] = np.nan

you can see that all Sine values below 0.75 have been omitted from the plot, but all the values from the Cosine column are left unchanged:

 

Complete example code:

import pandas as pd
import numpy as np
from matplotlib import pyplot as plt
plt.style.use("ggplot")

# Load pre-built time series example dataset
df = pd.read_csv("https://datasets.techoverflow.net/timeseries-example.csv", parse_dates=["Timestamp"])
df.set_index("Timestamp", inplace=True)

# Plot original code
df.plot()
plt.savefig("TimeSeries-Original.svg")

and this is the code to plot the filtered dataset:

df[df < -0.75] = np.nan
df.plot()
plt.savefig("TimeSeries-NaN.svg")

 

Posted by Uli Köhler in Data science, pandas, Python
This website uses cookies to improve your experience. We'll assume you're ok with this, but you can opt-out if you wish. Cookie settingsACCEPTPrivacy &amp; Cookies Policy