8

Using python3 and cartopy, having this code:

import matplotlib.pyplot as plt
import cartopy
import cartopy.io.shapereader as shpreader
import cartopy.crs as ccrs

ax = plt.axes(projection=ccrs.PlateCarree())
ax.add_feature(cartopy.feature.LAND)
ax.add_feature(cartopy.feature.OCEAN)
ax.add_feature(cartopy.feature.COASTLINE)
ax.add_feature(cartopy.feature.BORDERS, linestyle='-', alpha=.5)
ax.add_feature(cartopy.feature.LAKES, alpha=0.95)
ax.add_feature(cartopy.feature.RIVERS)

ax.set_extent([-150, 60, -25, 60])

shpfilename = shpreader.natural_earth(resolution='110m',
                                      category='cultural',
                                      name='admin_0_countries')

reader = shpreader.Reader(shpfilename)
countries = reader.records()

for country in countries:
    if country.attributes['SOVEREIGNT'] == "Bulgaria":
        ax.add_geometries(country.geometry, ccrs.PlateCarree(), facecolor=(0, 1, 0), label = "A")
    else:
        ax.add_geometries(country.geometry, ccrs.PlateCarree(), facecolor=(1, 1, 1), label = country.attributes['SOVEREIGNT'])
plt.rcParams["figure.figsize"] = (50,50)
plt.show()

I get this:

enter image description here

Question: What should I write, in order to get a red "A" over Bulgaria (or any other country, which I refer to in country.attributes['SOVEREIGNT'])? Currently the label is not shown at all and I am not sure how to change the font of the label. Thus, it seems that the following only changes the color, without adding the label:

ax.add_geometries(country.geometry, ccrs.PlateCarree(), facecolor=(0, 1, 0), label = "A")
Vityata
  • 39,812
  • 7
  • 40
  • 77

1 Answers1

9

You can retrieve the centroid of the geometry and plot the text at that location:

import matplotlib.patheffects as PathEffects

for country in countries:

    if country.attributes['SOVEREIGNT'] == "Bulgaria":
        g = ax.add_geometries(country.geometry, ccrs.PlateCarree(), facecolor=(0, 1, 0), label="A")

        x = country.geometry.centroid.x        
        y = country.geometry.centroid.y

        ax.text(x, y, 'A', color='red', size=15, ha='center', va='center', transform=ccrs.PlateCarree(), 
                path_effects=[PathEffects.withStroke(linewidth=5, foreground="k", alpha=.8)])

    else:
        ax.add_geometries(country.geometry, ccrs.PlateCarree(), facecolor=(1, 1, 1), label = country.attributes['SOVEREIGNT'])

With the extent focused on "Bulgaria" it looks like:

enter image description here

edit:

To get "dependencies" separate, consider using the admin_0_map_units instead of admin_0_map_countries, see the Natural Earth documentation .

To highlight small countries/regions you could add a buffer to the geometry with something like:

highlight = ['Singapore', 'Liechtenstein']

for country in countries:

    if country.attributes['NAME'] in highlight:

        if country.geometry.area < 2:
            geom = [country.geometry.buffer(2)]
        else:
            geom = [country.geometry]

        g = ax.add_geometries(geom, ccrs.PlateCarree(), facecolor=(0, 0.5, 0, 0.6), label="A", zorder=99)

        x = country.geometry.centroid.x        
        y = country.geometry.centroid.y

        ax.text(x, y+5, country.attributes['NAME'], color='red', size=14, ha='center', va='center', transform=ccrs.PlateCarree(), 
                path_effects=[PathEffects.withStroke(linewidth=3, foreground="k", alpha=.8)])

    else:
        ax.add_geometries(country.geometry, ccrs.PlateCarree(), facecolor=(1, 1, 1), label=country.attributes['NAME'])

enter image description here

You could split a specific country with something like this, It uses Shapely to perform an intersection at the middle of the geometry. Ultimately it might be "cleaner" to separate the plotting and spatial analysis (splitting etc) in to more distinct steps. Mixing it like this probably makes it harder to re-use the code for other cases.

from shapely.geometry import LineString, MultiLineString

for country in countries:

    if country.attributes['NAME'] in 'China':

        # line at the centroid y-coord of the country
        l = LineString([(-180, country.geometry.centroid.y), 
                        (180, country.geometry.centroid.y)])

        north_poly = MultiLineString([l, north_line]).convex_hull
        south_poly = MultiLineString([l, south_line]).convex_hull

        g = ax.add_geometries([country.geometry.intersection(north_poly)], ccrs.PlateCarree(), facecolor=(0.8, 0.0, 0.0, 0.4), zorder=99)
        g = ax.add_geometries([country.geometry.intersection(south_poly)], ccrs.PlateCarree(), facecolor=(0.0, 0.0, 0.8, 0.4), zorder=99)

        x = country.geometry.centroid.x        
        y = country.geometry.centroid.y

        ax.text(x, y, country.attributes['NAME'], color='k', size=16, ha='center', va='center', transform=ccrs.PlateCarree(), 
                path_effects=[PathEffects.withStroke(linewidth=5, foreground="w", alpha=1)], zorder=100)

    else:
        ax.add_geometries(country.geometry, ccrs.PlateCarree(), facecolor=(1, 1, 1), label=country.attributes['NAME'])

enter image description here

Rutger Kassies
  • 47,359
  • 12
  • 97
  • 92
  • Definitely a good answer, thank you for the efforts! And as a completely unrelated question - if I want to map Singapore or Lichtenstein in green, what should I do? – Vityata Feb 12 '20 at 13:49
  • 1
    I've made an edit to the answer. The image shows the polygon with a buffer of 2 degrees if the area is less than 2 squared degrees (not the best unit for mapping). I added an offset to the label to avoid overplotting the country itself. – Rutger Kassies Feb 12 '20 at 14:16
  • I managed to filter the dependent territories away by checking the `country.attributes['NAME_EN'] == "France"`. Anyway, concerning "Singapore" and "Lichtenstein", I guess I should try to give manual coordinates? – Vityata Feb 12 '20 at 14:41
  • 1
    Those countries are going to be very small on a large map, but I think you could still use the same concept to get it done. – Rutger Kassies Feb 12 '20 at 15:20
  • I am impressed, thank you :) Without trying to abuse your friendliness, what if I wanted to divide China into 2 parts - Northern and Southern China, is there a quick way to do it? And where exactly in the documentation are the type of units? I found this one "https://www.naturalearthdata.com/downloads/110m-cultural-vectors/", but I did not see any reference to `admin_1_map_units` or `admin_0_countries` and their usage. Although there should be some system, as https://www.naturalearthdata.com/downloads/10m-cultural-vectors/10m-admin-0-details/ presents something. – Vityata Feb 12 '20 at 15:50
  • 1
    I noticed this remark in the docs: "If you want to see the dependent overseas regions broken out (like in ISO codes, see France for example), use map units instead." and thought it might be more appropriate in this case. It's in the "about" part at: https://www.naturalearthdata.com/downloads/10m-cultural-vectors/10m-admin-0-countries/ – Rutger Kassies Feb 13 '20 at 07:42
  • Impressed again for the country splitting, based on the centroid. It definitely makes lots of sense. – Vityata Feb 13 '20 at 08:50
  • I have built upon the answer, including its examples to an article - https://www.vitoshacademy.com/python-making-maps-with-cartopy/ The working code is in GitHub here - https://github.com/Vitosh/Python_personal/blob/master/JupyterNotebook/cartopy_countries.ipynb – Vityata Feb 23 '20 at 16:14