General questions about using Satpy

[gersmit]

Dear Mr. Hoese, hello DJ,

My name is Ger from the Netherlands and I follow your tutorial-satpy-half-day course via Github. All okay but a small problem question: When I run HRV via the scn_cropped I get the message "ValueError: Can't compare areas of different types". Is it possible to explain to me how to deal with HRV images.

Many thanks in advance for an answer

Regards Ger

[djhoese]

Hi @gersmit. Are you loading anything else besides HRV? Do you have the full traceback for the error message? This typically comes up when there are more than one type of data in the Scene and the Scene can't determine which dataset is highest or lowest resolution. I can provide more detail when I know exactly which part of the code is failing (via the traceback you will provide).

[gersmit]

DJ, thanks for the fast replay. All images works fine with scn_cropped, only HRV not. Script is not to match and i think the problem is the scn_cropped, but how do i work with HRV because of the resolution? This is the test i do with all the datasets : dataset_names = scn.all_dataset_names() print('\n'.join(map(str, dataset_names))) scn.load(['HRV']) scn_cropped = scn.crop(ll_bbox=(-10, 34, 25, 54)) -> (europe) local_scene = scn_cropped local_scn.show('HSV')

ERROR = Traceback (most recent call last): File "/home/pi/eumetview/test_HRV.py", line 83, in scn_cropped = scn.crop(ll_bbox=(-10, 34, 25, 54)) File "/home/pi/.local/lib/python3.7/site-packages/satpy/scene.py", line 604, in crop min_area, area, ll_bbox, xy_bbox) File "/home/pi/.local/lib/python3.7/site-packages/satpy/scene.py", line 486, in _slice_area_from_bbox x_slice, y_slice = src_area.get_area_slices(dst_area) File "/usr/local/lib/python3.7/dist-packages/pyresample-1.15.0-py3.7-linux-armv7l.egg/pyresample/geometry.py", line 428, in get_area_slices raise NotImplementedError NotImplementedError

[djhoese]

Just after the .load call, what does print(scn['HRV'].attrs['area']) show you? Is there a chance you aren't providing all of the SEVIRI segment files to the Scene object when you create it?

[gersmit]

This is wat i get after lines files = {"seviri_l1b_native" : file} scn = Scene(filenames=files)

print load = : <pyresample.geometry.StackedAreaDefinition object at 0x937066b0> print('\n'.join(map(str, dataset_names))) give all dataset names including HSV

[djhoese]

This would suggest you aren't providing all of the segment files for this data. Or that you are providing too many (more than one time step). Could you show me the list of filenames you are providing to the Scene?

[gersmit]

I think you mean these names? HRV IR_016 IR_039 IR_087 IR_097 IR_108 IR_120 IR_134 VIS006 VIS008 WV_062 WV_073

[djhoese]

No, sorry. I mean the files in:

files = {"seviri_l1b_native" : file}
scn = Scene(filenames=files)

[gersmit]

files : {'seviri_l1b_native': ['MSG4-SEVI-MSG15-0100-NA-20201022155743.481000000Z-NA.nat']}

[djhoese]

Sorry @gersmit, I'm not very familiar with seviri data. I've asked some of the other pytroll developers about this (the ones from Europe) and they reminded me that HRV channel is a special case. It is made up of two separate image portions. If you used the HRIT reader it has the functionality to automatically join these two portions together in a single large image, but the native reader you are using doesn't have that functionality yet.

I'd recommend using the HRIT data if at all possible since it will behave much better. If that isn't an option then you'll unfortunately need to resample the HRV data before trying to crop it. You can try to do new_scn = scn.resample('seviri_0deg') which will resample the data using nearest neighbor to a 3km version of the seviri full disk sector. Depending on your use case you may want to resample to an area that is 1km resolution (the native HRV resolution). Unfortunately we don't have a builtin AreaDefinition made for this so you'd have to make your own. From my understanding, HRIT is the typical format used by Satpy's users for this reason.

Let me know how else I can help. Sorry this probably wasn't the answer you wanted.

[gersmit]

First of all, I would like to thank you very much for taking some time out to help me. The technique can cause problems especially if you don't use it every day.

Do you have examples for me how I can deal with creating an AreaDefinition or how I can possibly make it work or tackle the resolution. The technique you use has only been known to me for a month because Eumetsat has switched from MSG files via FTP to data via URL. So I know little about it except for what I learned from your course.

I hope you have some examples so I can experiment a bit and hope I can ask you again later.

Regards Ger PS are there other web pages where I can learn more about this topic,

[djhoese]

I will ask some of the other satpy devs to take a look at this and see if they can help. As for other web pages:

The pyresample documentation has everything about AreaDefinitions (see Geometry definitions and Geometry Utilities on the left): https://pyresample.readthedocs.io/en/latest/

The Satpy documentation: https://satpy.readthedocs.io/en/stable/

Various satpy notebooks that may help: https://github.com/pytroll/pytroll-examples/tree/master/satpy

A tutorial that was taught at EUMETSAT: https://github.com/pytroll/pytroll-examples/tree/master/pytroll-exercises-two-day

Otherwise, to help you further, I'd like to know what your end goal is exactly? Are you trying to make images? What will you use them for? Do you need to crop the Scene? There may be other ways to do what you are trying to do, we just need to know what that is.

[gersmit]

My goal is a hobby, I am retired and I love computers and technology to see the world from above. I already have heat maps and use Python to create 15-minute moving images of Europe. But I want more, every time I am busy and read articles like yours I want more, more beautiful and better.

I now say good night and keep in touch and wish you a pleasant working day. Talk to you you later

[pnuu]

The "correct" way would be to resample the data to an area definition. To create the exact are you'd get with the .crop() method you can take the width, height and extent values from the cropped "normal" channel area, and modify the existing seviri_0deg area @djhoese mentioned and save it (with a different name) to your own $PPP_CONFIG_DIR/areas.yaml file.

Please check the values with the native format data, I have a feeling the data has different orientation than the HRIT data I have (which has South up), so I need to reorient it in the code below.

from satpy import Scene
import glob

fnames = glob.glob("/home/lahtinep/data/satellite/geo/0deg/*202003201800*")
glbl = Scene(reader="", filenames=fnames)
glbl.load(["IR_108"])
lcl = glbl.crop(ll_bbox=(-10, 34, 25, 54))
adef = lcl["IR_108"].attrs["area"]
# Rotate the area definition. Do not do this if Native format has North up already
adef = adef[::-1, ::-1]
print(adef.area_extent)
# -> (-871617.152427787, 3385955.09987868, 2128786.126153232, 4745137.785075882)
print(adef.width)
# -> 1000
print(adef.height)
# -> 453

As the HRV channel has triple the resolution, width and height will need to be multiplied by 3 to get the full resolution. With these we get the new area:

seviri_0deg_europe_crop_1km:
  description: Full globe MSG image 0 degrees cropped to Europe, 1 km resolution
  projection:
    proj: geos
    lon_0: 0.0
    a: 6378169.0
    b: 6356583.8
    h: 35785831.0
  shape:
    height: 1359
    width: 3000
  area_extent:
    lower_left_xy: [-871617.152427787, 3385955.09987868]
    upper_right_xy: [2128786.126153232, 4745137.785075882]

Save this to your $PPP_CONFIG_DIR/areas.yaml, and then you can do this:

glbl.load(["HRV"])
lcl = glbl.resample("seviri_0deg_europe_crop_1km")
lcl.show("HRV")

[gersmit]

Dear Pnuu, Thanks for the explanation and the piece of script. This gives more clarity and try this out in the coming days.

Dear JD and Pnuu,

What I do achieve is to get a better and clearer picture of the Netherlands. In the appendix you see Europe with the Netherlands above (red circle). Now I think, if I can use the HRV images, the Netherlands will also become more visible. Or do you have another solution for me to achieve this? Working with the other images is perfect.

Thanks so far and look forward to hearing from you. Regards Ger

[sjoro]

just thinking about this... we do have fci_0deg_1km area defined in the areas.yaml. this is basically what the SEVIRI HRV channel is: 1km geos-projection with a 0 degree subsatellite longitude. so, a solution for this with "default" Satpy-installation, and hence skipping the $PPP_CONFIG_DIR, would be:

scn = Scene(filenames=[fname], reader="seviri_l1b_native")
scn.load(["HRV"])
res = scn.resample('fci_0deg_1km', radius_of_influence=3000)
eur = res.crop(ll_bbox=(-10, 34, 25, 54))
eur.show("HRV")

[gersmit]

Hi Sjoro,

I have learned a lot from previous posts, thanks for that.

Now I come to the conclusion that I am getting a memory error with HRV. I run on a Raspberry pi 4, this to make everything run automatically. Can I also use fci_0deg_1km with VIS008 (this go good) to get a good picture of the Netherlands. Or can you explain to me how to work to beautifully project the Netherlands.

Thanks in advance for message, Regards Ger

[sjoro]

if you want to stay in geos-projection, the code above would work for any SEVIRI channel. otherwise you can always define your own AreaDefinition in areas.yaml as @pnuu explained. i would encourage you to have a look at the eurol-area and increase it's resolution. it seems to be defined as 3km grid.

[gersmit]

Dear Sjono,

It all works fine with a small question: What are the values for in areas.yaml projection: lon_0: 0.0 a: 6378169.0 b: 6356583.8 h: 35785831.0

And is it possible to make the result from areas.yaml sharper after increase resolution?

Thanks Ger

[djhoese]

@gersmit The "projection" defines the coordinate space for your new image. This is kind of a big topic that I try to cover in the resampling portion of this tutorial. You can find more information specific to the geostationary projection here: https://proj.org/operations/projections/geos.html

Your options for making an image sharper:

1. Increase the resolution of the AreaDefinition being resampled to.
2. Change the resampling algorithm used. Nearest neighbor is the default ('nearest'). There are also a couple others. See the Satpy resampling documentation for more information.
3. Apply some sharpening algorithm on the data/image to make it look sharper. This isn't something we really provide and isn't likely to produce what you're looking for.
4. Use higher quality/resolution data. Obviously you don't have much choice here, but at the end of the day you can't make an image look way better without better quality input imagery. If there is something you're really hoping to see out of this imagery let us know.

[gersmit]

Hi DJ,

As I have written, I have learned a lot from all of you in recent days.

You have understood that HRV is present but the Raspberry Pi does not handle it. So back to Hrit and work okay. What I want to achieve is an image with the Netherlands in the middle (see attachment). I've got gets with Python Opencv for better picture but what you say the input determines it. I will continue to study what you explain with regard to AreaDefinition Nearest neighbor and should be able to come up with that.

Greetings and thanks Ger

[gersmit]

Dear DJ, I managed to make a map of the Netherlands (appendix), very nice. Only the measurements are not correct, difference in height and width, can I update that?

Ger

[djhoese]

How are you making the AreaDefinition? What does the AreaDefinition look like (the YAML text if you're doing it in the YAML file)?

[gersmit]

As you explane from the crop. Make Height and Width double value

seviri_0deg_nederland_crop_2km: description: Full globe MSG image 0 degrees cropped to Europe, 2 km resolution projection: proj: geos lon_0: 0.0 a: 6378169.0 b: 6356583.8 h: 35785831.0 shape: height: 848 width: 566 area_extent: lower_left_xy: [-220529.63268756866, 4163059.3925714493] upper_right_xy: [1051641.30961895, 5012173.488497734]

[djhoese]

You said:

Only the measurements are not correct, difference in height and width, can I update that?

Do you mean the number of pixels are not the same? Different number of columns and rows (you'd like a square image)?

This all ends up being with the math of the extents and the number of pixels.

pixel_size_x = (upper_right_xy[0] - lower_left_xy[0]) / width
pixel_size_y = (upper_right_xy[1] - lower_left_xy[1]) / height

So if you want the pixel size to be 2km (as your "description" says) then you'll need to do the math to figure out what the other values need to be (using the seviri_0deg as a starting point I guess). Or you could look at the create_area_def function described in the pyresample documentation I mentioned before.

[gersmit]

👍 I find out, thanks

[gersmit]

Dear DJHOESE One small other question if you can help me: All crop and other action i can plot with Matplotlib and get images with coastlines etc. Please can you tell me how i can manege that with composite? composite = compositor([local_scene[0.6], ... local_scene[0.8], ... local_scene[10.8]]) Can print to : from satpy.writers import to_image img = to_image(composite)

How can i make it possible the result from the composite plot to Mlatplotlib with lines.

Thanks for an answer Ger

[djhoese]

Where did you get that code that you are using? This tutorial goes through how you can configure a custom composite if that's what you want to do. The code you have doesn't look like anything from this tutorial, but if you assuming composite is an xarray.DataArray object (run print(type(composite)) to check) then:

from satpy.writers import get_enhanced_image
img = get_enhanced_image(img)
img.data.plot.imshow(rgb='bands')

Should be a good starting point. You can use this tutorial to fill in the rest of the information for how to get coastlines using cartopy.

[gersmit]

I have this: composite = 'snow' local_scn.load([composite]) local_scn.show(composite) This works, So far, so good. The problem is, how can i put the composite in the row: crs = local_scn["snow"].attrs["area"].to_cartopy_crs() I try all possiblities but always get error, only now, not with other plots.

And this: **from satpy.composites import GenericCompositor

compositor = GenericCompositor("overview") composite = compositor([local_scene[0.6], ... local_scene[0.8], ... local_scene[10.8]])** The same, how do i import composite in matplotlib

[djhoese]

how can i put the composite in the row:

What do you mean by "row"?

I try all possiblities but always get error

What error do you get?

[gersmit]

This is wat i wil make is make snow chart: **from satpy.writers import to_image composite = 'snow' local_scn.load([composite]) from satpy.writers import get_enhanced_image img = get_enhanced_image('snow') img.data.plot.imshow(img)

crs = local_scn["snow"].attrs["area"].to_cartopy_crs() ax = plt.axes(projection=crs) ax.coastlines() ax.gridlines() ax.set_global() plt.imshow(local_scn['snow'], transform=crs, extent=crs.bounds, origin='upper') plt.show()

What i get == : MemoryError or AttributeError: 'str' object has no attribute 'squeeze'

[djhoese]

MemoryError means your system doesn't have enough memory. This is likely since you are using matplotlib instead of saving to something like a geotiff where the data can be compressed and tiled. With matplotlib the whole data array has to be in memory at the same time.

For the AttributeError, you should be using img = get_enhanced_image(local_scn[composite]) instead (note how I'm passing the DataArray local_scn[composite] instead of the string "snow").

[gersmit]

Thanks, make a lot clear. Is there an other way to put country lines between local_scn.load and local_scn. save? Then i wil direct save and not plot. Otherwise, is it possible to make lines on a png image after save? I just bought new sd-card for the RPI and have to clear the RPI memory because al test are doing on this.

Thanks Ger

[djhoese]

You could reduce the resolution of the image by resampling or something like the aggregate method on the Scene (https://satpy.readthedocs.io/en/stable/api/satpy.html#satpy.scene.Scene.aggregate). Since you are saving to a PNG matplotlib plot you aren't actually going to get the full resolution of the data in the output image (a 800x600 pixel matplotlib image isn't going to show every pixel of a 10000x10000 input array).

Otherwise, there are satpy based solutions but they don't produce a plot like matplotlib, they "burn" the lines into the image. See the overlay and decorate options of get_enhanced_image: https://satpy.readthedocs.io/en/stable/api/satpy.writers.html?highlight=get_enhanced_image#satpy.writers.get_enhanced_image

However, if you are using the overlay and decorate kwargs then you should maybe just be using Scene.save_datasets (which also accepts these kwargs) to save your image as a PNG.

Hope this helps.

[gersmit]

Thanks DJ, i going to try and hope it wil give result. You hear from me

Ger

[gersmit]

Hi DJ, Everything works fine now for me as i wanted. Thanks

But in your good writen Satpy documents i found: From satpy.composites import StaticImageCompositor I save create a TIF file with img.save('image.tif') from : satpy.composites backgroundcompositor, but cannot load it with StaticImageCompositor. Error : No Supported Files Found.

Can you explain, or show me where i find artikel, how to handle with images?

Thanks Ger

[djhoese]

This is a known pain point and we are working on a better solution. The compositor that that is trying to use (StaticImageCompositor) relies on an image to exist on your computer before it can be loaded. See https://github.com/pytroll/satpy/issues/1362 for more information. If you download https://neo.sci.gsfc.nasa.gov/archive/blackmarble/2016/global/BlackMarble_2016_3km_geo.tif and put it in your current working directory then it should work. Alternatively, if you download it and put it in another directory, you can set the environment variable SATPY_ANCPATH=/path/to/image_dir and it should work. If this doesn't work please provide the full traceback/error you are seeing.

[gersmit]

That's a shame, because I want to import an image to use as a background in the CloudCompositor. The explanation with the CloudCompositor is namely: composites can be used as an overlay on top of e.g. static terrain images. Can you explain to me how I can do this differently with a png, tif or other kind of image? I have a routine in Python but I prefer to use your nice Satpy software because then everything stays together.

Thanks in advance for an answer, Regards Ger

[djhoese]

I'm not sure what you're asking. My previous comment told you how to do it with the StaticImageCompositor. If you want to see an example in the actual Satpy configs:

https://github.com/pytroll/satpy/blob/master/satpy/etc/composites/visir.yaml#L424-L427

You could do the same for your own custom composites if that's what you want to do. Put the PNG/tif in a directory and use SATPY_ANCPATH as I've described above to point to where you put the images.

[gersmit]

DJ, thanks for your explanation.

Sorry to ask for the known road and good that you still make it clear to me. Excuse me but for me it is like a whole new matter and I have been busy for a few weeks now and with your help I am starting to understand it quite well. As for the yaml files, that is searching for me. For you as super software techies easy but for me and possibly other amateurs it is figuring out the yaml files and how they should function in the whole of the scripts.

Sorry again and I will continue to investigate.

Regards Ger

[djhoese]

No problem. Some of the documentation for how these YAML files work is a little lacking.

The main issue @gersmit is I'm not sure what you are trying to do. Are you trying to create your own composite that doesn't exist in Satpy or are you trying to use one that already exists in Satpy? What are you trying to produce?

[gersmit]

I now know that just like areas.yaml also for compsites a yaml file should be used. I just compared it and now find out how it works. Fortunately, there is Google translate, which makes it a lot easier as technical English is difficult for some of us. What I want is the following: With "Maps with Cartopy", https://rabernat.github.io/research_computing_2018/maps-with-cartopy.html, it is easy to create backgroud maps in tif. I want to merge this with BackgroundCompositor. So in summary: create a map, import the tif with StaticImageCompositor and then merge in BackgroundCompositor. I don't think this can be difficult but does it make sense?

I hope you understand my thoughts.

[djhoese]

Ok. I think I got it. So you'll be creating your own custom composite. I would start with:

https://satpy.readthedocs.io/en/stable/composites.html#creating-composite-configuration-files

Then you can look at existing composite configurations for inspiration. For example this SEVIRI one:

https://github.com/pytroll/satpy/blob/master/satpy/etc/composites/seviri.yaml#L498-L503

This creates an image with an RGB composite in the foreground and a "night_background" composite in the background. The night background is defined here:

https://github.com/pytroll/satpy/blob/master/satpy/etc/composites/visir.yaml#L424-L427

This composite is the same as I linked to above and can use the SATPY_ANCPATH environment variable to point to any image you want. So following the instructions in the top link, you can create your own custom YAML files, then you can create your own foreground/background composite (in your own composites/visir.yaml or your own composites/seviri.yaml):

composites:
  your_map_background:
    compositor: !!python/name:satpy.composites.StaticImageCompositor
    standard_name: night_background_hires
    filename: your_background_image.tif
  your_rgb_name_with_background:
    compositor: !!python/name:satpy.composites.BackgroundCompositor
    standard_name: night_ir_with_background
    prerequisites:
      - true_color
      - your_map_background

You then can run the python:

import os
os.environ['PPP_CONFIG_DIR'] = '/path/to/your_custom_config_dir'
os.environ['SATPY_ANCPATH'] = '/path/to/your_background_images'

scn = Scene(reader='seviri', filenames=filenames)
scn.load(['your_rgb_name_with_background'])
# the rest of the code you had as usual for either saving to a geotiff or creating a plot with matplotlib

Hope this helps.

[gersmit]

Thanks @djhoese , lot of things to try. Keeps me from the street a couple of days. No problem, weather is bad in the Neterlands. 👍

[gersmit]

@djhoese , I am starting to understand and keep getting further, thanks for this. Smal question: Is there a posibility to rotate image for 180 degrees before ' img.show' or 'img.save'

from satpy.writers import to_image img = to_image(composite) img.invert([False, False, True]) img.stretch("linear") img.gamma(1.7) img.show() img.save('image.tif')

[djhoese]

No, but it depends what you need. Why does it need to be rotated? Does it need to be flipped (not rotated)? That sounds like a change in your area definition could make that work (swap the two "y" extents).

[gersmit]

With Matplotlib, plt.show is good, north is north With Satpy, from satpy.writers import to_image, img.show, north is south, strange but trough. So i have to rotate the image perhaps there is a , what you say FLIP or an rotate?

[gersmit]

@djhoese Very slowly I start to understand how everything should work. Thanks for that. Small question about a part of generic.yaml: Can you tell me how I can take the variable min_value and max_value from enhancements: colorized_ir_clouds to Python to use it in Matplotlob plot as value in colorbar.

**colorized_ir_clouds: standard_name: colorized_ir_clouds operations:

• name: colorize method: !! python / name: satpy.enhancements.colorize kwargs: palettes:
  • {colors: rdylgn, min_value: 233, max_value: 263}
  • {colors: grays, min_value: 264, max_value: 303.15}**

[djhoese]

There is no Satpy interface for doing this. You could load the yaml yourself with the yaml library. Something like:

import yaml
yaml_dict = yaml.load(open('/path/to/enhancements/generic.yaml', 'r').read(), loader=yaml.UnsafeLoader)
yaml_dict['enhancements']['colorized_ir_clouds']['palettes'][0]['min_value']

But I'm not sure this makes sense. It all depends on what you're trying to accomplish. Are you using this enhancement and then putting the result into matplotlib?

Edit: I think normally, you would create your own matplotlib colormap and plot (imshow) the regular un-enhanced data.

[gersmit]

Dear @djhoese For making plot in Matplotlib, cartopy, as ax1 = fig.add_subplot(111), i get the error: ValueError: Axes should be an instance of GeoAxes, got <class 'matplotlib.axes._axes.Axes'>

Do you no how i can solve that problem. Thanks @gersmit

[gersmit]

@djhoese in response to the previous question: When i use this script in Python: fig_1 = plt.figure(figsize=(10,10)) # create new image and assign the variable "fig_1" to it ax1 = fig_1.add_subplot(211) # add subplot to "fig_1" and assign another name to it ax1.set_title("My image1") # set title --> ax1.coastlines() # I GET ERROR ax1.imshow(clouds) # print image in grayscale

This is working and i can make more plots on one Figure. But can you explane to me what i have to do to put coastlines to it. When i use ax1.coastline i get the error : Traceback (most recent call last): File "/home/pi/eumetview/make_plots_in_figure.py", line 202, in ax1.gridlines() AttributeError: 'AxesSubplot' object has no attribute 'gridlines'

I hoop you have an answer

Ger