pytroll/satpy

Resample HRV channel of Eumetsat image

Opened this issue · 4 comments

Describe the bug
I am trying to do a resample of the Eumetsat image ("EO:EUM:DAT:MSG:HRSEVIRI") in channel HRV with the aim of crop an area and I got a RuntimeWarning that makes my code goes slowly and I can not fix:

/opt/miniconda3/envs/helios/lib/python3.10/site-packages/dask/core.py:127: RuntimeWarning: invalid value encountered in sin

I try to create a mask if the image have nan values but the result is the same

To Reproduce

def load_and_crop_image(image_path, channel_to_load, area):
    if image_path.endswith(".nat"):
        reader = "seviri_l1b_native"
    elif image_path.endswith(".grb"):
        reader = "seviri_l2_grib"
    else:
        raise ValueError(f"Unknown file extension for {image_path}")
    
    with warnings.catch_warnings(record=True) as caught_warnings:
        warnings.simplefilter("always")
        scn = Scene(filenames={reader: [image_path]})
        scn.load(channel_to_load)
    for warning in caught_warnings:
        if "The quality flag for this file indicates not OK" in str(warning.message):
            log.info("The quality flag for this file indicates not OK")
            exit()
    if channel_to_load == ['HRV']:
        scene = resample_geostationary_hrv(scn, channel_to_load)
        cropped_scene = scene.crop(ll_bbox=GEOGRAPHIC_BOUNDS[area])
    else:
        cropped_scene = scn.crop(ll_bbox=GEOGRAPHIC_BOUNDS[area])

    return cropped_scene

def resample_geostationary_hrv(scn, channel_to_load):
    orb_params = scn[channel_to_load[0]].attrs['orbital_parameters']

    # Define the geostationary projection
    area_id = "Geostationary"
    description = "Geostationary Projection"
    proj_id = "geos"
    proj_dict = {"proj": "geos", "lon_0": orb_params['projection_longitude'], "h": orb_params['projection_altitude']}
    width = scn[channel_to_load[0]].shape[1]
    height = scn[channel_to_load[0]].shape[0]
    area_extent = AREA_EXTENT['HRV']

    area_def = geometry.AreaDefinition(area_id, description, proj_id, proj_dict, width, height, area_extent)

    # Resample the image to the geostationary projection
    scn_resampled = scn.resample(area_def)
    return scn_resampled


**Environment Info:**
 - OS: [e.g. OSX, Windows, Linux]
 - Satpy Version: '0.48.0'
 - PyResample Version: '1.30.0'
pnuu commented

Hi,

Your code isn't complete so it can't be run.

Do you resample the cropped scene? That should not be necessary and you can play around with the reduce_data and mask_area arguments in .resample() call to see if they affect the processing time. So just switch between True and False for different combinations. And resample directly the original Scene in which you load the channel data:

scn_resampled = original_scn.resample(area_def, reduce_data=True, mask_area=True)

In general the RuntimeWarning shouldn't affect the processing speed, it's just a warning that there are nan values in some call to sin(). Like in the space.

For me loading HRV, resampling to the built-in euro4 area and saving to geotiff takes roughly 5 seconds.

thank you for you quick response and your explanation. What I need to do is to crop the image for different countries, for example Spain, Great Britain, etc. With all channel except HRV I can just load the scene and crop it using the coordinates of the place I want to crop. However in the HRV channel, I think because of the grid (which is different from the other channels), I need to resample de image first and the crop it and it is in this process when I get the warning.
Again thank you very much for your time.

pnuu commented

Indeed. One option is to create an area definition from the other cropped channels and resample the HRV channel directly to that/those area/areas.

cropped_area = cropped_scene["IR_108"].attrs["area"]

Thank you very much for your help, I receive the same warning but I like more the way you proposed so I changed it :)