Skip to content

Instantly share code, notes, and snippets.

@jpwspicer
Last active February 25, 2020 22:06
Show Gist options
  • Star 2 You must be signed in to star a gist
  • Fork 1 You must be signed in to fork a gist
  • Save jpwspicer/bd738886467c56c5d029 to your computer and use it in GitHub Desktop.
Save jpwspicer/bd738886467c56c5d029 to your computer and use it in GitHub Desktop.
Examples of using Basemap with Julia

Julia Basemap Examples

Last Update: 04.16.2019 (Updated to Julia 1.0 & PyPlot 2.8)

N.B. The Basemap library will be maintained until 2020, after which it will be replaced by Cartopy. Gist on how to use Cartopy with Julia here. Announcement on Matplotlib.org here.

Contents

Basemap

Installation

To install and begin using Basemap, run the following commands in the Julia command line:

using Conda
Conda.add("basemap")

--

Python -> Julia PyPlot Translation

When using Python modules (of which Basemap is an example) in Julia rather than in Python, the syntax generally changes as follows:

  • Python: map.fillcontinents(color='coral')
  • Julia 1.0: map.fillcontinents(color="coral")
  • Julia ≤0.7: map[:fillcontinents](color="coral")

Some Python functions are not replicated in Julia and must be either substitued or calculated by hand:

Python Julia
X, Y = meshgrid(x,y) X, Y = repeat(x', length(y), 1), repeat(y, 1, length(x))
parallels = np.arange(0., 90, 10.) parallels = collect(0:10:90)
np.indices((nlats, nlons))[0,:,:] repeat([i for i = 1:nLats], 1, nLons)
np.indices((nlats, nlons))[1,:,:] repeat([j for j = 1:nLons]', nLats)

Basemap Examples

Contours
using PyPlot, PyCall
basemap = pyimport("mpl_toolkits.basemap")

# Set up orthographic map projection with perspective of satellite looking down at 45N, 100W.
# Use low resolution coastlines.
map = basemap.Basemap(projection="ortho", lat_0=45, lon_0=-100, resolution="l")

# Draw coastlines, country boundaries, fill continents.
map.drawcoastlines(linewidth=0.25)
map.drawcountries(linewidth=0.25)
map.fillcontinents(color="coral", lake_color="aqua")

# Draw the edge of the map projection region (the projection limb)
map.drawmapboundary(fill_color="aqua")

# Draw lat/lon grid lines every 30 degrees.
map.drawmeridians(collect(0:30:360))
map.drawparallels(collect(-90:30:90))

# Make up some data on a regular lat/lon grid.
nLats = 73; nLons = 145; δ = 2π/(nLons-1)
lats = repeat(0.5π .- δ*[i for i = 1:nLats], 1, nLons)
lons = repeat*[j for j = 1:nLons]', nLats)

wave = 0.75(sin.(2lats).^8).*cos.(4lons)
mean = 0.5cos.(2lats).*(sin.(2lats).^2 .+ 2)

# Compute native map projection coordinates of lat/lon grid.
x, y = map(rad2deg.(lons), rad2deg.(lats))

# Contour data over the map.
cs = map.contour(x, y, wave+mean, 15, linewidths=1.5)
title("Global Contour Lines")

Contour

--

Great Circles
using PyPlot, PyCall
# basemap = pyimport("mpl_toolkits.basemap")

lightBlue = (220/255, 220/255, 255/255)
lightGreen = (230/255, 255/255, 230/255)

# Set up mercator map projection.
map = basemap.Basemap(llcrnrlon=-100, llcrnrlat=20, urcrnrlon=20, urcrnrlat=60,
            rsphere=(6378137.00, 6356752.3142),
            resolution="l", projection="merc",
            lat_0=40, lon_0=-20, lat_ts=20)
map.drawcoastlines(linewidth=0.5)
map.fillcontinents(color=lightGreen, lake_color=lightBlue)
map.drawcountries(linewidth=0.25)
map.drawmapboundary(fill_color=lightBlue)

# Draw parallels & meridians
map.drawmeridians(collect(-180:30:180), labels=[1,1,0,1])
map.drawparallels(collect(10:20:90), labels=[1,1,0,1])

# Start, end coordinates (New York & London)
ny  = (40.78, -73.98)
lon = (51.53, 0.08)

# Draw great circle route between New York and London
map.drawgreatcircle(ny[2], ny[1], lon[2], lon[1], linewidth=3, color="m")
for (lat, lon) in [ny, lon]
    x, y = map(lon, lat)
    map.scatter(x, y, zorder=3, s=60, edgecolor="k", color="yellow")
end

title("Great Circle from New York to London")

Great Circle

--

Day / Night
using PyPlot, PyCall, Dates
basemap = pyimport("mpl_toolkits.basemap")

# Miller projection
map = basemap.Basemap(projection="mill", lon_0=180)

# Plot coastlines, draw label meridians and parallels
map.drawcoastlines(linewidth=0.5)
map.drawcountries(linewidth=0.25)
map.fillcontinents(color="coral", lake_color="aqua")
map.drawmeridians(collect(map.lonmin:60:map.lonmax+30), labels=[0,0,0,1])
map.drawparallels(collect(-90:30:90), labels=[1,0,0,0])
map.drawmapboundary(fill_color="aqua")

# Shade the night areas using current time in UTC
CS = map.nightshade(now(Dates.UTC))
title("Day/Night Map for " * Libc.strftime(time()) * " (UTC)")

Day / Night

@mafla
Copy link

mafla commented Dec 7, 2018

Hi, thanks for the examples! Any suggestion how I could do something like described here

from mpl_toolkits.basemap import Basemap, maskoceans in Julia?

@johannspies
Copy link

Thanks from my side also.

I have a problem with the following example in Julia 1.0.3:

nLats = 73; nLons = 145; delta = 2*pi/(nLons-1)
lats = convert(Array{Float64}, repmat(0.5*pi - delta*[i for i=1:nLats], 1, nLons))

MethodError: no method matching -(::Float64, ::Array{Float64,1})
Closest candidates are:
  -(::Float64, !Matched::Float64) at float.jl:397
  -(::Float64) at float.jl:387
  -(!Matched::PyObject, ::Any) at /home/js/.julia/packages/PyCall/0jMpb/src/pyoperators.jl:14
  ...

Stacktrace:
 [1] top-level scope at In[9]:3

Also

julia> Pkg.checkout("PyPlot")
ERROR: UndefVarError: checkout not defined

@jpwspicer
Copy link
Author

Apologies for the time taken to respond, GitHub does not notify authors of Gist comments.

@mafla, the Basemap and maskoceans functions can be used in Julia as follows:

using PyPlot, PyCall
basemap = pyimport("mpl_toolkits.basemap")

map = basemap.Basemap(projection="ortho", lat_0=45, lon_0=-100, resolution="l")

mdata = basemap.maskoceans(lons2, lats2, data2, resolution="c", grid=10, inlands=true)

@johannspies, the errors are due to syntax changes in Julia 1.0, to which I have updated the examples.
The first error is fixed by replacing - with .-, and the checkout command has been replaced by using Pkg; Pkg.update("PyPlot").

@ianfiske
Copy link

To fix the error

KeyError: 'PROJ_LIB'

I needed

using PyCall

py"""
import os
import conda

conda_file_dir = conda.__file__
conda_dir = conda_file_dir.split('lib')[0]
proj_lib = os.path.join(os.path.join(conda_dir, 'share'), 'proj')
os.environ["PROJ_LIB"] = proj_lib
"""

basemap = pyimport("mpl_toolkits.basemap")

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment