satpy: resize images


R. Alblas
 

Ernst,

Do you know if it  is possible to overwrite the shape numbers as defined in the areas.yaml?

Currently an extra entry in area.yaml is defined with different width/height of the same area. If this could be overwritten then no extra entries are needed, and with just one boolean in python scripts a lower resolution for movies and the right subdir (images versus frames) could be set. I couldn't find a function in the satpy docs to change the 'shape' numbers as defined in areas.yaml.

The other way, using -resize in convert, is much slower, as you will know.

Cheers,
Rob


Ernst Lobsiger
 

Rob,

you can read in some 'area' from the areas.yaml file and change width and height.
See the screen copy below what I just tested at the prompt not in a script so far:

(pytroll) eumetcast@kallisto:~$ python
Python 3.10.8 | packaged by conda-forge | (main, Nov 22 2022, 08:26:04) [GCC 10.4.0] on linux
Type "help", "copyright", "credits" or "license" for more information.
>>> from satpy.resample import get_area_def
>>> my_area=get_area_def('germ')
>>> print(my_area)
/home/eumetcast/miniconda3/envs/pytroll/lib/python3.10/site-packages/pyproj/crs/crs.py:1296: UserWarning: You will likely lose important projection information when converting to a PROJ string from another format. See: https://proj.org/faq.html#what-is-the-best-format-for-describing-coordinate-reference-systems
  proj = self._crs.to_proj4(version=version)
Area ID: germ
Description: Germany
Projection: {'a': '6378144', 'lat_0': '90', 'lat_ts': '50', 'lon_0': '5', 'no_defs': 'None', 'proj': 'stere', 'rf': '298.253168108487', 'type': 'crs', 'units': 'm', 'x_0': '0', 'y_0': '0'}
Number of columns: 1024
Number of rows: 1024
Area extent: (-155100.4363, -4441495.3795, 868899.5637, -3417495.3795)
>>> my_area.width=2048
>>> my_area.height=2048
>>> print(my_area)
Area ID: germ
Description: Germany
Projection: {'a': '6378144', 'lat_0': '90', 'lat_ts': '50', 'lon_0': '5', 'no_defs': 'None', 'proj': 'stere', 'rf': '298.253168108487', 'type': 'crs', 'units': 'm', 'x_0': '0', 'y_0': '0'}
Number of columns: 2048
Number of rows: 2048
Area extent: (-155100.4363, -4441495.3795, 868899.5637, -3417495.3795)
>>>

There are also different functions to define the whole area within the code.

Hope that helps,
Ernst


Ernst Lobsiger
 

On Fri, Jan 27, 2023 at 03:08 AM, R. Alblas wrote:
and with just one boolean in python scripts a lower resolution for movies and the right subdir (images versus frames) could be set.
Rob,

with what I proposed this will not work easily. In my scripts area is a string. If given to the resampler it reads in the area (object).
You AFAIK can give the resampler the area object but area as string is also used in the legend and the file names. So this will ask
for some extra coding in GEOstuff.py The use of subdirectory sat/frames versus subdirectory sat/images is not well defined either.

Cheers,
Ernst


R. Alblas
 

Ernst,

Thanks, this works, my_area can easily be inserted in global_scene.resample:

        my_area=get_area_def(area)
        my_area.width=1080
        my_area.height=1080
        new_scene = global_scene.resample(my_area, **resampler_kwargs)  # **defined @ top of module

And the generated image is indeed smaller

But... it's still much slower than defining width/height in areas.yml:

3km: 38 secs
rss: 15 secs
3km with adapted area ('my_area'): 29 secs

In this test I left in the extra code also for msg_fes_rss.

This time difference is especially in new_scene.save_dataset(), where I wouldn't expect this (but probably the actual resample computation is done in here?).

As far as I can see the only place where 'area' is not used just as text is in resample(area,...) . Strange...

Rob.


On 27-01-2023 12:53, Ernst Lobsiger via groups.io wrote:

Rob,

you can read in some 'area' from the areas.yaml file and change width and height.
See the screen copy below what I just tested at the prompt not in a script so far:

(pytroll) eumetcast@kallisto:~$ python
Python 3.10.8 | packaged by conda-forge | (main, Nov 22 2022, 08:26:04) [GCC 10.4.0] on linux
Type "help", "copyright", "credits" or "license" for more information.
>>> from satpy.resample import get_area_def
>>> my_area=get_area_def('germ')
>>> print(my_area)
/home/eumetcast/miniconda3/envs/pytroll/lib/python3.10/site-packages/pyproj/crs/crs.py:1296: UserWarning: You will likely lose important projection information when converting to a PROJ string from another format. See: https://proj.org/faq.html#what-is-the-best-format-for-describing-coordinate-reference-systems
  proj = self._crs.to_proj4(version=version)
Area ID: germ
Description: Germany
Projection: {'a': '6378144', 'lat_0': '90', 'lat_ts': '50', 'lon_0': '5', 'no_defs': 'None', 'proj': 'stere', 'rf': '298.253168108487', 'type': 'crs', 'units': 'm', 'x_0': '0', 'y_0': '0'}
Number of columns: 1024
Number of rows: 1024
Area extent: (-155100.4363, -4441495.3795, 868899.5637, -3417495.3795)
>>> my_area.width=2048
>>> my_area.height=2048
>>> print(my_area)
Area ID: germ
Description: Germany
Projection: {'a': '6378144', 'lat_0': '90', 'lat_ts': '50', 'lon_0': '5', 'no_defs': 'None', 'proj': 'stere', 'rf': '298.253168108487', 'type': 'crs', 'units': 'm', 'x_0': '0', 'y_0': '0'}
Number of columns: 2048
Number of rows: 2048
Area extent: (-155100.4363, -4441495.3795, 868899.5637, -3417495.3795)
>>>

There are also different functions to define the whole area within the code.

Hope that helps,
Ernst


Ernst Lobsiger
 

Rob,

there are fcomposites and tcomposites. The resampling is normally delayed
and done under fcomposites IIRC in compute_writer_results(save_objects).

Ernst


R. Alblas
 

A, I see, overlooked the tcomposite part. Now it works

Thanks,

Rob.

On 27-01-2023 14:24, Ernst Lobsiger via groups.io wrote:

Rob,

there are fcomposites and tcomposites. The resampling is normally delayed
and done under fcomposites IIRC in compute_writer_results(save_objects).

Ernst


Ernst Lobsiger
 

On Fri, Jan 27, 2023 at 06:51 AM, R. Alblas wrote:
A, I see, overlooked the tcomposite part. Now it works
Rob and All,

Yes, tuning Satpy for top speed in SPS 4.x makes things less obvious:

tcomposites means generate=True (default, the EMC 3.0 way). If you
have a list of composites, all channels of a composite are resampled
and the composite is generated. Then the next composite is treated the
same way. This means that common channels are resampled multiple
times, which of course is much slower with long lists of composites.

fcomposites means generate=False, if you have a list of composites
(that use common channels), every channel is resampled only once.
The composites are generated delayed, just mixed up at the end.

That said, if you just do ONE composite it's probably the same speed
whether you evaluate this with generate=True or with generate=False.

Unfortunately not all composites can be treated as fcomposites.
GEO composites with background files do only work as tcomposites.
And there are even some LEO composites that need generate=True.

Whether a certain composite is treated with True=1 or False=0
stands as a last number of the  instrument_abr{}  dictionaries.

New and unknown composites are treated with generate=True.

Cheers,
Ernst


R. Alblas
 

There is something very strange going on with resizing images.

To make sure that this is not because something is wrong in my scripts I tried the original scripts of Ernst:

MSG2-msg_ioc_rss.py

That gives the expected image, with reduced size (width/height defined 1080x1080 in areas.yaml).

Now I did just adapt in areas.yaml the size to 1081x1081 and I get following result:


This is the result before magick, so the png in tmpdirs/xmsg2.

A size of 1079x1079 is OK. Looks like everything between 1081x1081 and 1200x1200 gives this result; above 1200 and under 1081 it's OK.

Note that all numbers in areas.yaml for msg_ioc_3km and msg_ioc_rss are the same, except the shape-numbers.

Looks like resizing using the shape numbers in areas.yaml not always gives the expected results, or maybe it's not the right thing to resize this way???

Rob.


On 27-01-2023 17:50, Ernst Lobsiger via groups.io wrote:

On Fri, Jan 27, 2023 at 06:51 AM, R. Alblas wrote:
A, I see, overlooked the tcomposite part. Now it works
Rob and All,

Yes, tuning Satpy for top speed in SPS 4.x makes things less obvious:

tcomposites means generate=True (default, the EMC 3.0 way). If you
have a list of composites, all channels of a composite are resampled
and the composite is generated. Then the next composite is treated the
same way. This means that common channels are resampled multiple
times, which of course is much slower with long lists of composites.

fcomposites means generate=False, if you have a list of composites
(that use common channels), every channel is resampled only once.
The composites are generated delayed, just mixed up at the end.

That said, if you just do ONE composite it's probably the same speed
whether you evaluate this with generate=True or with generate=False.

Unfortunately not all composites can be treated as fcomposites.
GEO composites with background files do only work as tcomposites.
And there are even some LEO composites that need generate=True.

Whether a certain composite is treated with True=1 or False=0
stands as a last number of the  instrument_abr{}  dictionaries.

New and unknown composites are treated with generate=True.

Cheers,
Ernst


Ernst Lobsiger
 

On Sat, Jan 28, 2023 at 04:08 AM, R. Alblas wrote:
Looks like resizing using the shape numbers in areas.yaml not always gives the expected results, or maybe it's not the right thing to resize this way???
Rob,

I have resized msg_fes_3km in my ../userconfig/areas.yaml to msg4_s1081 = 1081x1081, msg4_s1082 = 1082x1082, msg4_s1083 = 1083x1083,
msg4_s1084 1084x1084. Everything worked as expected with exact heights of the resulting images and increased width due to the added legend.

I wonder whether you are messing with the cache. Please turn off all caching. Maybe you also have to change the IDs (ID=name (?) as I did).
When resizing areas inside your scripts you probably also have to take care that you use calculated and rounded integers for width and height.

Cheers,
Ernst


R. Alblas
 


I wonder whether you are messing with the cache. Please turn off all caching. Maybe you also have to change the IDs (ID=name (?) as I did).
When resizing areas inside your scripts you probably also have to take care that you use calculated and rounded integers for width and height.

Cheers,
Ernst

Ernst,

I have now switched off caching:

GEOcache = False
OVRcache = False

And also cleaned the cache: SPSdata/cache

Now adapting width/height in areas.yaml does work properly.

But resizing in script still doesn't work properly. What I did: in MSG2-msg_ioc_rss.py, I added:


  area = 'msg_ioc_3km'
  from satpy.resample import get_area_def
  sarea=get_area_def(area)
  print(sarea)
  sarea.width=1100
  sarea.height=1100
  sarea.area_id='resized_1100'
  area=sarea.area_id
  print(sarea)

sarea goes to global_scene.resample(), 2x (for fcomposite and tcomposite)

area is going to the rest (only used for text, creating filenames etc. as far as I can see)

Result:

Overlays: borders etc. are resized properly.
But satellite data is not resized, see image below.

This image is from SPSdata/tmpdirs/xmsg2 (for this email png converted to jpg to reduce file size). I have made sure that this is indeed the file generated, it's easy to mix-up things...

And again checked that cache is off, and nothing is generated in SPSdata/cache.

Comparing the change by doing print(sarea) before and after change (see above) shows the correct numbers, i.e.:

Area ID: resized_1100
Number of columns: 1100
Number of rows: 1100

and all other numbers are the same.

What's going on here...???




Ernst Lobsiger
 

On Mon, Jan 30, 2023 at 01:33 AM, R. Alblas wrote:

What's going on here...???

 

Rob,

My first idea was that when you read in the area its ID is retained and the projection cache is scanned for this ID.  Then resampling 1100x1100 would just show the top left part of the original 3712x3712. Now if this also happens with both caches off, then changing the shape this way for a global inside GEOstuff.py might not be the right thing to do (or reading in with get_are_def() already does something behind the scene). There is a difference with your first picture where everything was contained in a circle probably due to the application of a cached OVR file not transparent in free space. A last idea is to try different resamplers changing resampler_kwargs on top of the module.

Maybe that's a question for Martin or Dave on the google list. The developers are usually ready to help rather fast.

Regards,
Ernst


Ernst Lobsiger
 


R. Alblas
 


My first idea was that when you read in the area its ID is retained and the projection cache is scanned for this ID.  Then resampling 1100x1100 would just show the top left part of the original 3712x3712. Now if this also happens with both caches off, then changing the shape this way for a global inside GEOstuff.py might not be the right thing to do (or reading in with get_are_def() already does something behind the scene). There is a difference with your first picture where everything was contained in a circle probably due to the application of a cached OVR file not transparent in free space. A last idea is to try different resamplers changing resampler_kwargs on top of the module.

Ernst,

I don't think it has something to do with globals:

sarea=get_area_def('...')
sarea goes into geo_images() as parameter
in geo_images() sarea goes into global_scene.resample() as parameter


Btw, in geo_images() I did also do a print of height and width after:

            layers, height, width = new_scene[fcomposites[0]].shape

and it shows the right adapted values.


Maybe that's a question for Martin or Dave on the google list. The developers are usually ready to help rather fast.

I'll give it a try in the pytroll group.

Cheers,

Rob.



Regards,
Ernst


Ernst Lobsiger
 

Rob,

I can reproduce your problem here . The developers usually ask for a slimmed down script (<= 20 lines) to reproduce the issue.

Good luck,
Ernst


R. Alblas
 

Ernst,

Solution is:

  sarea=get_area_def(area).copy(width=1100,height=1100)
and that works.

Apparently something is cached, and with an explicit copy that's overwritten.

Another suggestion they gave:

  sarea=get_area_def(area).aggregate(x=3,y=3)

to resize with a factor 3. Which also works

Thanks,

Rob.

On 30-01-2023 14:25, Ernst Lobsiger via groups.io wrote:

Rob,

I can reproduce your problem here . The developers usually ask for a slimmed down script (<= 20 lines) to reproduce the issue.

Good luck,
Ernst


Ernst Lobsiger
 

On Mon, Jan 30, 2023 at 06:30 AM, R. Alblas wrote:
  sarea=get_area_def(area).copy(width=1100,height=1100)
and that works.
Rob,

yes, Martin was your man. As I suspected something (area caching)  is going on behind the scene. It works here.

Cheers,
Ernst


Ernst Lobsiger
 

Rob,

I see two possibilities to handle that without breaking all SPS 4.1 GEO scripts:

1)
you can define a last (optional) positional parameter resize=false in function
geo_images(...) which when specified would scale to some defined monitor height.

2)
you can define 'resize' as an additional keyword argument of function
geo_images(...) and treat it accordingly when specified.

Cheers,
Ernst


R. Alblas
 


On 30-01-2023 16:05, Ernst Lobsiger via groups.io wrote:
Rob,

I see two possibilities to handle that without breaking all SPS 4.1 GEO scripts:......

Ernst,

I do it a bit different. In (e.g.) MSG.py a boolean 'for_movie' is defined to resize or not (which is driven by a command line argument).

In case of a resize, sarea is created, and at the same time subdir is changed from 'images' to 'frames'.

sarea goes into geo_images() the same way as 'area' in the original; in geo_images() it is detected if the parameter 'area' is a string or not:

    if isinstance(sarea,str):

If it's not a string then sarea.area_id is going to get_overlays() etc. and sarea is going to resample().


That way it's backwards compatible and an extra argument isn't even necessary. (I keep GEOstuff.py as close to the original as possible).
Of course there are many possibilities to solve this; thanks for sharing yours.

Cheers,

Rob.