3535import iris
3636import iris .exceptions
3737import iris .unit
38- from iris .coord_systems import GeogCS
38+ from iris .coord_systems import GeogCS , RotatedGeogCS , TransverseMercator
3939from iris .fileformats .rules import is_regular , regular_step
4040from iris .fileformats .grib import grib_phenom_translation as gptx
4141from iris .fileformats .grib ._load_convert import (_STATISTIC_TYPE_NAMES ,
@@ -118,13 +118,16 @@ def ensure_set_int32_value(grib, key, value):
118118# Reference Regulation 92.1.6
119119_DEFAULT_DEGREES_UNITS = 1.0e-6
120120
121+ # Units scaling for Grid Definition Template 3.12
122+ _TRANSVERSE_MERCATOR_DEGREES_UNITS = 1.0e-2
121123
122124###############################################################################
123125#
124126# Identification Section 1
125127#
126128###############################################################################
127129
130+
128131def centre (cube , grib ):
129132 # TODO: read centre from cube
130133 gribapi .grib_set_long (grib , "centre" , 74 ) # UKMO
@@ -198,8 +201,7 @@ def shape_of_the_earth(cube, grib):
198201 # Spherical earth.
199202 if ellipsoid .inverse_flattening == 0.0 :
200203 gribapi .grib_set_long (grib , "shapeOfTheEarth" , 1 )
201- gribapi .grib_set_long (grib ,
202- "scaleFactorOfRadiusOfSphericalEarth" , 0 )
204+ gribapi .grib_set_long (grib , "scaleFactorOfRadiusOfSphericalEarth" , 0 )
203205 gribapi .grib_set_long (grib , "scaledValueOfRadiusOfSphericalEarth" ,
204206 ellipsoid .semi_major_axis )
205207 # Oblate spheroid earth.
@@ -256,7 +258,7 @@ def scanning_mode_flags(x_coord, y_coord, grib):
256258 int (y_coord .points [1 ] - y_coord .points [0 ] > 0 ))
257259
258260
259- def latlon_grid_common (cube , grib ):
261+ def horizontal_grid_common (cube , grib ):
260262 # Grib encoding of the sequences of X and Y points.
261263 y_coord = cube .coord (dimensions = [0 ])
262264 x_coord = cube .coord (dimensions = [1 ])
@@ -333,7 +335,7 @@ def grid_definition_template_0(cube, grib):
333335 """
334336 # Constant resolution, aka 'regular' true lat-lon grid.
335337 gribapi .grib_set_long (grib , "gridDefinitionTemplateNumber" , 0 )
336- latlon_grid_common (cube , grib )
338+ horizontal_grid_common (cube , grib )
337339 latlon_points_regular (cube , grib )
338340
339341
@@ -354,7 +356,7 @@ def grid_definition_template_1(cube, grib):
354356 rotated_pole (cube , grib )
355357
356358 # Encode the lat/lon points.
357- latlon_grid_common (cube , grib )
359+ horizontal_grid_common (cube , grib )
358360 latlon_points_regular (cube , grib )
359361
360362
@@ -372,16 +374,17 @@ def grid_definition_template_5(cube, grib):
372374 # Without this, setting "gridDefinitionTemplateNumber" = 5 causes an
373375 # immediate error.
374376 # See: https://software.ecmwf.int/issues/browse/SUP-1095
375- # This is acceptable, as the subsequent call to 'latlon_grid_common' will
376- # set these to the correct horizontal dimensions (by calling 'grid_dims').
377+ # This is acceptable, as the subsequent call to 'horizontal_grid_common'
378+ # will set these to the correct horizontal dimensions
379+ # (by calling 'grid_dims').
377380 gribapi .grib_set (grib , "Ni" , 1 )
378381 gribapi .grib_set (grib , "Nj" , 1 )
379382 gribapi .grib_set (grib , "gridDefinitionTemplateNumber" , 5 )
380383
381384 # Record details of the rotated coordinate system.
382385 rotated_pole (cube , grib )
383386 # Encode the lat/lon points.
384- latlon_grid_common (cube , grib )
387+ horizontal_grid_common (cube , grib )
385388 latlon_points_irregular (cube , grib )
386389
387390
@@ -402,30 +405,41 @@ def grid_definition_template_12(cube, grib):
402405
403406 # Set some keys specific to GDT12.
404407 # Encode the horizontal points.
405- M_TO_CM = 100
406408 x_step = regular_step (x_coord )
407409 y_step = regular_step (y_coord )
408- gribapi .grib_set (grib , "Di" , float (abs (x_step ))* M_TO_CM )
409- gribapi .grib_set (grib , "Dj" , float (abs (y_step ))* M_TO_CM )
410- latlon_grid_common (cube , grib )
410+ gribapi .grib_set (grib , "Di" ,
411+ float (abs (x_step )) / _TRANSVERSE_MERCATOR_DEGREES_UNITS )
412+ gribapi .grib_set (grib , "Dj" ,
413+ float (abs (y_step )) / _TRANSVERSE_MERCATOR_DEGREES_UNITS )
414+ horizontal_grid_common (cube , grib )
411415
412416 # GRIBAPI expects unsigned ints in X1, X2, Y1, Y2 but it should accept
413417 # signed ints, so work around this.
414418 # See https://software.ecmwf.int/issues/browse/SUP-1101
415- ensure_set_int32_value (grib , "Y1" , int (y_coord .points [0 ]* M_TO_CM ))
416- ensure_set_int32_value (grib , "Y2" , int (y_coord .points [- 1 ]* M_TO_CM ))
417- ensure_set_int32_value (grib , "X1" , int (x_coord .points [0 ]* M_TO_CM ))
418- ensure_set_int32_value (grib , "X2" , int (x_coord .points [- 1 ]* M_TO_CM ))
419+ ensure_set_int32_value (
420+ grib , "Y1" ,
421+ int (y_coord .points [0 ] / _TRANSVERSE_MERCATOR_DEGREES_UNITS ))
422+ ensure_set_int32_value (
423+ grib , "Y2" ,
424+ int (y_coord .points [- 1 ] / _TRANSVERSE_MERCATOR_DEGREES_UNITS ))
425+ ensure_set_int32_value (
426+ grib , "X1" ,
427+ int (x_coord .points [0 ] / _TRANSVERSE_MERCATOR_DEGREES_UNITS ))
428+ ensure_set_int32_value (
429+ grib , "X2" ,
430+ int (x_coord .points [- 1 ] / _TRANSVERSE_MERCATOR_DEGREES_UNITS ))
419431
420432 # Lat and lon of reference point are measured in millionths of a degree.
421433 gribapi .grib_set (grib , "latitudeOfReferencePoint" ,
422- cs .latitude_of_projection_origin * 1e6 )
434+ cs .latitude_of_projection_origin / _DEFAULT_DEGREES_UNITS )
423435 gribapi .grib_set (grib , "longitudeOfReferencePoint" ,
424- cs .longitude_of_central_meridian * 1e6 )
436+ cs .longitude_of_central_meridian / _DEFAULT_DEGREES_UNITS )
425437
426438 # False easting and false northing are measured in units of (10^-2)m.
427- gribapi .grib_set (grib , "XR" , cs .false_easting * M_TO_CM )
428- gribapi .grib_set (grib , "YR" , cs .false_northing * M_TO_CM )
439+ gribapi .grib_set (grib , "XR" ,
440+ cs .false_easting / _TRANSVERSE_MERCATOR_DEGREES_UNITS )
441+ gribapi .grib_set (grib , "YR" ,
442+ cs .false_northing / _TRANSVERSE_MERCATOR_DEGREES_UNITS )
429443
430444 # GRIBAPI expects a signed int for scaleFactorAtReferencePoint
431445 # but it should accept a float, so work around this.
@@ -449,23 +463,23 @@ def grid_definition_section(cube, grib):
449463 cs = x_coord .coord_system # N.B. already checked same cs for x and y.
450464 regular_x_and_y = is_regular (x_coord ) and is_regular (y_coord )
451465
452- if isinstance (cs , iris . coord_systems . GeogCS ):
466+ if isinstance (cs , GeogCS ):
453467 if not regular_x_and_y :
454468 raise iris .exceptions .TranslationError (
455469 'Saving an irregular latlon grid to GRIB (PDT3.4) is not '
456470 'yet supported.' )
457471
458472 grid_definition_template_0 (cube , grib )
459473
460- elif isinstance (cs , iris . coord_systems . RotatedGeogCS ):
474+ elif isinstance (cs , RotatedGeogCS ):
461475 # Rotated coordinate system cases.
462476 # Choose between GDT 3.1 and 3.5 according to coordinate regularity.
463477 if regular_x_and_y :
464478 grid_definition_template_1 (cube , grib )
465479 else :
466480 grid_definition_template_5 (cube , grib )
467481
468- elif isinstance (cs , iris . coord_systems . TransverseMercator ):
482+ elif isinstance (cs , TransverseMercator ):
469483 # Transverse Mercator coordinate system (template 3.12).
470484 grid_definition_template_12 (cube , grib )
471485
0 commit comments