Conversion of atmosphere hybrid sigma pressure coordiate

[lguo-uk]

I try to cf.read a .nc file, in which there are a variabe BC(time,lev,lat,lon), lev is atmosphere hybrid sigma pressure coordinate, and associated variables for converting vertical coordiante.

When I ncdump the .nc file, BC and parameters are stored separately (see attchment 1). However, once cf.read, all these parameters become domain ancils to BC of cf.field (see attachment 2).

My questions:

1. Is there any function to automatically convert vertical coodinate to pressure coodinate?
2. If not, how can I extract parameters, so that I can convert it myself? I've tried BC.get_domain_ancillary("longname:vertical coordinate formula term: ap(k)"). However, I've told that 'Field' object has no attribute 'get_domainancillary'.

Reason for asking this question is to plot the vertical profile of the variable.

Many Thanks, Liang

Attachment 1. ncdump -h N2000_AEROSLO_ensmean.2050-2099_BC.nc

netcdf N2000_AEROSLO_ensmean.2050-2099_BC { dimensions: time = UNLIMITED ; // (600 currently) lev = 26 ; lat = 96 ; lon = 144 ; bnds = 2 ; ilev = 27 ; variables: float BC(time, lev, lat, lon) ; BC:long_name = "Concentration" ; BC:units = "kg m-3" ; BC:cell_methods = "time: mean" ; double P0 ; P0:long_name = "reference pressure" ; P0:units = "Pa" ; float PS(time, lat, lon) ; PS:long_name = "Surface pressure" ; PS:units = "Pa" ; PS:cell_methods = "time: mean" ; double a(lev) ; a:long_name = "vertical coordinate formula term: ap(k)" ; a:units = "Pa" ; double a_bnds(lev, bnds) ; a_bnds:long_name = "vertical coordinate formula term: ap(k+1/2)" ; a_bnds:units = "Pa" ; double b(lev) ; b:long_name = "vertical coordinate formula term: b(k)" ; b:units = "1" ; double b_bnds(lev, bnds) ; b_bnds:long_name = "vertical coordinate formula term: b(k+1/2)" ; b_bnds:units = "1" ; double ilev(ilev) ; ilev:standard_name = "atmosphere_hybrid_sigma_pressure_coordinate" ; ilev:axis = "Z" ; ilev:positive = "down" ; ilev:long_name = "hybrid sigma pressure coordinate" ; ilev:units = "level" ; ilev:formula_terms = "ap: ap b: b ps: ps" ; ilev:bounds = "ilev_bnds" ; double ilev_bnds(ilev, bnds) ; ilev_bnds:standard_name = "atmosphere_hybrid_sigma_pressure_coordinate" ; ilev_bnds:units = "level" ; ilev_bnds:formula_terms = "ap: ap_bnds b: b_bnds ps: ps" ; double lat(lat) ; lat:standard_name = "latitude" ; lat:long_name = "latitude" ; lat:units = "degrees_north" ; lat:axis = "Y" ; double lev(lev) ; lev:standard_name = "atmosphere_hybrid_sigma_pressure_coordinate" ; lev:axis = "Z" ; lev:positive = "down" ; lev:long_name = "hybrid sigma pressure coordinate" ; lev:units = "level" ; lev:formula_terms = "a: a b: b p0: P0 ps: PS" ; lev:bounds = "lev_bnds" ; double lev_bnds(lev, bnds) ; lev_bnds:standard_name = "atmosphere_hybrid_sigma_pressure_coordinate" ; lev_bnds:units = "level" ; lev_bnds:formula_terms = "a: a_bnds b: b_bnds p0: P0 ps: PS" ; double lon(lon) ; lon:standard_name = "longitude" ; lon:long_name = "longitude" ; lon:units = "degrees_east" ; lon:axis = "X" ; double time(time) ; time:standard_name = "time" ; time:long_name = "time" ; time:bounds = "time_bnds" ; time:units = "days since 2000-01-01 00:00:00" ; time:calendar = "365_day" ; time:axis = "T" ; double time_bnds(time, bnds) ;

Attachment 2. print(BC)

Field: long_name=Concentration (ncvar%BC) Data                  : long_name=Concentration(time(600), atmosphere_hybrid_sigma_pressure_coordinate(26), latitude(96), longitude(144)) kg m-3 Cell methods                  : time(600): mean Dimension coords                  : time(600) = [2050-02-01 00:00:00, ..., 2100-01-01 00:00:00] 365_day                  : atmosphere_hybrid_sigma_pressure_coordinate(26) = [3.5446380000000097, ..., 992.5560999999998] level                  : latitude(96) = [-90.0, ..., 90.0] degrees_north                  : longitude(144) = [0.0, ..., 357.5] degrees_east Coord references                  : standard_name:atmosphere_hybrid_sigma_pressure_coordinate Domain ancils
                 : long_name=vertical coordinate formula term: ap(k)(atmosphere_hybrid_sigma_pressure_coordinate(26)) = [0.0035446380000000097, ..., 0.0] Pa                  : long_name=vertical coordinate formula term: b(k)(atmosphere_hybrid_sigma_pressure_coordinate(26)) = [0.0, ..., 0.9925560999999998] 1                  : long_name=reference pressure() = 100000.0 Pa                  : long_name=Surface pressure(time(600), latitude(96), longitude(144)) = [[[68659.515625, ..., 101580.4921875]]] Pa

[sadielbartholomew]

Hi @lguo-uk, thanks for your question. I'm just quickly commenting to inform you that the cf-python maintainers, David Hassell and myself, are on leave (for at least most of) this week, but we will will strive to answer your question above upon our return. I might get an opportunity to provide a partial or full answer tomorrow or Friday. Thanks for your patience and apologies for the wait.

[davidhassell]

Hi @lguo-uk, there is not yet a builtin function to do this, although it is on out list of items we'd like to add. It's not been prioritised as no one has asked for it until now, so we'll review things and see if we can give it a higher priority.

Until then, you can indeed do it yourself by retrieving the required domain ancillaries and coordinates and applying the appropriate formula. You nearly had the right command. Try:

# Get the domain ancillaries
ak = BC.domain_ancillary("long_name:vertical coordinate formula term: ap(k)")
bk = BC.domain_ancillary("long_name=vertical coordinate formula term: b(k)")
p0 = BC.domain_ancillary("long_name=reference pressure()")
ps = BC.domain_ancillary("long_name=Surface pressure")

# Add a size 1 Z dimension to the surface pressure so that broadcasting works
ps = ps.insert_dimension(-1)

# Calculate air pressure
p = ak * p0 + bk * ps

# Create the new coordinate object
c = cf.AuxiliaryCoordinate()
c.standard_name = 'air_pressure'
c.set_data(p.data)

# Insert the new coordinates into the field
BC.set_construct(c)

I haven't tested this, not having the file, so I hope that I've not made any typos. Anyway, the idea should be clear. Do let us know if this works, or otherwise.

All the best, David

[lguo-uk]

Hi @davidhassell @sadielbartholomew ,

Thanks for the reply.

I've given it a try. However, performing calculation of p using DomainAncillary seems problematic.

(Pdb) p = ak p0+bk ps SyntaxError: invalid syntax (Pdb) p = ak p0 * SyntaxError: invalid syntax (Pdb) p = ak p0 SyntaxError: invalid syntax (Pdb) p = ak * SyntaxError: invalid syntax

Any idea?

Liang

[davidhassell]

Hi, Not sure why that doesn't work, here is an example that does - you should be able to replicate it:

import cf

f = cf.example_field(1)

a = f.domain_ancillary('ncvar%a')
b = f.domain_ancillary('ncvar%b')
orog = f.domain_ancillary('surface_altitude')

orog = orog.insert_dimension(-1)

z = a + b * orog

c = cf.AuxiliaryCoordinate()
c.standard_name = 'altitude'
c.set_data(z.data)

print(repr(c))

Produces:

<CF AuxiliaryCoordinate: altitude(10, 9, 1) m>

[lguo-uk]

Thanks, @davidhassell, for the example. It works fine with my JASMIN account.

One possible reason for causing the error might be that units of ap is wrong. It should be 1 instead of Pa. After overriding ap's units, operation runs successfully.

BC=cf.read('N2000_5xBCemEA_ensmean.2050-2099_BC_original.nc').select_field("long_name=Concentration")

# Get the domain ancillaries
ak=BC.domain_ancillary("long_name=vertical coordinate formula term: ap(k)")
ak.override_units("1",inplace=True) # original units is "Pa". However, ak shows values of scale, needs to be multiplied with p0 to form ap.
bk=BC.domain_ancillary("long_name=vertical coordinate formula term: b(k)")
p0=BC.domain_ancillary("long_name=reference pressure")
ps=BC.domain_ancillary("long_name=Surface pressure")

# Add a size 1 Z dimension to the surface pressure so that broadcasting works
# Z dimention has to be appended on the back for calculation
ps=ps.insert_dimension(-1)

# Calculate air pressure
p=ak*p0+bk*ps
p.transpose([0,3,2,1],inplace=True)

# Create the new coordinate object
c=cf.AuxiliaryCoordinate()
c.standard_name = 'air_pressure'
c.set_data(p.data)

# Insert the new coordinates into the field
BC.set_construct(c)

(Pdb) print(BC)
Field: long_name=Concentration (ncvar%BC)
-----------------------------------------
Data            : long_name=Concentration(time(600), atmosphere_hybrid_sigma_pressure_coordinate(26), latitude(96), longitude(144)) kg m-3
Cell methods    : time(600): mean
Dimension coords: time(600) = [2050-02-01 00:00:00, ..., 2100-01-01 00:00:00] 365_day
                : atmosphere_hybrid_sigma_pressure_coordinate(26) = [3.5446380000000097, ..., 992.5560999999998] level
                : latitude(96) = [-90.0, ..., 90.0] degrees_north
                : longitude(144) = [0.0, ..., 357.5] degrees_east
Auxiliary coords: air_pressure(time(600), atmosphere_hybrid_sigma_pressure_coordinate(26), longitude(144), latitude(96)) = [[[[354.46380000000096, ..., 99975.9110635078]]]] Pa
Coord references: standard_name:atmosphere_hybrid_sigma_pressure_coordinate
Domain ancils   : long_name=vertical coordinate formula term: ap(k)(atmosphere_hybrid_sigma_pressure_coordinate(26)) = [0.0035446380000000097, ..., 0.0] 1
                : long_name=vertical coordinate formula term: b(k)(atmosphere_hybrid_sigma_pressure_coordinate(26)) = [0.0, ..., 0.9925560999999998] 1
                : long_name=reference pressure() = 100000.0 Pa
                : long_name=Surface pressure(time(600), latitude(96), longitude(144)) = [[[68686.3203125, ..., 100725.703125]]] Pa

My motivation of converting the vertical coordiante from hybrid sigma from pressure is to visualise the variable (black carbon concentration) in a concentration-pressure profile plot. Therefore, I've calculated an area mean using BC_mean=BC.subspace(latitude=cf.wi(10,30),longitude=cf.wi(100,120)).collapse("mean","area"). However, after the mean, the ancillary coordinates air_pressure is gone:

(Pdb) print(BC_mean)
Field: long_name=Concentration (ncvar%BC)
-----------------------------------------
Data            : long_name=Concentration(time(600), atmosphere_hybrid_sigma_pressure_coordinate(26), latitude(1), longitude(1)) kg m-3
Cell methods    : time(600): mean area: mean
Dimension coords: time(600) = [2050-02-01 00:00:00, ..., 2100-01-01 00:00:00] 365_day
                : atmosphere_hybrid_sigma_pressure_coordinate(26) = [3.5446380000000097, ..., 992.5560999999998] level
                : latitude(1) = [19.89473684210526] degrees_north
                : longitude(1) = [110.0] degrees_east
Coord references: standard_name:atmosphere_hybrid_sigma_pressure_coordinate
Domain ancils   : long_name=vertical coordinate formula term: ap(k)(atmosphere_hybrid_sigma_pressure_coordinate(26)) = [0.0035446380000000097, ..., 0.0] 1
                : long_name=vertical coordinate formula term: b(k)(atmosphere_hybrid_sigma_pressure_coordinate(26)) = [0.0, ..., 0.9925560999999998] 1
                : long_name=reference pressure() = 100000.0 Pa

Here are questions:

1. How do I turn ancillary coordinates air_pressure from cf.domainancillary.DomainAncillary to cf.Filed? So that, .subspace and .collapse could work on it.
2. Or, will it be possble to keep an area meam ancillary coordinates air_pressure with BC_mean? So that, the information can be used futher analysis or visualisation.

Cheers, Liang

[davidhassell]

Hi Liang,

Great - sounds like we're making progress (as I am on #142).

The auxiliary coordinates are removed in general because I didn't think that it was clear how to process them to meaningful values. It may be in this case that an area average of the coordinates is appropriate, but in general this is not obvious.

However, I'm thinking now of an extension to collapse that allows you to request a collapse (using "mean") of such coordinates so that they may be retained.

But in the mean time, the easiest way forward is probably to do as you suggest - turn the altitudes into a field in their own right. This is easy:

altitude = BC.convert('altitude')

E.g.

In [37]: import cf

In [38]: f = cf.example_field(1)

In [39]: print(f)
Field: air_temperature (ncvar%ta)
---------------------------------
Data            : air_temperature(atmosphere_hybrid_height_coordinate(1), grid_latitude(10), grid_longitude(9)) K
Cell methods    : grid_latitude(10): grid_longitude(9): mean where land (interval: 0.1 degrees) time(1): maximum
Field ancils    : air_temperature standard_error(grid_latitude(10), grid_longitude(9)) = [[0.76, ..., 0.32]] K
Dimension coords: atmosphere_hybrid_height_coordinate(1) = [1.5]
                : grid_latitude(10) = [2.2, ..., -1.76] degrees
                : grid_longitude(9) = [-4.7, ..., -1.18] degrees
                : time(1) = [2019-01-01 00:00:00]
Auxiliary coords: latitude(grid_latitude(10), grid_longitude(9)) = [[53.941, ..., 50.225]] degrees_N
                : longitude(grid_longitude(9), grid_latitude(10)) = [[2.004, ..., 8.156]] degrees_E
                : long_name=Grid latitude name(grid_latitude(10)) = [--, ..., b'kappa']
Cell measures   : measure:area(grid_longitude(9), grid_latitude(10)) = [[2391.9657, ..., 2392.6009]] km2
Coord references: grid_mapping_name:rotated_latitude_longitude
                : standard_name:atmosphere_hybrid_height_coordinate
Domain ancils   : ncvar%a(atmosphere_hybrid_height_coordinate(1)) = [10.0] m
                : ncvar%b(atmosphere_hybrid_height_coordinate(1)) = [20.0]
                : surface_altitude(grid_latitude(10), grid_longitude(9)) = [[0.0, ..., 270.0]] m

In [40]: orog = f.convert('surface_altitude')

In [41]: print(orog)
Field: surface_altitude (ncvar%surface_altitude)
------------------------------------------------
Data            : surface_altitude(grid_latitude(10), grid_longitude(9)) m
Dimension coords: grid_latitude(10) = [2.2, ..., -1.76] degrees
                : grid_longitude(9) = [-4.7, ..., -1.18] degrees
Auxiliary coords: latitude(grid_latitude(10), grid_longitude(9)) = [[53.941, ..., 50.225]] degrees_N
                : longitude(grid_longitude(9), grid_latitude(10)) = [[2.004, ..., 8.156]] degrees_E
                : long_name=Grid latitude name(grid_latitude(10)) = [--, ..., b'kappa']
Cell measures   : measure:area(grid_longitude(9), grid_latitude(10)) = [[2391.9657, ..., 2392.6009]] km2
Coord references: grid_mapping_name:rotated_latitude_longitude

[lguo-uk]

Great, @davidhassell. cf.Field.convert() is handy 👍. Thank you for sharing that.

I've successfully coverted the vertical coordinate from Hybrid Sigma to Pressure 😃. Looking forward to the new functionality #142.

[davidhassell]

Thanks, @lguo-uk, that's good news. The new functionality won't be long - I've been working on it this week.