Simple HiSPARC station plot with Basemap

[twhyntie]

Not quite sure where to put this but here's some simple code for getting the station lats and lons and plotting them on a map (provided by Basemap):

# Do:
# $ conda install basemap
# to get this - ideally create a new Conda environment for it.
from mpl_toolkits.basemap import Basemap

#...for the MATH.
import numpy as np

#...for the plotting.
import matplotlib.pyplot as plt

# Get the HiSPARC station information.
from sapphire import HiSPARCStations

def get_latlontext(cluster):
    """Create list of latitude, longitudes, and legend text for a cluster

    Make a list of locations for each of station and its detectors, for the
    stations the number is included.

    :param cluster: HiSPARCStations object.

    """

    lats = []
    lons = []
    for station in cluster.stations:
        latitude, longitude, _ = station.get_lla_coordinates()
        lats.append(latitude)
        lons.append(longitude)
        for detector in station.detectors:
            latitude, longitude, _ = detector.get_lla_coordinates()
            lats.append(latitude)
            lons.append(longitude)
    return lats, lons

## List of the latitudes.
lats = []

## List of the longitudes.
lons = []

# Loop over thae stations (trial and error - FIXME!!!).
for i in range(100,200):
    #lats, lons = get_latlontext(HiSPARCStations([102], force_stale=True))
    try:
        mylats, mylons = get_latlontext(HiSPARCStations([i], force_stale=True))
    except KeyError:
        continue # skip if station number not found...

    # FIXME: just awful, I apologise.
    for j, lat in enumerate(mylats):
        lats.append(lat)
        lons.append(mylons[j])

# Create a new figure.
fig=plt.figure()

# Create the axes.
ax=fig.add_axes([0.1,0.1,0.8,0.8])

# Setup a mercator map projection.
m = Basemap(llcrnrlon=-20.,llcrnrlat=45.,urcrnrlon=20.,urcrnrlat=60.,\
            rsphere=(6378137.00,6356752.3142),\
            resolution='l',projection='merc',\
            lat_0=40.,lon_0=-20.,lat_ts=20.)

# Draw the coastlines.
m.drawcoastlines()

# Draw the continents.
m.fillcontinents()

# Draw the parallels.
m.drawparallels(np.arange(10,90,20),labels=[1,1,0,1])

# Draw the meridians.
m.drawmeridians(np.arange(-180,180,30),labels=[1,1,0,1])

# Convert to x and y coordinates suitable for plotting.
x,y = m(lons, lats)

# Plot the points.
m.plot(x, y, 'bo', markersize=5)

# Set a title.
ax.set_title('HiSPARC Stations')

# Show the plot!
plt.show()

# Save the figure.
fig.savefig("map.png")

It produces this:

Not amazing but at least I've got something out with my CernVM Sapphire. Thanks again for the Conda tip - it's brilliant!

[153957]

Some hints:

Get station numbers (docs)

from sapphire import Network
Network().station_numbers()

Get cluster object with all stations (docs)

from sapphire import HiSPARCNetwork
HiSPARCNetwork(force_stale=True)

You should end up with all stations that can be seen here: http://data.hisparc.nl/maps/

[twhyntie]

Thanks! :smile: