Simply cartopy example - Global plot but data not present.

[Leighton Regayre]

Hi all,

I'm attempting to follow a simple Cartopy example and am running into these problems:

Firstly I can generate a Robinson plot with coastlines, but cannot plot the data onto it.

Secondly, if I plot the data without transformation, it is rotated by 90 degrees and replicated.

The code I'm using is:
import numpy as np
from math import *
import os
import matplotlib.pyplot as plt
from netCDF4 import Dataset as netcdf_dataset
import iris
import iris.plot as iplt
import cartopy.crs as ccrs
from mpl_toolkits.basemap import Basemap
import pylab as plb


SO2data=np.loadtxt('data.dat', skiprows=0, dtype='float32')

lat = SO2data[:,0]
long = SO2data[:,1]
so2_rt = SO2data[:,2]

datalen = so2_rt.size

# Check for index of the first value in lat.
long_first = long[0]
index_long = np.argwhere(long==long_first)

# Here is where the 1st repeat of the original longitudinal value occurs:
long_len = index_long[1]
# This is the number of different longitudinal values for each lat. (180)
# The emissions array will be this wide.

# Take a slice of the longitude.
long_universal = long[:long_len]

lat_universal = lat[0:datalen:long_len]

# Now emissions can be converted into 2D:
so2_rt_2d = so2_rt.reshape(np.asscalar(long_len), lat_universal.size)

log_so2 = np.log(so2_rt_2d)

# Here's the construction of the grid for contour:
X,Y = plb.meshgrid(lat_universal,long_universal)

# so some attempts to plot:

ax = plt.axes(projection=ccrs.Robinson())

plt.contour(X, Y, log_so2, transform=Robinson())
ax.coastlines()

plt.show()


#Thanks

[Andrew Dawson]

Hi Leighton,

There are actually a few different problems with the script you attached, the most fundamental is that you mixed up latitude and longitude in your input file, which is why your data is rotated and doubled. As for the plotting, the transform keyword to a cartopy plotting function specifies the source transform, rather than the destination. The data you attached is on a regular lat-lon grid so the correct transform is ccrs.PlateCarree(), this tells cartopy it has to transform data from PlateCarree to Robinson (or something essentially like that, my explanation might not quite hold up to rigorous scrutiny from the experts!).

I've pasted in a simplified script and attached the output so you can tell if it does what you want. The script assumes knowledge of the order of the columns in the input file and assumes the input data is on a regular grid. The way I have written it may not be the best or most practical way, but it should give you the right idea:

import numpy as np
import matplotlib.pyplot as plt
import cartopy.crs as ccrs


# Read data from a text file. Data is in three columns: lon, lat, SO2.

SO2data=np.loadtxt('data.dat', skiprows=0, dtype='float32')

longitude_all = SO2data[:, 0]
latitude_all = SO2data[:, 1]
so2_rt = SO2data[:,2]

# Find the number of latitudes and longitudes assuming a regular grid.
datalen = so2_rt.shape[0]
nlat = len(np.where(longitude_all == longitude_all[0])[0])
nlon = datalen // nlat

# Extract 1D latitude and longitude coordinates.
latitude = latitude_all[:nlat]
longitude = longitude_all[::nlat]

# Re-shape SO2 to 2D and lat-lon order.
so2_rt_2d = so2_rt.reshape([nlon, nlat]).transpose()
log_so2 = np.log(so2_rt_2d)

# Plot log(SO2) on a Robinson projection.
ax = plt.axes(projection=ccrs.Robinson())
ax.contour(longitude, latitude, log_so2, transform=ccrs.PlateCarree())
ax.coastlines()
plt.show()

[Leighton Regayre]

Hi Andrew,

Thanks for taking the time to rewrite the code. I was surprised by the usage of the transform keyword. I haven't seen any examples that indicate how to transform a regularly gridded dataset onto an alternative axis, so that's really useful.

Your code is a lot snappier than my original as well, which helps.

Thanks again,

Leighton.

[Phil Elson]

Awesome answer Andrew!

I'm glad it helped Leighton - the galley still needs a lot more examples, but the use of transform exists in the following example: http://scitools.org.uk/cartopy/docs/latest/examples/waves.html and indeed in the documentation here: http://scitools.org.uk/cartopy/docs/latest/matplotlib/advanced_plotting.html#contour-plots

Cheers,

[bblay]

Awesome answer Andrew!

yep

my explanation might not quite hold up to rigorous scrutiny from the experts!

other experts ;)