Turn bounding_box class to function

[atteggiani]

The bounding_box class:

replace_landsurface/src/replace_landsurface_with_ERA5land_IC.py

Lines 37 to 115 in 2dcc1d7

	class bounding_box():
	""" Container class to hold spatial extent information."""
	def __init__(self, ncfname, maskfname, var):
	"""
	Initialization function for bounding_box class
	Parameters
	----------
	ncfname : POSIX string
	POSIX path to input data (from the NetCDF archive)
	maskfname : POSIX string
	POSIX path to mask information to define the data to be cut out
	var : string
	The name of the mask variable that defines the spatial extent
	Returns
	-------
	None. The variables describing the spatial extent are definied within bounding box object.
	"""
	# Read in the mask and get the minimum/maximum latitude and longitude information
	if Path(maskfname).exists():
	d = iris.load(maskfname, var)
	d = d[0]
	d = xr.DataArray.from_iris(d)
	lons = d['longitude'].data
	lonmin = np.min(lons)
	lonmax = np.max(lons)
	if lonmax > 180.:
	lonmax = lonmax-360.
	lats = d['latitude'].data
	latmin = np.min(lats)
	latmax = np.max(lats)
	d.close()
	else:
	print(f'ERROR: File {maskfname} not found', file=sys.stderr)
	raise
	# Read in the file from the high-res netcdf archive
	if Path(ncfname).exists():
	d = xr.open_dataset(ncfname)
	else:
	print(f'ERROR: File {ncfname} not found', file=sys.stderr)
	sys.exit(1)
	# Get the longitude information
	lons = d['longitude'].data
	# Coping with numerical inaccuracy
	adj = (lons[1] - lons[0])/2.
	# Work out which longitudes define the minimum/maximum extents of the grid of interest
	lonmin_index = np.argwhere((lons > lonmin - adj) & (lons < lonmin + adj))[0][0]
	lonmax_index = np.argwhere((lons > lonmax - adj) & (lons < lonmax + adj))[0][0]
	# Get the latitude information
	lats = d['latitude'].data
	# Work out which longitudes define the minimum/maximum extents of the grid of interest
	# use the same adjustment as for longitude
	latmin_index = np.argwhere((lats > latmin - adj) & (lats < latmin + adj))[0][0]
	latmax_index = np.argwhere((lats > latmax - adj) & (lats < latmax + adj))[0][0]
	# Swap the latitude min/max if upside down (is upside down for era5-land)
	if latmax_index < latmin_index:
	tmp_index=latmin_index
	latmin_index=latmax_index
	latmax_index=tmp_index
	# Set the boundaries
	self.lonmin=lonmin_index
	self.lonmax=lonmax_index
	self.latmin=latmin_index
	self.latmax=latmax_index

ultimately, is used to get the extents of the input file in terms of min and max longitudes and latitudes (lonmin, lonmax, latmin and latmax).

The class does not currently hold any methods (apart from the __init__) or attributes apart from the lat/lon extents it needs to produce.

I suggest turning it into a function to:

• simplify the code
• simplify testing

[truth-quark]

Indeed, the single method/__init__ only class is a sign of unclear design. Also, __init__ is doing multiple things: file I/O, lat/long calcs, control flow, array ops & setting member vars. It's a sign to break down the problem into smaller pieces.

Shapely may suit for providing generic geometry data structures too:
https://shapely.readthedocs.io/en/stable/manual.html