Search code examples
pythonmatplotlibgeospatialtopography

How to change background and outline color on topographic map?

I followed this excellent guide by Adam Symington and successfully created the following topographic map of Sabah (a state in Malaysia, which is a Southeast Asian nation). The awkward blob of black in the upper left corner is my attempt to plot certain coordinates on the map.

enter image description here

I would like to improve this diagram in the following ways:

EDIT: I have figured item (1) out and posted the solution below. (2) and (3) pending.

  1. [SOLVED] The sch dataframe contains coordinates of all schools in the state. I would like to plot these on the map. I suspect that it is currently going wonky because the axes are not "geo-axes" (meaning, not using lat/lon scales) - you can confirm this by setting ax.axis('on'). How do I get around this? [SOLVED]

  2. I'd like to set the portion outside the actual territory to white. Calling ax.set_facecolor('white') isn't working. I know that the specific thing setting it to grey is the ax.imshow(hillshade, cmap='Greys', alpha=0.3) line (because changing the cmap changes the background); I just don't know how to alter it while keeping the color within the map as grey.

  3. If possible, I'd like the outline of the map to be black, but this is just pedantic.

All code to reproduce the diagram above is below. The downloadSrc function gets and saves the dependencies (a 5.7MB binary file containing the topographic data and a 0.05MB csv containing the coordinates of points to plot) in a local folder; you need only run that once.

import rasterio
from rasterio import mask as msk

import matplotlib.pyplot as plt
from matplotlib.colors import LinearSegmentedColormap
from matplotlib.colors import ListedColormap

import numpy as np
import pandas as pd
import geopandas as gpd
import earthpy.spatial as es
from shapely.geometry import Point

def downloadSrc(dl=1):
    if dl == 1:
        import os
        os.mkdir('sabah')

        import requests
        r = requests.get('https://raw.githubusercontent.com/Thevesh/Display/master/sabah_tiff.npy')
        with open('sabah/sabah_topog.npy', 'wb') as f: f.write(r.content)

        df = pd.read_csv('https://raw.githubusercontent.com/Thevesh/Display/master/schools.csv')
        df.to_csv('sabah/sabah_schools.csv')

# Set dl = 0 after first run; the files will be in your current working directory + /sabah
downloadSrc(dl=1)

# Load topography of Sabah, pre-saved from clipped tiff file (as per Adam Symington guide)
value_range = 4049
sabah_topography = np.load('sabah/sabah_topog.npy')

# Load coordinates of schools in Sabah
crs={'init':'epsg:4326'}
sch = pd.read_csv('sabah/sabah_schools.csv',usecols=['lat','lon'])
geometry = [Point(xy) for xy in zip(sch.lon, sch.lat)]
schools = gpd.GeoDataFrame(sch, crs=crs, geometry=geometry)

# Replicated directly from guide, with own modifications only to colours
sabah_colormap = LinearSegmentedColormap.from_list('sabah', ['lightgray', '#e6757b', '#CD212A', '#CD212A'], N=value_range)
background_color = np.array([1,1,1,1])
newcolors = sabah_colormap(np.linspace(0, 1, value_range))
newcolors = np.vstack((newcolors, background_color))
sabah_colormap = ListedColormap(newcolors)
hillshade = es.hillshade(sabah_topography[0], azimuth=180, altitude=1)

# Plot
plt.rcParams["figure.figsize"] = [5,5]
plt.rcParams["figure.autolayout"] = True
fig, ax = plt.subplots()
ax.imshow(sabah_topography[0], cmap=sabah_colormap)
ax.imshow(hillshade, cmap='Greys', alpha=0.3)
schools.plot(color='black', marker='x', markersize=10,ax=ax)

ax.axis('off')
plt.show()
Solution

  • As it turns out, I had given myself the hint to answering point (1), and also managed to solve (2).

    For (1), the points simply needed to be rescaled, and we get this:

    enter image description here

    I did so by getting the max/min points of the map from the underlying shapefile, and then scaling it based on the max/min points of the axes, as follows:

    # Get limit points
l = gpd.read_file('param_geo/sabah.shp')['geometry'].bounds
lat_min,lat_max,lon_min,lon_max  = l['miny'].iloc[0], l['maxy'].iloc[0], l['minx'].iloc[0], l['maxx'].iloc[0]
xmin,xmax = ax.get_xlim()
ymin,ymax = ax.get_ylim()

# Load coordinates of schools in Sabah and rescale
crs={'init':'epsg:4326'}
sch = pd.read_csv('sabah/sabah_schools.csv',usecols=['lat','lon'])
sch.lat = ymin + (sch.lat - lat_min)/(lat_max - lat_min) * (ymax - ymin)
sch.lon = xmin + (sch.lon - lon_min)/(lon_max - lon_min) * (xmax - xmin)

    For (2), the grey background is coming from the fact that the hillshade array has values outside the map area which are being mapped to grey. To remove the grey, we need to nullify these values.

    In this specific case, we can leverage on the fact that we know the top right corner of this map is "outside" the map (every country in the world will have at least one corner for which this is true, because no country is a perfect square):

    top_right = hillshade[0,-1]
hillshade[hillshade == top_right] = np.nan

    And voila, a beautiful white background:

    enter image description here

    For (3), I suspect it requires us to rescale the Polygon from the shapefile in a manner similar to how we rescaled the coordinates.