In the previous post we created an interactive map using two GeoJSON files. In this post we will show how to create those files.

First we had to read the original Excel tables. We known that pandas+xlrd can take care of that easily. The next challenge was to create a GeoJSON Polygon, for the HF-radar range radius, based on the angle and range information. The range radius looks like a wedge, so we used the matplolib's Wedge patch and exported the points as a polygon using the geojson module. Piece of cake ;-)

The last two problems we need to solve were: (1) we needed to ensure that people on Windows and Linux could run the script using the Anaconda distribution and obtain the same results; (2) use fiona to save shapefiles. The problem is that the fiona module shipped in the default channel is broken. They ship it with gdal 2.0.0 and fiona does not work with gdal 2.0 yet.

We solved both problems with one tool: conda-execute. With conda-execute we can ensure that anyone running the script will end up with the exact same environment we originally planned. And we can specify the gdal version for that environment "hard-coding" the fiona work-around in the script itself.

Take a look at the final script below. You can also see the GeoJSONs files hfradar.geojson and stations.geojson rendered on GitHub, and check the final result out in the previous blog post.

PS: Be sure to read the conda_execute blog post.

#!/usr/bin/env python

# conda execute
# env:
#  - fiona
#  - gdal <2.0.0
#  - matplotlib
#  - geojson
#  - pandas
#  - xlrd
# channels:
#  - ioos
# run_with: python

"""
Create GeoJSON Features for SECOORA Assets.

To test the features copy-and-paste the GeoJSON file onto:

http://geojson.io/

Run with,
$ conda execute -v data_frame2gis.py
to ensure same results.

"""

import json
import fiona
import pandas as pd
from matplotlib.patches import Wedge
from geojson import FeatureCollection, Feature, Polygon, Point

url = ("https://raw.githubusercontent.com/ocefpaf/"
       "secoora_assets_map/gh-pages/secoora_icons/")


status_colors = dict(Planned="orange",
                     Operational="green",
                     Permitting="yellow",
                     Construction="yellow")

platforms_icons = dict({"Fixed Surface Buoy": "buoy",
                        "Fixed Bottom Station": "circ",
                        "Fixed Bottom Mount Mooring": "tri",
                        "Fixed Coastal Station": "shore_station",
                        "HFRadar": "hfradar"})

icon = (url + "{platform}-{status}.png").format
ranges = dict({5: 190,
               8: 160,
               12: 130,
               16: 100})


def wedge(radar):
    """
    Make a HF-radar `matplotlib.patches.Wedge` from a StartAngle, SpreadAngle,
    Center (Longitude, Latitude), and Radius (radar range).

    """
    r = ranges[int(radar["MHz"])] / 111.1  # deg2km
    theta1, theta2 = radar["StartAngle"], radar["SpreadAngle"]
    center = radar["Longitude"], radar["Latitude"]
    try:
        return Wedge(center, r, theta1+theta2, theta1)
    except ValueError:
        return None


def mpl_patch2geo_polygon(patch):
    """Close the matplotlib`patch` and return a GeoJSON `Polygon`."""
    p = patch.get_path()
    vertices = p.vertices.tolist()
    return Polygon([vertices + [vertices[-2]]])


def parse_hfradar(df, **kw):
    """
    Expect a pandas.DataFrame with the following columns:
    - ResponsibleParty
    - Type
    - DisplayTitle
    - Abrreviated ID
    - Latitude
    - Longitude
    - MHz Status
    - StartAngle
    - SpreadAngle

    """
    # https://github.com/mapbox/simplestyle-spec/
    defaults = dict(stroke="#aeccae",
                    stroke_width=1,
                    stroke_opacity=0.5,
                    fill="#deffde",
                    fill_opacity=0.25)
    defaults.update(kw)

    polygons, points = [], []
    for status, group in df.groupby("Status"):
        kw = dict(platform="hfradar", status=status_colors[status])
        for name, row in group.iterrows():
            popupContent = "{} ({} MHz)".format(row["DisplayTitle"],
                                                row["MHz"])
            properties = dict(icon=icon(**kw),
                              name=name,
                              popupContent=popupContent)
            patch = wedge(row)
            if patch:
                polygon = mpl_patch2geo_polygon(patch)
            point = Point([row["Longitude"], row["Latitude"]])

            points.append(Feature(geometry=point, properties=properties))
            polygons.append(Feature(geometry=polygon, properties=defaults))
    return FeatureCollection(points+polygons)


def parse_stations(df):
    """
    Group SECOORA stations spreadsheet into PlatformType+Status.

    Expect a pandas.DataFrame with the following columns:
    - PlatformType
    - Status
    - Longitude
    - Latitude
    - LocationDescription
    - DisplayTitle (name)

    """
    features = []
    for (platformtype, status), group in df.groupby(["PlatformType",
                                                     "Status"]):
        kw = dict(status=status_colors[status],
                  platform=platforms_icons[platformtype])
        for name, row in group.iterrows():
            properties = dict(icon=icon(**kw),
                              name=row["Name"],
                              popupContent=row["LocationDescription"])
            geometry = Point([row["Longitude"], row["Latitude"]])
            feature = Feature(geometry=geometry, properties=properties)
            features.append(feature)
    return FeatureCollection(features)


def save_geojson(geojson, fname):
    """Save to GeoJSON."""
    kw = dict(sort_keys=True, indent=2, separators=(",", ": "))
    with open(fname, "w") as f:
        json.dump(geojson, f, **kw)


def save_shapefile(geojson, fname, geometry="Point"):
    """
    Save one `geometry` type from a geojson of a __geo_interface__ as a
    shapefile`.

    CAVEAT: this is a lossy conversion! I am passing along only the name
    property.

    """
    schema = {"geometry": geometry,
              "properties": {"name": "str:80"}}

    with fiona.open(fname, "w", "ESRI Shapefile", schema) as f:
        for k, feature in enumerate(geojson["features"]):
            if feature["geometry"]["type"] == geometry:
                try:
                    name = feature["properties"]["name"]
                except KeyError:
                    name = k
                f.write({
                    "geometry": feature["geometry"],
                    "properties": {"name": name},
                    })

if __name__ == "__main__":
    import os

    directory = "spreadsheets"
    save = "data"

    # Stations.
    fname = os.path.join(directory, "secoora_station_assets.xlsx")
    df = pd.read_excel(fname, index_col=2)

    geojson = parse_stations(df)
    save_geojson(geojson, fname=os.path.join(save, "stations.geojson"))
    save_shapefile(geojson, fname=os.path.join(save, "stations.shp"),
                   geometry="Point")

    # HFRadar.
    fname = os.path.join(directory, "secoora_hfradar_sites.xlsx")
    df = pd.read_excel(fname, index_col=3)

    geojson = parse_hfradar(df)
    save_geojson(geojson, fname=os.path.join(save, "hfradar.geojson"))
    save_shapefile(geojson, fname=os.path.join(save, "hfradar_point.shp"),
                   geometry="Point")
    save_shapefile(geojson, fname=os.path.join(save, "hfradar_polygon.shp"),
                   geometry="Polygon")