flipped seviri image using pyresample

[Bruno Picard]

Hi all,

I hope this is the right channel to ask this question.

If it is not, my apologies, just let me know where can I ask this question in a more appropriate location.

I'm Bruno Picard, I'm familiar with satpy, pyresample and the pytroll group.

I'm currently involved in the SWOT mission commissioning phase, using some precipitation product to assess the impact of rain

Note that this may be a well-known question but I was not able to find any solution.

but it may be related to

https://github.com/pytroll/satpy/issues/1160

also described here

https://gist.github.com/sfinkens/a07ac7982699ed8fc0fb433bb75818b7

Objective

The goal is to map a H60 gridded product distributed by Eumetsat Hydro-SAF

the precipitation products are listed here
- https://hsaf.meteoam.it/Products/ProductsList?type=precipitation
the H60 products are detailed here
- https://hsaf.meteoam.it/Products/Detail?prod=H60B
 

H60 precipitation gridded products are described as
Instantaneous precipitation maps generated combining geostationary (GEO) IR images from operational geostationary satellites 'calibrated' by precipitation measurements from MW images on Low Earth Orbit (LEO) satellites, processed soon after each acquisition of a new image from GEO.

They are based on Seviri 0deg

Issue

Using pyresample to compute x,y coordinates and longitude/latitudes coordinates results in a fliplr / flipud image

Description

the product does not include the coordinates. lon/lat are available in a separate file.

I try to compute lon/lat using pyresample but failed to have the correct sorting or columns/lines

# read the dataset

fname = 'h60_20230207_0730_fdk.nc'

ds = xr.open_dataset(hdle)  

# load the area def

area_def = get_area_def('seviri_0deg')

# get native coordinate

x, y = area_def.get_proj_coords()

# load crs

source_crs = CRS.from_proj4(ds.attrs['gdal_projection'])

# =  '+proj=geos +a=6378169.000000 +b=6356583.800000 +lon_0= 0.000000 +h=35785831.000000 +sweep=y'

# compute lon lat from x, y using Proj

geos = Proj(source_crs, preserve_units=False)
lat, lon = geos(x, y,inverse=True)

but the lat/lon grid are just messed up.

I would need to flipud/fliplr the data grid (precipitation) to have correct map.

I would appreciate any help on this,

thank you

Bruno

[Raspaud Martin]

Salut Bruno, 

I'm not familiar with the h60 product, but here is my take on the problem. The seviri_0deg area you are loading is describing the area of level 1 seviri data, which is flipped ud and lr by design, with the x and y coordinates being reversed. My suspicion is that the h60 product corrected this to present a product in the upright orientation. 

By the way, you can directly do area_def.get_lonlats() to get longitude and latitude of an area definition. 

Now, onto a solution: I would suggest you try to create an area definition from the data you have in the file. Apparently you have the proj string/crs, and if you have the extends (lower left and upper right corners in projection coordinates) you can construct the area definition you need. 

Do you think that is possible with this product? 

Best regards 

Martin 

Sent from my Galaxy

--
You received this message because you are subscribed to the Google Groups "pytroll" group.
To unsubscribe from this group and stop receiving emails from it, send an email to pytroll+u...@googlegroups.com.
To view this discussion on the web, visit https://groups.google.com/d/msgid/pytroll/a3cbcfe9-2fcc-4ade-a6cb-950baa2f3663n%40googlegroups.com.

[Bruno Picard]

Salut Martin,

As you can see, one of my one-post-every-two-years brought me here :-D

You probably have it right, thank you for the info.

The thing is, there is not so many information in the netcdf file

Attributes:
    hsaf60_algorithm_version:  2.0
    parallax_correction:       Mode_on
    satellite_identifier:      MSG4
    sub-satellite_longitude:    0.0f
    satellite_altitude:        35785831
    satellite_altitude_unit:   m
    r_eq:                      6378.169000
    r_eq_unit:                 Km
    r_pol:                     6356.583800
    r_pol_unit:                Km
    cgms_projection:           +proj=geos +coff=1856.000000 +cfac=13642337.00...
    gdal_projection:           +proj=geos +a=6378169.000000 +b=6356583.800000...

I did try some  by-hand area definition but je tourne en rond since my main source of information is satpy L1B area configuration

# from https://holmdk.github.io/2020/01/19/geostationary_satellite_img.html
lower_left_xy  = [-5570248.686685662, -5567248.28340708]
upper_right_xy = [ 5567248.28340708,   5570248.686685662]
# consistency with satpy area definition
lower_left_xy  = [-5570248.4773392612, -5567248.0741734440]
upper_right_xy = [+5567248.0741734440, +5570248.4773392612]
# ~ https://gist.github.com/sfinkens/a07ac7982699ed8fc0fb433bb75818b7
# ~ from satpy.resample import get_area_def

 area_def = get_area_def('seviri_0deg')

# ~ Area extent: (-5570248.6867, -5567248.2834, 5567248.2834, 5570248.6867)

every thing above is very similar, and then I tried

source_crs = CRS.from_proj4(ds.attrs['gdal_projection'])

crs_dict = source_crs.to_dict()

projection_name = crs_dict['proj']
lon_0 = crs_dict['lon_0']
a     = crs_dict['a']
b     = crs_dict['b']
h     = crs_dict['h']
proj  = crs_dict['proj']

area_def = AreaDefinition('areaD', projection_name, 'areaD',
                          {'lon_0': lon_0,
                           'a': a,
                           'b': b,
                           'h': h,
                           'proj': proj},
                          height, width,
                          (lower_left_xy[0], lower_left_xy[1],
                          upper_right_xy[0], upper_right_xy[1]))

but it does not solve anything.

Should I swap the area_extend somehow ? any clue ?

merci !

Bruno

[Raspaud Martin]

Yes, just swap the lower left and upper right corners in the area extent, that should do the trick.

aex = (aex[2], aex[3], aex[0], aex[1])

But I see you have the cfac and coff (and probably lfac and loff too) in your file, you could use that to derive the area extent. There is a description in the cgms 03 document I believe, but you can also look here:

https://github.com/pytroll/satpy/blob/main/satpy/readers/_geos_area.py#L63

Maybe that function can help? 

Anyway, glad you are still using satpy :) 

Best regards, 

To view this discussion on the web, visit https://groups.google.com/d/msgid/pytroll/90844e48-242a-4a94-8565-9b0a873ef2fan%40googlegroups.com.

[Bruno Picard]

Yes, it was perfect guidance, see the code below.

If it can help others, I was asked to share the code so here is the gitlab public repository

https://gitlab.com/bpicard_fluctus/h60_howto

and also a towards science note (I wanted to give it a try since long)

https://medium.com/@bpicard_43779/mapping-and-clipping-eumetsat-h-saf-h60-precipitation-products-defined-on-a-geostationary-grid-6656c6a60e1f

Thanks again !

# maybe there is already a parser somewhere in satpy or pyproj ?
parse_dict = {}
for keyval in ds.attrs['cgms_projection'].split():
    try:
        key,val = keyval.split('=')
        parse_dict[key] = val
    except:
        pass
pdict = {}
pdict['scandir'] = 'N2S'
pdict['h']       = float(parse_dict['+h'])*1000 - float(parse_dict['+r_eq'])*1000
pdict['loff']    = float(parse_dict['+loff'])
pdict['coff']    = float(parse_dict['+coff'])
pdict['lfac']    = float(parse_dict['+lfac'])
pdict['cfac']    = float(parse_dict['+cfac'])
pdict['ncols']   = ds.x.size
pdict['nlines']  = ds.y.size
# ~ (-5567248.28340708, -5567248.28340708, 5570248.686685662, 5570248.686685662)
area_extent =   get_area_extent(pdict)
# but this is for seviri L1B which is swapped
# so it should be
# (area_extent[2],area_extent[3],area_extent[0],area_extent[1])
# (5570248.686685662, 5570248.686685662, -5567248.28340708, -5567248.28340708)

area_def = AreaDefinition('areaD', projection_name, 'areaD',
                          {'lon_0': lon_0,
                           'a': a,
                           'b': b,
                           'h': h,
                           'proj': proj},
                          height, width,

                          (area_extent[2],area_extent[3],area_extent[0],area_extent[1])
                          )
# get coordinates
x, y     = area_def.get_proj_coords()
lon, lat = area_def.get_lonlats(nprocs=7)
lon[np.isinf(lon)] = np.nan
lat[np.isinf(lat)] = np.nan
# set lat/lon to ds
ds['latitude']  = (('x','y'),lat)
ds['longitude'] = (('x','y'),lon)