Fix projection handling for NWC GEO v2025 data

[pnuu]

This PR fixes the projection handling of NWC SAF GEO v2025 data in nwcsaf-geo reader. The a, b and h parameters of the proj_string and the area extent are now scaled (de-scaled?) with the hard-coded value of projection_definition:perspective_point_height, which can't be accessed from the attributes given by xr.open_dataset(). The struct(?) can't be accessed with netCDF4, h5netcdf nor h5py either (in a simple way in any case) with the use.

@strandgren can you give this a go to see how well the geolocation actually matches?

• Closes nwcsaf-geo reader generates wrong AreaDefinition for NWC-SAF v2025 FCI data #3366
• Tests added
• Fully documented
• Add your name to AUTHORS.md if not there already

[pnuu]

Oh, almost forgot: the code completion from GitHub Copilot helped with the first iteration of _get_v2025_proj_str_and_scale() method.

[codecov]

Codecov Report

❌ Patch coverage is 97.87234% with 1 line in your changes missing coverage. Please review.
✅ Project coverage is 96.32%. Comparing base (6c571f2) to head (5ac3ee6).
⚠️ Report is 77 commits behind head on main.

Files with missing lines	Patch %	Lines
satpy/readers/nwcsaf_nc.py	96.55%	1 Missing ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##             main    #3367      +/-   ##
==========================================
- Coverage   96.35%   96.32%   -0.04%     
==========================================
  Files         466      466              
  Lines       59083    59119      +36     
==========================================
+ Hits        56931    56945      +14     
- Misses       2152     2174      +22     

Flag	Coverage Δ
behaviourtests	3.58% <0.00%> (-0.01%)	⬇️
unittests	96.41% <97.87%> (-0.04%)	⬇️

Flags with carried forward coverage won't be shown. Click here to find out more.

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[strandgren]

Thanks for fixing! I tested the code and it works as expected 🙂

See one minor comment/suggestion in-line. Also, any particular reason why you now decided to round the numbers to three digits?

[pnuu]

Also, any particular reason why you now decided to round the numbers to three digits?
Mainly for consistency, that's what the other proj strings had. Also it'll lessen the changes of floating point precision errors in CRS equality comparisons.

[strandgren]

Looks good to me now :) I leave it to you to decide whether we should wait for a reply from NWC-SAF regarding the consistency across missions before merging or not. Probably that could be dealt with in a new PR if needed.

[pnuu]

They actually replied "please create a ticket at NWC SAF website after login". I don't have an account there, and vote for merging.

[strandgren]

I have created the ticket on the nwcsaf login page now, but I agree that this PR can be merged independently of this.

[strandgren]

Got response from NWC-SAF now.

My question:

Dear NWC-SAF SW team,

I noticed that the GEO satellite projection information and their units stored for FCI in the products generated using the v2025 software is quite different compared to the projection information stored in SEVIRI products generated with the v2021 software:

FCI v2025:
:gdal_projection = "+proj=geos +a=0.1782279581 +b=0.1776303847 +lon_0=0.000000 +h=1.000000 +sweep=y" ;
:gdal_geotransform_table = -0.1555898f, 5.588715e-05f, 0.f, 0.1555898f, 0.f, -5.588715e-05f ;
:gdal_xgeo_up_left = -0.1555898f ;
:gdal_ygeo_up_left = 0.1555898f ;
:gdal_xgeo_low_right = 0.1555898f ;
:gdal_ygeo_low_right = -0.1555898f ;

SEVIRI v2021
:gdal_projection = "+proj=geos +a=6378137.000000 +b=6356752.300000 +lon_0=0.000000 +h=35785863.000000 +sweep=y" ;
:gdal_geotransform_table = -5570249.f, 3000.403f, 0.f, 5570249.f, 0.f, -3000.403f ;
:gdal_xgeo_up_left = -5570249.f ;
:gdal_ygeo_up_left = 5570249.f ;
:gdal_xgeo_low_right = 5567248.f ;
:gdal_ygeo_low_right = -5567248.f ;

So it seems like the projection parameters in the FCI product have been normalized by the satellite height over ground in meters? Is this an expected behavior and something that is also applicable to SEVIRI and other GEO sensors in the v2025 software?

I also see that in the FCI products we have the new dataset called projection_definition :

int projection_definition ;
????????????projection_definition:grid_mapping_name = "geostationary" ;
????????????projection_definition:latitude_of_projection_origin = 0. ;
????????????projection_definition:longitude_of_projection_origin = 0. ;
????????????projection_definition:semi_major_axis = 6378137. ;
????????????projection_definition:semi_minor_axis = 6356752. ;
????????????projection_definition:perspective_point_height = 35786400. ;
????????????projection_definition:proj4 = "+proj=geos +a=0.1782279581 +b=0.1776303847 +lon_0=0.000000 +h=1.000000 +sweep=y" ;
????????????projection_definition:sweep_angle_axis = "y" ;

Is this also something which is available for all GEO sensors in the v2025 software?

Finally, I see some minor differences between the SEVIRI and FCI files (v2021 vs v2025) regarding the semi-minor axis and height, is that also expected?

The reason for my enquiry is that we need to adapt the reader in the Pytroll/Satpy library in order to properly support NWC-SAF SW FCI data and therefore need to understand the expected behavior and scope of the changes.

Many thanks in advance!

Kind regards,
Johan Strandgren
EUMETSAT

Their response:

Dear Johan
Your appreciations are correct.
We have changed the geostationary projection units from meters to radians. So, the projection parameter units have changed, and so have the corner units. We did this to follow the CF convention recommendations. In these units, h=1.0. If you want to go back to meters you need to multiply by the value of h in meters, which you find in projection_definition:perspective_point_height.

There is also a new variable called projection_definition which has all the values of the projection. As you can see, for convenience, some of them in meters, others in radians. This variable is linked to the product parameter variables (for example cma: cloud mask) via an attribute called "grid_mapping". So, the variable, cma for example, has an attribute called grid_mapping, which has as content a string with the name of a variable that defines the projection (in our case we called it projection_definition). Again, this is a requirement from the CF conventions, in order to be 100 per cent CF compliant.

The size of the Earth for SEVIRI and FCI _SHOULD_ come from the size of the Earth adopted by EUMETSAT ground processing when generating the level 1.5 data. This size of the Earth (semi minor and semi major axis) are needed to have the exact same projection as the EUMETSAT product produced by the ground segment.

Finally, remind you that the variables have a new dimension added to the variables which is time. So cma, for example, would have shape 1 x 1200 x 2300 for v2025 and shape 1200 x 2300 for previous versions. Again, following 100 per cent CF conventions.
And, yes, all these changes _SHOULD_ apply to all NWC SAF GEO products regardless of the satellite input (MTG, MSG, GOES, Himawari).
Any other questions, please let us know.

So a few takeaways relevant for this PR:

1. Since they are aligning with CF this will probably not change back to the previous behavior and therefore it's probably not unique to v2025. Hence, it may be better to check if the version is less than 2025 and do a legacy derivation in that case? Rather than having a method for v2025 specifically? Otherwise we'll face the same issue with the next major version update.

2. @pnuu you said you had issues accessing the projection_definition and ended up hard-coding projection_definition:perspective_point_height. But from their reply it seems like we can access it from the dataset instead via the attribute grid_mapping (e.g. ct:grid_mapping) - perhaps worth testing?

[pnuu]

The version is given in global attributes in source variable. In the case of v2025 data I have, the string is "NWC/GEO version v2025.1". I'll have to check the other file version I have if they follow the same pattern.

With xr.open_dataset() that we use currently to read the data, I see the grid_mapping of a dataset to be just a string:

In [7]: nc["ct"].attrs["grid_mapping"]
Out[7]: 'projection_definition'

Same with netcDF4.File():

In [15]: nc["ct"].grid_mapping
Out[15]: 'projection_definition'

Also ncdump gives only the same string:

$ ncdump -h S_NWC_CT_MTI1_MSG-N-NR_20250923T130000Z.nc | grep "ct:grid_mapping"
		ct:grid_mapping = "projection_definition" ;

[pnuu]

The version identifier with different versions of the data:

$ ncdump -h v2018/S_NWC_CT_MSG4_MSG-N-VISIR_20220203T104500Z.nc | grep source
		:source = "NWC/GEO version v2018.1" ;

$ ncdump -h v2021/S_NWC_CT_MSG3_MSG-N-VISIR_20231211T073000Z.nc | grep source
		:source = "NWC/GEO version v2021.2" ;

$ ncdump -h v2021.3/S_NWC_CT_MSG3_MSG-N-VISIR_20250206T130000Z.nc | grep source
		:source = "NWC/GEO version v2021.3" ;

So something like str(source).split(" ")[-1] >= "v2025" should work.