Unrotate latitude longitude coordinates plot mismatch

[Peter Willetts]

Hi

I've tried to make something close to a minimum working example, to illustrate the problem I am having (but needs a pp cube of data).

The top figure shows a plot of mean rainfall, as plotted from the original rotated grid latitudes and longitudes, and the second plot is after the coordinates have been unrotated and the cube updated with them, as in the function in the script.

I can't work out if this is something to do with my method of unrotating and updating the cube coordinates, or with the plotting.

Thanks in advance

Peter

from iris.coords import DimCoord
from iris.coord_systems import RotatedGeogCS, GeogCS
from iris.analysis.cartography import unrotate_pole

from iris.exceptions import CoordinateNotFoundError
from numpy import meshgrid, array

import pdb

def unrotate_pole_update_cube(cube):
    try:

        lat = cube.coord('grid_latitude').points
        lon = cube.coord('grid_longitude').points

        lat_coord_name='grid_latitude'
        lon_coord_name='grid_longitude'

    except CoordinateNotFoundError:

        print 'Uses latitude/longitude coords, not grid_latitude/grid_longitude'
        pass

    try:
        cs = cube.coord_system()
    except TypeError:
        cs=[]
 
    if isinstance(cs, RotatedGeogCS):
       
        lons, lats = meshgrid(lon, lat)
       
        lons,lats = unrotate_pole(lons,lats, cs.grid_north_pole_longitude, cs.grid_north_pole_latitude)

        lon=lons[0]
        lat=lats[:,0]
               
        #pdb.set_trace()

        for i, coord in enumerate (cube.coords()):
            if coord.standard_name==lat_coord_name:
                lat_dim_coord_cube = i
            if coord.standard_name==lon_coord_name:
                lon_dim_coord_cube = i
             
        # IRIS unrotate_pole uses the default cartopy.crs.Geodetic, which is
        # WGS84 by default
        #wgs84_cs = GeogCS(semi_major_axis= 6378137.,
          #                               inverse_flattening=298.257223563)

        wgs84_cs = GeogCS(6371229.0)
       
        cube.remove_coord(lat_coord_name)
        cube.remove_coord(lon_coord_name)

        cube.add_dim_coord(DimCoord(points=lat, standard_name='latitude', units='degrees', coord_system=wgs84_cs), lat_dim_coord_cube)
        cube.add_dim_coord(DimCoord(points=lon, standard_name='longitude', units='degrees', coord_system=wgs84_cs), lon_dim_coord_cube)

    return cube

import iris
import cartopy.crs as ccrs
import matplotlib.pyplot as plt

import numpy as np

import iris.plot as iplt
diag = 'rain_mean'

pp_file_path='/projects/cascade/pwille/moose_retrievals/'

experiment_id = 'djzny'
expmin1 = experiment_id[:-1]

lon_high = 101
lon_low = 64.
lat_high = 33.
lat_low =-6.

min_contour = 0.
max_contour = 3.
tick_interval=0.2

clevs = np.linspace(min_contour, max_contour,32)

fu = '/projects/cascade/pwille/moose_retrievals/%s/%s/%s.pp' % (expmin1, experiment_id, diag)

cube  = iris.load_cube(fu)*3600

# Mean Rain

fig = plt.figure()
ax = plt.axes(projection=ccrs.PlateCarree(), extent=(lon_low,lon_high,lat_low,lat_high))

pdb.set_trace()

iplt.contourf(cube, clevs)
ax.coastlines(resolution='110m', color='k', linewidths=1.5)

plt.show()

# Unrotated Mean Rain

rain_mean = unrotate_pole_update_cube(cube)

fig = plt.figure()
ax = plt.axes(projection=ccrs.PlateCarree(), extent=(lon_low,lon_high,lat_low,lat_high))

#iplt.contourf(rain_mean, clevs, cmap=cmap)
plt.contourf(rain_mean.coord('longitude').points, rain_mean.coord('latitude').points, rain_mean.data, clevs)
ax.coastlines(resolution='110m', color='k', linewidths=1.5)

plt.show()
#plt.savefig()

[Peter Willetts]

Should have said! Rainfall location doesn't match up in plots.

[Andrew Dawson]

Your method is wrong, the basic problem is with this part:

lons, lats = meshgrid(lon, lat)

lons, lats = unrotate_pole(lons,lats, cs.grid_north_pole_longitude, cs.grid_north_pole_latitude)
lon=lons[0]    # DON'T DO THIS
lat=lats[:,0]  # OR THIS

The rotated latitudes and longitudes (coordinates grid_latitude and grid_longitude) are 1D and orthogonal, meaning you can describe a 2d space with the 2 1-d coordinates. The unrotated coordinates are 2D, meaning that each row of longitude is different, and likewise each column of latitude is different, which means you cannot choose just one slice of each to represent the full coordinate space, you need all of them.

You can create 2D AuxCoord instances instead of DimCoords and add these to the cube. The iris plotting functions will work just fine with 2D AuxCoords so there should be no need to give them and special treatment. Not you also do no need to remove the original grid_latitude and grid_longitude dimension coordinates as they are still valid, instead you can just specify which coordinates to plot against.

However, I'm not sure why you want to do this at all... Cartopy/Iris does this for you when it plots rotated pole data in a Plate Carree projection (and indeed any other CRS and projection combination). If you really want the data on a Plate Carree projection with 1D dimension coordinates (lat/lon) then you will need to regrid/reproject to move the actual data points onto a regular grid in lat/lon space.

[Phil Elson]

As ever, awesome answer Andrew!

To second what was already stated, if the reason you are doing this is purely for visualisation, cartopy will deal with your data just fine. Though given you have already written some pretty in-depth cartopy code above, I assume you already know this, so it would be really interesting to hear your use case and we can see if there is already functionality to get to your desired destination.

Cheers,

--
You received this message because you are subscribed to the Google Groups "Iris" group.
To unsubscribe from this group and stop receiving emails from it, send an email to scitools-iri...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[Andrew Dawson]

Peter responded offline:

OK, thanks
I need to unrotate onto a lat lon grid to do an area weighted regrid.