Missing support for lambert azimuthal projection

[uandrae]

Is your feature request related to a problem? Please describe.

In the DE330 DestinE project we output hydrological data in the lambert azimuthal projection. When trying to plot this with earthkit.qickmap it fails as shown below. Is this projection not yet supported?

[snh@ac6-101 ~]$ python3 ./ek.py 
Traceback (most recent call last):
  File "/etc/ecmwf/nfs/dh1_home_b/snh/./ek.py", line 8, in <module>
    qmap.block(data)
  File "/home/snh/.local/lib/python3.10/site-packages/earthkit/plots/quickmap.py", line 9, in wrapper
    getattr(subplot, function.__name__)(*args, **kwargs)
  File "/home/snh/.local/lib/python3.10/site-packages/earthkit/plots/components/subplots.py", line 225, in wrapper
    return self._extract_plottables(
  File "/home/snh/.local/lib/python3.10/site-packages/earthkit/plots/components/subplots.py", line 371, in _extract_plottables
    x_values = source.x_values
  File "/usr/local/apps/python3/3.10.10-01/lib/python3.10/functools.py", line 981, in __get__
    val = self.func(instance)
  File "/home/snh/.local/lib/python3.10/site-packages/earthkit/plots/sources/earthkit.py", line 178, in x_values
    return self.extract_x()
  File "/home/snh/.local/lib/python3.10/site-packages/earthkit/plots/sources/earthkit.py", line 143, in extract_x
    return self.extract_xy()[0]
  File "/home/snh/.local/lib/python3.10/site-packages/earthkit/plots/sources/earthkit.py", line 136, in extract_xy
    points = self.data.to_points(flatten=False)
  File "/home/snh/.local/lib/python3.10/site-packages/earthkit/data/core/fieldlist.py", line 386, in to_points
    raise ValueError("to_points(): geographical coordinates in original CRS are not available")
ValueError: to_points(): geographical coordinates in original CRS are not available

Describe the solution you'd like

No response

Describe alternatives you've considered

No response

Additional context

No response

Organisation

DE330@DestinE

[JamesVarndell]

Hi @uandrae, thanks for opening this issue.

Could you please share the code that you ran which triggered this error, and if possible the data that you were trying to plot? The Lambert azimuthal equal area projection is supported in earthkit-plots, but in this case it looks like earthkit is unable to extract that original Lambert gridpoints from your data, so it would be very useful to see the original file you're trying to plot. Many thanks!

[uandrae]

Thanks for getting back to me. Data is extracted from MARS as

retrieve,
  class=ce,
  date=2023-08-10,
  expver=1,
  levtype=sfc,
  model=lisflood,
  origin=ecmf,
  param=240023,
  step=6,
  stream=efas,
  time=00:00:00,
  type=fc,
  target="efas_2023 .grib"

using the earthkit plotting

import earthkit.data
import earthkit.maps
infile = "efas_2023.grib"
data = earthkit.data.from_source("file", infile)
earthkit.maps.quickplot(data)

Fails as above. Whereas if I use a date = 2024-08-10 where the data is stored in latlon in MARS the plotting works.

[uandrae]

Did you manage to reproduce my error?

[JamesVarndell]

Hi @uandrae,

Many thanks for sharing your code - I have managed to reproduce the error, and I'm working on a fix which should be released early this week. Thanks for your patience! :)

[JamesVarndell]

Hi @uandrae,

I've released a fix for this at earthkit-plots version 0.2.3 - you should be able to get it with:

pip install earthkit-plots --upgrade

Note also that this issue was migrated from the earthkit-maps project to earthkit-plots; these used to be two separate libraries but for many good reasons (and while earthkit is still in beta) we decided to merge earthkit-maps into earthkit-plots. Apologies for the inconvenience!

If you want to run the equivalent of earthkit.maps.quickplot(data) in earthkit-plots, you'll need to do the following:

import earthkit.plots.quickmap as qmap

qmap.plot(data)

Give that a go and let me know if it works as expected. Many thanks!

[uandrae]

I'm afraid that

import earthkit.data
import earthkit.plots.quickmap as qmap

infile = "/scratch/snh/txt2nc/efas_2024.grib"
data = earthkit.data.from_source("file", infile)
qmap.plot(data)

just ends without errors, or any generated plots.

[uandrae]

Tried again with

[snh@ac6-102 ~]$ module list

Currently Loaded Modules:
  1) gcc/8.5.0   2) prgenv/gnu   3) python3/3.11.10-01

 

[snh@ac6-102 ~]$ pip list installed | grep earthkit
earthkit                      0.8.0
earthkit-data                 0.10.4
earthkit-geo                  0.2.0
earthkit-meteo                0.1.1
earthkit-plots                0.2.3
earthkit-plots-default-styles 0.0.1
earthkit-regrid               0.3.4
earthkit-transforms           0.3.3

Running

import earthkit.data
import earthkit.plots.quickmap as qmap

infile = "/scratch/snh/txt2nc/efas_2023.grib"
data = earthkit.data.from_source("file", infile)
qmap.plot(data)

Still gives the same error:

[snh@ac6-102 ~]$ python3 ./ek.py 
Traceback (most recent call last):
  File "/etc/ecmwf/nfs/dh1_home_b/snh/./ek.py", line 6, in <module>
    qmap.plot(data)
  File "/usr/local/apps/python3/3.11.10-01/lib/python3.11/site-packages/earthkit/plots/quickmap.py", line 9, in wrapper
    getattr(subplot, function.__name__)(*args, **kwargs)
  File "/usr/local/apps/python3/3.11.10-01/lib/python3.11/site-packages/earthkit/plots/components/subplots.py", line 596, in plot
    return method(*args, style=style, **kwargs)
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/usr/local/apps/python3/3.11.10-01/lib/python3.11/site-packages/earthkit/plots/components/subplots.py", line 227, in wrapper
    return self._extract_plottables(
           ^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/usr/local/apps/python3/3.11.10-01/lib/python3.11/site-packages/earthkit/plots/components/subplots.py", line 374, in _extract_plottables
    x_values = source.x_values
               ^^^^^^^^^^^^^^^
  File "/usr/local/apps/python3/3.11.10-01/lib/python3.11/functools.py", line 1001, in __get__
    val = self.func(instance)
          ^^^^^^^^^^^^^^^^^^^
  File "/usr/local/apps/python3/3.11.10-01/lib/python3.11/site-packages/earthkit/plots/sources/earthkit.py", line 178, in x_values
    return self.extract_x()
           ^^^^^^^^^^^^^^^^
  File "/usr/local/apps/python3/3.11.10-01/lib/python3.11/site-packages/earthkit/plots/sources/earthkit.py", line 143, in extract_x
    return self.extract_xy()[0]
           ^^^^^^^^^^^^^^^^^
  File "/usr/local/apps/python3/3.11.10-01/lib/python3.11/site-packages/earthkit/plots/sources/earthkit.py", line 136, in extract_xy
    points = self.data.to_points(flatten=False)
             ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/usr/local/apps/python3/3.11.10-01/lib/python3.11/site-packages/earthkit/data/core/fieldlist.py", line 386, in to_points
    raise ValueError("to_points(): geographical coordinates in original CRS are not available")
ValueError: to_points(): geographical coordinates in original CRS are not available

[JamesVarndell]

Hi @uandrae, it looks like you're not actually using earthkit-plots 0.2.3, even though your module list suggests that you are - the telltale sign is that your error references line 136 in earthkit/plots/sources/earthkit.py, which differs to line 136 in 0.2.3. Could you be using a virtual environment or somthing similar that differs to your loaded modules?

[uandrae]

OK, this is for python3/3.11.10-01, for python3/3.10.10-01 I expect a plot to appear, but nothing happens?

[JamesVarndell]

Hi @uandrae, the earthkit.plots.quickmap.plot() function generates a new Map object, but does not save or show this object without further steps (unless you're in a Jupyter notebook, in which case the plot should appear as the output of the cell in which qmap.plot() was called).

If you're running this in a script, you need to call .figure.save() on the result of qmap.plot():

qmap.plot(data).figure.save("my_image.png")

Side note: if you upgrade to earthkit-plots 0.2.4, you can skip the .figure part, as there is a new convenience save method on the Map object returned from qmap.plot():

qmap.plot(data).save("my_image.png")