cf-python icon indicating copy to clipboard operation
cf-python copied to clipboard

Field.subspace fails to work properly for cyclic field

Open sws98slg opened this issue 6 months ago • 4 comments

On upgrading from a very old version of cfpython (3.7.0) to the latest version (3.18.0) field.subspace no longer seems to work properly for a cyclic field. As an example using the datasets provided with cfplot

f=cf.read('cfplot_data/tas_A1.nc')[0] f.subspace(X=cf.wi(175,-175)) Gives for the old version: Out[139]: <CF Field: air_temperature(time(1680), latitude(73), longitude(3)) K> (correct answer) But for the new version: ValueError: No indices found from: X=<CF Query: (wi [175, -175])>

f.subspace(X=cf.wi(-5,5)) Gives for both versions Out[140]: <CF Field: air_temperature(time(1680), latitude(73), longitude(3)) K> (correct answer)

f.subspace(X=cf.wi(355,5)) Gives for the old version: Out[142]: <CF Field: air_temperature(time(1680), latitude(73), longitude(3)) K> But for the new version: ValueError: No indices found from: X=<CF Query: (wi [355, 5])>

I've tried explicitly setting the field to be cyclic

f.cyclic('X',iscyclic=True) Before doing the subspacing but this doesn't help.

I hope I've explained this fully. I can see there was a previous issue reported (#828) that seems related, but this is marked as closed. For further info my environments are below.

My old environment is Platform: Linux-3.10.0-1160.119.1.el7.x86_64-x86_64-with-centos-7.9.2009-Core HDF5 library: 1.10.4 netcdf library: 4.6.3 udunits2 library: libudunits2.so.0 python: 3.7.7 netCDF4: 1.5.3 cftime: 1.2.1 numpy: 1.19.2 psutil: 5.6.3 scipy: 1.4.1 matplotlib: 3.3.2 ESMF: 7.1.0r cfdm: 1.8.7.0 cfunits: 3.3.0 cfplot: 3.0.38 cf: 3.7.0

And my new one is Platform: Linux-4.18.0-553.34.1.el8_10.x86_64-x86_64-with-glibc2.28 Python: 3.12.11 packaging: 25.0 numpy: 2.3.2 cfdm.core: 1.12.2.0 udunits2 library: /home/users/sws98slg/.conda/envs/myenv3_RACC2/lib/libudunits2.so.0 HDF5 library: 1.14.6 netcdf library: 4.9.2 netCDF4: 1.7.2 h5netcdf: 1.6.4 h5py: 3.14.0 zarr: 3.1.1 s3fs: 2025.7.0 scipy: 1.16.0 dask: 2025.7.0 cftime: 1.6.4 cfunits: 3.3.7 cfdm: 1.12.2.0 esmpy/ESMF: 8.8.1 psutil: 7.0.0 matplotlib: 3.10.5 cfplot: 3.4.0 cf: 3.18.0

sws98slg avatar Aug 11 '25 17:08 sws98slg

Hi, thanks for the report.

The answer is to order the arguments cf.wi so that they are increasing:

>>> import cf
>>> f = cf.read('tas_A1.nc')[0]
>>> f.subspace(X=cf.wi(-175, 175))
<CF Field: air_temperature(time(1680), latitude(73), longitude(93)) K>
>>> f.subspace(X=cf.wi(-5, 5))
<CF Field: air_temperature(time(1680), latitude(73), longitude(3)) K>
>>> f.subspace(X=cf.wi(355, 365))
<CF Field: air_temperature(time(1680), latitude(73), longitude(3)) K>

OK - that works, but I agree that it should be clever enough, in the cyclic case, to map cf.wi(355, 5) to cf.wi(-5, 5)) for you. I don't know why it stopped working sometime after v3.7.0 (2020) - that's far too long ago for me to remember or want to work out why :) - but I'll put in a fix so your original code works as you expected, for the next release (not the one that's coming later this week, but the one due on early October).

davidhassell avatar Aug 12 '25 09:08 davidhassell

I have another test case (based on a report on the JASMIN helpdesk), which seems to be somewhat related to this issue, but the exact solution/workaround proposed above does not seem to hold in this case because the longitudes used in this case (-21, -19) are already in increasing order.

import cf
field_data = cf.read("test3.nc")[0]
field_data.dim('dimensioncoordinate0').set_property('standard_name', 'latitude')
field_data.dim('dimensioncoordinate1').set_property('standard_name', 'longitude')
field_data.cyclic('X', period=360)
sub = field_data.subspace(latitude=cf.wi(34,36), longitude=cf.wi(-21,-19))

This uses the input file that is wrapped in this zip file (because GitHub doesn't like .nc attachments). In fact it is subsetted from an original 4d file, as you will see in the history attribute, for sake of a more minimal example.

It fails with cf-python 3.18.0 but succeeds with cf-python 3.16.2. (These are as installed in different versions of Jaspy.)

In the failure case, it will nonetheless work for the longitude range (360-21, 360-19).

alaniwi avatar Aug 15 '25 11:08 alaniwi

For convenience, here is a dump of the test file (with edits to the list of coordinate values for conciseness):

$ ncdump -c test3.nc 
netcdf test3 {
dimensions:
	latitude = 180 ;
	longitude = 360 ;
variables:
	float latitude(latitude) ;
		latitude:units = "degrees_north" ;
		latitude:long_name = "latitude" ;
	float longitude(longitude) ;
		longitude:units = "degrees_east" ;
		longitude:long_name = "longitude" ;
	double msl(latitude, longitude) ;
		msl:_FillValue = -32767. ;
		msl:missing_value = -32767s ;
		msl:units = "Pa" ;
		msl:long_name = "Mean sea level pressure" ;
		msl:standard_name = "air_pressure_at_mean_sea_level" ;
		msl:cell_methods = "time, number: mean" ;

// global attributes:
		:Conventions = "CF-1.6" ;
		:history = "Fri Aug 15 12:11:26 2025: ncks -x -v number,time /tmp/test2.nc /tmp/test3.nc\n",
			"Fri Aug 15 12:10:55 2025: ncwa -O -a time,number /tmp/test.nc /tmp/test2.nc\n",
			"Fri Aug 15 12:10:14 2025: ncks -O -d time,0,0 -d number,0,0 /gws/nopw/j04/ncas_climate_vol1/users/lbaker/c3s_data/reconverted_data/ukmo_mslp_members_oct_sys602_hindcast.nc /tmp/test.nc\n",
			"2025-07-07 13:26:49 GMT by grib_to_netcdf-2.36.0: grib_to_netcdf ukmo_mslp_members_oct_sys602_hindcast.grib -o reconverted_data/ukmo_mslp_members_oct_sys602_hindcast.nc" ;
		:NCO = "netCDF Operators version 5.3.3 (Homepage = http://nco.sf.net, Code = http://github.com/nco/nco, Citation = 10.1016/j.envsoft.2008.03.004)" ;
data:

 latitude = 89.5, 88.5, 87.5, [...etc...] -88.5, -89.5 ;

 longitude = 0.5, 1.5, 2.5, [...etc...], 357.5, 358.5, 359.5 ;
}

alaniwi avatar Aug 15 '25 12:08 alaniwi

Hi @alaniwi,

Thanks for this - another, related bug is causing the behaviour you are seeing. It's related to the fact that the X coordinates do not have bounds (which is fine, but it led to some dodgy logic in cf.Field.cyclic).

PR to fix both issues to follow ...

David

(P.S. you don't need to set the standard names in this case - the units are sufficient to define 'X' and 'Y': https://ncas-cms.github.io/cf-python/attribute/cf.DimensionCoordinate.X.html)

davidhassell avatar Oct 15 '25 11:10 davidhassell