Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Resample HRV channel of Eumetsat image #2906

Open
PaulaPerezL opened this issue Sep 19, 2024 · 4 comments
Open

Resample HRV channel of Eumetsat image #2906

PaulaPerezL opened this issue Sep 19, 2024 · 4 comments

Comments

@PaulaPerezL
Copy link

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
Copy link
Member

pnuu commented Sep 19, 2024

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.

@PaulaPerezL
Copy link
Author

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
Copy link
Member

pnuu commented Sep 20, 2024

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"]

@PaulaPerezL
Copy link
Author

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

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

2 participants