[gdal-dev] question re netcdf grid_mapping

[Alan Snow]

Hi Michael,

I don't claim to be an expert in the projection you are using or what the
GDAL team would like to do in this scenario, but hopefully the information
presented here will be helpful to you.

The grid mapping name "Oblique Mercator (LV95 - CH1903+)" in the netCDF
file does not appear to be in the list of supported grid mappings:
http://cfconventions.org/cf-conventions/cf-conventions.html#appendix-grid-mappings

Using pyproj to dump EPSG;2056 to CF parameters (
https://pyproj4.github.io/pyproj/latest/build_crs_cf.html), the grid
mapping name is "oblique_mercator":
>>> from pyproj import CRS
>>> crs = CRS("EPSG:2056")
>>> crs.to_cf()
{'azimuth_of_central_line': 90.0,
 'false_easting': 2600000.0,
 'false_northing': 1200000.0,
 'geographic_crs_name': 'CH1903+',
 'grid_mapping_name': 'oblique_mercator',
 'horizontal_datum_name': 'CH1903+',
 'inverse_flattening': 299.1528128,
 'latitude_of_projection_origin': 46.95240555555556,
 'longitude_of_prime_meridian': 0.0,
 'longitude_of_projection_origin': 7.439583333333333,
 'prime_meridian_name': 'Greenwich',
 'projected_crs_name': 'CH1903+ / LV95',
 'reference_ellipsoid_name': 'Bessel 1841',
 'scale_factor_at_projection_origin': 1.0,
 'semi_major_axis': 6377397.155,
 'semi_minor_axis': 6356078.962818189}

If the pyproj version is correct, I am wondering if the appropriate action
in this scenario is to contact the provider of the netCDF file to look into
and potentially address the issue?

If the swiss oblique mercator differs from the regular oblique mercator and
is not correctly represented by the current grid mappings in the CF
conventions, it may be worth opening an issue here to add it to the CF
conventions: https://github.com/cf-convention/cf-conventions. Once the
change is accepted there, I imagine a PR to GDAL will likely be welcome.

Best,
Alan

On Tue, Mar 1, 2022 at 6:57 PM Michael Sumner <mdsumner at gmail.com> wrote:

> Hello, I'm looking at the constants for CF_PT_* projection families in
>
> frmts/netcdf/netcdf_cf_constants.h
>
> I have a grid_mapping variable with attributes:
>
>    # float swiss_lv95_coordinates ;
>        #          swiss_lv95_coordinates:_FillValue = -1.f ;
>        #          swiss_lv95_coordinates:grid_mapping_name = "Oblique
> Mercator (LV95 - CH1903+)" ;
>        #          swiss_lv95_coordinates:longitude_of_projection_center =
> 7.43958333 ;
>        #          swiss_lv95_coordinates:latitude_of_projection_center =
> 46.9524056 ;
>        #          swiss_lv95_coordinates:false_easting = 2600000. ;
>        #          swiss_lv95_coordinates:false_northing = 1200000. ;
>        #          swiss_lv95_coordinates:inverse_flattening = 299.1528128 ;
>        #          swiss_lv95_coordinates:semi_major_axis = 6377397.155 ;
>
> and by my understanding this would require something like
>
> #define    CF_PT_LV95     "Oblique Mercator (LV95 - CH1903+)"
>
> and then attendant logic in the driver for that case, this is "EPSG:2056".
>
> It seems like this is a very long project of collating specific values
> used for a crs name ... where the more general case would use those
> parameters to build an Oblique Mercator - but this is Swiss oblique
> mercator  with a simpler set of initialization params.
>
> Is it a case of needing a very specific else if to match this case, or is
> there some more to this?
>
> Is it worth contributing cases as PRs?
>
> (I'm inclined to advise "use -a_srs 'EPSG:2056'". )
>
> Cheers, Mike
>
> Michael Sumner
> Software and Database Engineer
> Australian Antarctic Division
> Hobart, Australia
> e-mail: mdsumner at gmail.com
> _______________________________________________
> gdal-dev mailing list
> gdal-dev at lists.osgeo.org
> https://lists.osgeo.org/mailman/listinfo/gdal-dev
>


-- 
Alan Snow
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/gdal-dev/attachments/20220301/bef4dc4b/attachment.html>