Figures for lectures in oceanography

In [1]:

import os
import geopandas


fname = os.path.join('data', '10m_admin_0_countries')
gdf = geopandas.read_file(fname)

south_america = gdf[gdf['CONTINENT'] == 'South America']

In [2]:

from shapely.ops import cascaded_union

geometries = [
    country['geometry'] for k, country in south_america.iterrows()
]
sa = cascaded_union(geometries)
sa.boundary

In [3]:

import geopandas

# NE data
fname = os.path.join('data', '10m_rivers_lake_centerlines')
rivers = geopandas.read_file(fname)

In [4]:

from shapely.geometry import box

sasq = box(*sa.bounds)

mask = [
    sasq.contains(river['geometry']) for k, river in rivers.iterrows()
]

sa_rivers = rivers[mask]

# I love ? unicode!
set(sa_rivers['name'])

{'Ajaj?',
 'Amapari',
 'Amazonas',
 'Apure',
 'Araguaia',
 'Aripuan?',
 'Atrato',
 'B?o-B?o',
 'Beni',
 'Bermejo',
 'Bois',
 'Braco Menor',
 'Branco',
 'Caquet?',
 'Carcara??',
 'Caron?',
 'Casiquiare',
 'Cauca',
 'Chico',
 'Chirrip del Atlntico',
 'Chubut',
 'Coari',
 'Coco',
 'Cojedes',
 'Colorado',
 'Contas',
 'Copiap?',
 'Corantijn',
 'Corumba',
 'Courantyne',
 'Cuiaba',
 'Culuene',
 'Cuyuni',
 'Desaguadero',
 'Deseado',
 'Doce',
 'Dulce',
 'Essequibo',
 'Gallegos',
 'Grande',
 'Gu?rico',
 'Guain?a',
 'Guapor?',
 'Gurupi',
 'Huallaga',
 'Ibicu?',
 'Indai?',
 'Iriri',
 'Itiquira',
 'Ituxi',
 'Iva?',
 'Ivinheima',
 'Jamanxim',
 'Japur?',
 'Jari',
 'Javari',
 'Jejui Guazu',
 'Jequitinhonha',
 'Jiparan?',
 'Juru?',
 'Juruena',
 'Juta?',
 'Kabalebo',
 'Lempa',
 'Limay',
 'Litani',
 'Madeira',
 'Madre de Dios',
 'Magdalena',
 'Maicuru',
 'Mamor?',
 'Mara??n',
 'Maroni',
 'Mearim',
 'Mendoza',
 'Meta',
 'Mico',
 'Mira',
 'Miranda',
 'Napo',
 'Negro',
 'Neuqu?n',
 'Nhamund?',
 None,
 'Oco?a',
 'Orinoco',
 'Panama Canal',
 'Paraiba',
 'Paran?',
 'Parana?ba',
 'Paranapanema',
 'Pardo',
 'Parna?ba',
 'Piau?',
 'Pilcomayo',
 'Pisco',
 'Pur?s',
 'Putumayo',
 'Ribeira',
 'Rio Rapel',
 'Rio das Mortes',
 'Ro Grande de Matagalpa',
 'S?o  Francisco',
 'Salado',
 'San Juan',
 'San Martin',
 'San Miguel',
 'Santa',
 'Santa Cruz',
 'Santa Mara',
 'Santiago',
 'Suriname',
 'Tambo',
 'Tapaj?s',
 'Taquari',
 'Teles Pires',
 'Tercero',
 'Tiet?',
 'Tocantins',
 'Trombetas',
 'Tumbes',
 'Uaup?s',
 'Ucayali',
 'Una',
 'Unini',
 'Uruguay',
 'Velhas',
 'Ventuari',
 'Verde',
 'Xingu',
 'Yapacani'}

In [5]:

include = {
    2: 'Amazonas',
    78: 'Negro',
    91: 'Madeira',
    99: 'Paraná',
    100: 'Araguaia',
    104: 'São Francisco',
    108: 'Paranaíba',
    115: 'Tocantins',
    197: 'Xingu',
    200: 'Tapajós',
    284: 'Jequitinhonha',
    293: 'Uruguay',
    335: 'Tocantins',
    400: 'Paranapanema',
    436: 'Parnaíba',
    515: 'Doce',
    688: 'Tietê',
}

In [6]:

import geopandas


# http://sigaceivap.org.br/map
fname = os.path.join('data', 'cam_rioparaibasul')
rps = geopandas.read_file(fname)

# https://gis.stackexchange.com/questions/98371/how-to-get-shapefile-of-river-from-openstreetmap
fname = os.path.join('data', 'paraguay_river.geojeon')
pry = geopandas.read_file(fname)

In [7]:

%matplotlib inline

import cartopy
import cartopy.crs as ccrs
from cartopy.feature import NaturalEarthFeature

import matplotlib.pyplot as plt


def make_map(projection=ccrs.PlateCarree(), figsize=(9, 9)):
    fig, ax = plt.subplots(
        figsize=figsize,
        subplot_kw=dict(projection=projection)
    )

    states = NaturalEarthFeature(
        category='cultural',
        scale='50m',
        facecolor='none',
        name='admin_1_states_provinces_shp'
    )

    ax.add_feature(cartopy.feature.LAND, alpha=0.25)
    ax.add_feature(states, edgecolor='gray', alpha=0.4)
    ax.add_feature(cartopy.feature.BORDERS, linestyle=':', alpha=0.2)
    return fig, ax

In [8]:

fig, ax = make_map()

bbox = (-76.9043, -33.7500, -36.0313, 6.6646)
ax.set_extent(bbox)

for rivernum, name in include.items():
    river = rivers[rivers['rivernum'] == rivernum]
    ax.add_geometries(
        [geom for geom in river['geometry']],
        ccrs.PlateCarree(),
        edgecolor='blue',
        facecolor='none',
        alpha=0.75,
    )

ax.add_geometries(
    [geom for geom in rps['geometry']],
    ccrs.PlateCarree(),
    edgecolor='blue',
    facecolor='none',
    alpha=0.75,
)


ax.add_geometries(
    [geom for geom in pry['geometry']],
    ccrs.PlateCarree(),
    edgecolor='blue',
    facecolor='none',
    alpha=0.75,
)

<cartopy.mpl.feature_artist.FeatureArtist at 0x7f3fc90ce3c8>

png

In [9]:

import folium
from folium.plugins import Fullscreen

lon = -(76.9043 + 33.7500) / 2.
lat = (-36.0313+6.6646) / 2.


def style_function(feature):
    return {
        'fillColor': '#ffaf00',
        'color': 'blue',
        'weight': 1.5,
    }


def highlight_function(feature):
    return {
        'fillColor': '#ffaf00',
        'color': 'blue',
        'weight': 4,
    }


m = folium.Map(location=[lat, lon], zoom_start=4)
Fullscreen(position='topright', force_separate_button=True).add_to(m)

for rivernum, name in include.items():
    river = rivers[rivers['rivernum'] == rivernum]
    line = folium.GeoJson(
        data=river['geometry'].__geo_interface__,
        style_function=style_function,
        highlight_function=highlight_function,
    )
    folium.Popup(name).add_to(line)
    line.add_to(m)

line = folium.GeoJson(
    data=rps['geometry'].__geo_interface__,
    style_function=style_function,
    highlight_function=highlight_function,
)

folium.Popup('Paraíba Sul').add_to(line)
line.add_to(m)

line = folium.GeoJson(
    data=pry['geometry'].__geo_interface__,
    style_function=style_function,
    highlight_function=highlight_function,
)

folium.Popup('Paraguay').add_to(line)
line.add_to(m)

m