2

I would like to plot data on a grid (which is in LCC projection) with Cartopy, so that the data fills the entire axis (and also axes, but that is not the issue here).

To make it clearer, here is what I do with Cartopy:

import cartopy.crs as ccrs
import numpy as np
import pyproj as p4
from mpl_toolkits.basemap import Basemap

lalo = #read latitudes and longitudes of my grid defined in a special LCC projection (below)

lat = np.reshape(lalo[:,1],(ny,nx))
lon = np.reshape(lalo[:,0],(ny,nx))
minlat = lat[0,0]
maxlat = lat[-1,-1]
minlon = lon[0,0]
maxlon = lon[-1,-1]
Z = np.ones((ny,nx)) #some data

#grid definition for cartopy:
myproj = ccrs.LambertConformal(central_longitude=13.3333, central_latitude=47.5,
                               false_easting=400000, false_northing=400000,
                               secant_latitudes=(46, 49))
fig = plt.figure()
ax = plt.axes(projection = myproj)    
plt.contourf(lon, lat, Z)#, transform=myproj) 
#no difference with transform option as lon,lat are already in myproj projection

The result is an image which does not fill the entire axis, but looks like this: Cartopy

When using Basemap like this:

a=6377397.155
rf=299.1528128
b= a*(1 - 1/rf)
m = Basemap(projection='lcc', resolution='h', rsphere=(a,b),
        llcrnrlon=minlon,llcrnrlat=minlat,urcrnrlon=maxlon,urcrnrlat=maxlat,
        llcrnrx=400000, llcrnry=400000,
        lat_1=46, lat_2=49, lat_0=47.5, lon_0=13.3333, ax=ax)
x,y = m(lon,lat)
m.contourf(x,y,Z)

I get the following (desired) image: Basemap

And finally, when using proj4 to convert lon and lat using this definition p4.Proj('+proj=lcc +lat_1=46N +lat_2=49N +lat_0=47.5N +lon_0=13.3333 +ellps=bessel +x_0=400000 +y_0=400000') I again get the desired image:

proj4

Is there any possibility to achieve this in cartopy as well?

In other words, I would like to have a plot where the data shows up in a perfect rectangle, and the background map is distorted accordingly, i.e. something like the opposite of this example (cannot install iris package, otherwise I would have tried with this example)

I have tried a few things like:

  • building a custom class for my projection as done here, just to be sure that the parameters are all set correctly (as in my proj4 definition).
  • played around with aspect ratios, but they only affect the axes not the axis, and a few more things.

Any help is greatly appreciated!

Community
  • 1
  • 1
yngwaz
  • 769
  • 1
  • 9
  • 17

1 Answers1

2

The important piece of information that is missing here is that your data is in lats and lons, not in the Cartesian transverse Mercator coordinate system. As a result you will need to use a Cartesian coordinate system which speaks lats and lons (spherical contouring has not been implemented at this point). Such a coordinate system exists in the form of the PlateCarree crs - so simply passing this as the transform of the contoured data should put your data in the right place.

plt.contourf(lon, lat, Z, transform=ccrs.PlateCarree())

This really highlights the fact that the default coordinate system of your data, is the same as that of the map, which in most cases is not longitudes and latitudes - the only way to change the CRS of your data is by passing the transform keyword.

HTH

pelson
  • 21,252
  • 4
  • 92
  • 99
  • Thank you, this works perfectly. My misconception was that I thought I had created a map in the same coordinate system as my data (ax = plt.axes(projection = myproj)). – yngwaz Oct 16 '14 at 08:55
  • Yep - that is the restriction that a tool like Basemap imposes. Cartopy allows you to completely separate the map your drawing from the coordinate system of your data, and it should all just work. – pelson Oct 16 '14 at 09:09