radar_tools API

domutils.radar_tools.get_instantaneous.get_instantaneous(valid_date=None, data_path=None, data_recipe=None, desired_quantity=None, median_filt=None, coef_a=None, coef_b=None, qced=True, missing=-9999.0, latlon=False, dest_lon=None, dest_lat=None, average=False, nearest_time=None, smooth_radius=None, odim_latlon_file=None, verbose=0)

Get instantaneous precipitation from various sources

Provides one interface for:

  • reading the following input formats:

    • Odim h5 composites

    • Stage IV files

    • MRMS data

    • CMC “standard files”

  • output to an arbitrary output grid

  • Consistent filtering (median and/or circular boxcar) on precip observations and the accompanying quality index

  • Various types of averaging

  • Find the nearest time where observations are available

The file extension present in data_recipe determines the file type of the source data Files having no extension are assumed to be in the ‘standard’ file format

Parameters:

  • valid_date (Optional[Any]) – datetime object with the validity date of the precip field

  • data_path (Optional[str]) – path to directory where data is expected

  • data_recipe (Optional[str]) – datetime code for constructing the file name eg: /%Y/%m/%d/qcomp_%Y%m%d%H%M.h5’ the filename will be obtained with data_path + valid_date.strftime(data_recipe)

  • desired_quantity (Optional[str]) – What quantity is desired in output dictionary precip_rate in [mm/h] and reflectivity in [dBZ] are supported

  • median_filt (Optional[int]) – If specified, a median filter will be applied on the data being read and the associated quality index. eg. medialFilter=3 will apply a median filter over a 3x3 boxcar If unspecified, no filtering is applied

  • coef_a (Optional[float]) – Coefficient a in Z = aR^b

  • coef_b (Optional[float]) – Coefficient b in Z = aR^b

  • qced (Optional[bool]) – Only for Odim H5 composites When True (default), Quality Controlled reflectivity (DBZH) will be returned. When False, raw reflectivity (TH) will be returned.

  • missing (Optional[float]) – Value that will be assigned to missing data

  • latlon (Optional[bool]) – Return latitudes and longitudes grid of the data

  • dest_lon (Optional[Any]) – Longitudes of destination grid. If not provided data is returned on its original grid

  • dest_lat (Optional[Any]) – Latitudes of destination grid. If not provided data is returned on its original grid

  • average (Optional[bool]) – Use the averaging method to interpolate data (see geo_tools documentation), this can be slow

  • nearest_time (Optional[float]) – If set, rewind time until a match is found to an integer number of nearestTime For example, with nearestTime=10, time will be rewinded to the nearest integer of 10 minutes

  • smooth_radius (Optional[float]) – Use the smoothing radius method to interpolate data, faster (see geo_tools documentation)

  • odim_latlon_file (Optional[str]) – file containing the latitude and longitudes of Baltrad mosaics in Odim H5 format

  • verbose (Optional[int]) – – Deprecated – Set >=1 to print info on execution steps

Returns:

None If no file matching the desired time is found

or

{

‘reflectivity’ 2D reflectivity on destination grid (if requested)

’precip_rate’ 2D precipitation rate on destination grid (if requested)

’total_quality_index’ Quality index of data with 1 = best and 0 = worse

’valid_date’ Actual validity date of data

’latitudes’ If latlon=True

’longitudes’ If latlon=True

}

domutils.radar_tools.get_accumulation.get_accumulation(end_date=None, duration=None, input_dt=None, desired_quantity=None, data_path=None, data_recipe=None, median_filt=None, coef_a=None, coef_b=None, qced=True, missing=-9999.0, latlon=False, dest_lon=None, dest_lat=None, average=False, nearest=None, smooth_radius=None, odim_latlon_file=None, verbose=0)

Get accumulated precipitation from instantaneous observations

This is essentially a wrapper around get_instantaneous. Data is read during the accumulation period and accumulated (in linear units of precipitation rates) taking the quality index into account. If interpolation to a different grid is desired, it is performed after the accumulation procedure.

If the desired quantity reflectivity or precip_rate is desired, then the returned quantity will reflect the average precipitation rate during the accumulation period.

With an endTime set to 16:00h and duration to 60 (minutes), data from:

  • 15:10h, 15:20h, 15:30h, 15:40h, 15:50h and 16:00h

will be accumulated

Parameters:

  • end_date (Optional[Any]) – datetime object representing the time (inclusively) at the end of the accumulation period

  • duration (Optional[Any]) – amount of time (minutes) during which precipitation should be accumulated

  • input_dt (Optional[int]) – time separation (minutes) of input data used for constructing accumulations

  • data_path (Optional[str]) – path to directory where data is expected

  • data_recipe (Optional[str]) – datetime code for constructing the file name eg: /%Y/%m/%d/qcomp_%Y%m%d%H%M.h5’ the filename will be obtained with data_path + valid_date.strftime(data_recipe)

  • desired_quantity (Optional[str]) –

    What quantity is desired in output dictionary. Supported values are:

    • accumulation Default option, the amount of water (in mm)

    • avg_precip_rate For average precipitation rate (in mm/h) during the accumulation period

    • reflectivity For the reflectivity (in dBZ) associated with the average precipitation rate

  • median_filt (Optional[int]) – If specified, a median filter will be applied on the data being read and the associated quality index. eg. medialFilter=3 will apply a median filter over a 3x3 boxcar If unspecified, no filtering is applied

  • coef_a (Optional[float]) – Coefficient a in Z = aR^b

  • coef_b (Optional[float]) – Coefficient b in Z = aR^b

  • qced (Optional[bool]) – Only for Odim H5 composites When True (default), Quality Controlled reflectivity (DBZH) will be returned. When False, raw reflectivity (TH) will be returned.

  • missing (Optional[float]) – Value that will be assigned to missing data

  • latlon (Optional[bool]) – Return latitudes and longitudes grid of the data

  • dest_lon (Optional[Any]) – Longitudes of destination grid. If not provided data is returned on its original grid

  • dest_lat (Optional[Any]) – Latitudes of destination grid. If not provided data is returned on its original grid

  • average (Optional[bool]) – Use the averaging method to interpolate data (see geo_tools documentation), this can be slow

  • nearest (Optional[float]) – If set, rewind time until a match is found to an integer number of nearest For example, with nearest=10, time will be rewinded to the nearest integer of 10 minutes

  • smooth_radius (Optional[float]) – Use the smoothing radius method to interpolate data, faster (see geo_tools documentation)

  • odim_latlon_file (Optional[str]) – file containing the latitude and longitudes of Baltrad mosaics in Odim H5 format

  • verbose (Optional[int]) – – Deprecated – Set >=1 to print info on execution steps

Returns:

None If no file matching the desired time is found

or

{

‘end_date’ End time for accumulation period

’duration’ Accumulation length (minutes)

’accumulation’ 2D accumulation (mm) on destination grid

’reflectivity’ 2D reflectivity on destination grid (if requested)

’precip_rate’ precipitation rate on destination grid (if requested)

’total_quality_index’ Quality index of data with 1 = best and 0 = worse

’latitudes’ If latlon=True

’longitudes’ If latlon=True

}

domutils.radar_tools.obs_process.obs_process(args=None)

Batch processing of precipitation observation

This function provides time and space interpolation of precipitiaton data. If desired, processing such as median filtering and boxcar-type filtering can be applied at the same time.

See for an example.

It is intended to be used as a command-line program, that would be called like:

python -m domutils.radar_tools.obs_process
          --input_data_dir   /space/hall4/sitestore/eccc/mrd/rpndat/dja001/data/radar_h5_composites/v8/
          --output_dir       /home/dja001/python/obs_process/outdir/
          --figure_dir       /home/dja001/python/obs_process/figdir/
          --fst_file_struc   %Y/%m/%d/%Y%m%d%H%M_mosaic.fst
          --input_file_struc %Y/%m/%d/qcomp_%Y%m%d%H%M.h5
          --h5_latlon_file   /home/dja001/shared_stuff/files/radar_continental_2.5km_2882x2032.pickle
          --t0               ${t_start}
          --tf               ${t_stop}
          --input_dt         10
          --sample_pr_file   /space/hall4/sitestore/eccc/mrd/rpndat/dja001/domains/hrdps_5p1_prp0.fst
          --ncores           40
          --complete_dataset True
          --median_filt      3
          --smooth_radius    4

Alternatively, is is possible to call this function directly from Python by defining a simple object whose attributes are the arguments. Such use is demonstrated in Interpolate a batch of observations.

Argument description:

read radar H5 files, interpolate/smooth and write to FST

usage: obs_process [-h] --input_data_dir INPUT_DATA_DIR --output_dir
                   OUTPUT_DIR --input_t0 INPUT_T0 --processed_file_struc
                   PROCESSED_FILE_STRUC --input_file_struc INPUT_FILE_STRUC
                   --input_dt INPUT_DT
                   [--tinterpolated_file_struc TINTERPOLATED_FILE_STRUC]
                   [--h5_latlon_file H5_LATLON_FILE] [--input_tf INPUT_TF]
                   [--fcst_len FCST_LEN] [--accum_len ACCUM_LEN]
                   [--output_t0 OUTPUT_T0] [--output_tf OUTPUT_TF]
                   [--output_dt OUTPUT_DT] [--t_interp_method T_INTERP_METHOD]
                   [--sample_pr_file SAMPLE_PR_FILE]
                   [--output_file_format OUTPUT_FILE_FORMAT] [--ncores NCORES]
                   [--complete_dataset COMPLETE_DATASET]
                   [--median_filt MEDIAN_FILT] [--smooth_radius SMOOTH_RADIUS]
                   [--figure_dir FIGURE_DIR] [--cartopy_dir CARTOPY_DIR]
                   [--figure_format FIGURE_FORMAT] [--log_level LOG_LEVEL]

Required named arguments

--input_data_dir

path of source radar mosaics files

--output_dir

directory for output fst files

--input_t0

yyyymmsshhmmss begining time; datestring

--processed_file_struc

strftime syntax for constructing fst filenames for output of obsprocess

--input_file_struc

strftime syntax for constructing H5 filenames

--input_dt

interval (minutes) between input radar mosaics

Optional named arguments

--tinterpolated_file_struc

strftime syntax for constructing time interpolated filenames

Default: “None”

--h5_latlon_file

Pickle file containing the lat/lons of the Baltrad grid

Default: “None”

--input_tf

yyyymmsshhmmss end time; datestring

Default: “None”

--fcst_len

duration of forecast (hours)

Default: None

--accum_len

duration of accumulation (minutes)

Default: “None”

--output_t0

yyyymmsshhmmss begining time; datestring

Default: “None”

--output_tf

yyyymmsshhmmss end time; datestring

Default: “None”

--output_dt

interval (minutes) between output radar mosaics

Default: “None”

--t_interp_method

time interpolation method

Default: “None”

--sample_pr_file

File containing PR to establish the domain

Default: “None”

--output_file_format

File format of processed files

Default: “npz”

--ncores

number of cores for parallel execution

Default: 1

--complete_dataset

Skip existing files, default is to clobber them

Default: “False”

--median_filt

box size (pixels) for median filter

Default: “None”

--smooth_radius

radius (km) where radar data be smoothed

Default: “None”

--figure_dir

If provided, a figure will be created for each std file created

Default: “no_figures”

--cartopy_dir

Directory for cartopy shape files

Default: “None”

--figure_format

File format of figure

Default: “gif”

--log_level

minimum level of messages printed to stdout and in log files

Default: “INFO”

domutils.radar_tools.exponential_zr.exponential_zr(data, coef_a=None, coef_b=None, dbz_to_r=False, r_to_dbz=False, missing=-9999.0, undetect=-3333.0)

Conversion between dBZ and precip rate (R) using exponential Z-R

The relation being used is Z = aR^b

When no coefficients are provided, use WSR-88D Z-R coefficients.

  • a=300, b=1.4

Use

  • a=200, b=1.6

for the original Marshall-Palmer.

To get dBZ from R we have:

\begin{align*} dBZ &= 10 \log_{10}(Z) \\ &= 10 \log_{10}(a) + 10 b \log_{10}(R) \\ &= U + V \log_{10}(R) \end{align*}

To get R from dBZ:

\begin{equation*} R = \frac{10^{\frac{dBZ}{(10 b)}} }{ 10^{\frac{\log_{10}(a)}{b}}} = \frac{10^{\frac{dBZ}{W}} }{ X} \end{equation*}
Parameters:

  • data (Any) – (array-Like) data to be converted

  • coef_a (Optional[float]) – Coefficient a in Z = aR^b

  • coef_b (Optional[float]) – Coefficient b in Z = aR^b

  • dbz_to_r (Optional[bool]) – When True input data is assumed in dBZ and output is R in mm/h

  • r_to_dbz (Optional[bool]) – When True, input data is assumed in mm/h and is converted to dBZ One of dBZtoR or RtodBZ must be set to True

Returns:

A numpy array of the size of input data containing the converted data.

Example

>>> import domutils.radar_tools as radar_tools
>>> import numpy as np
>>> #reflectivity values
>>> ref_dbz  = np.array([-10., -5., 0., 5, 10, 20, 30, 40, 50, 60])
>>> #convert to precip rate im mm/h
>>> rate_mmh = radar_tools.exponential_zr(ref_dbz, dbz_to_r=True)
>>> with np.printoptions(precision=3, suppress=True):
...     print(rate_mmh)
[  0.003   0.007   0.017   0.039   0.088   0.456   2.363  12.24   63.395
 328.354]
>>>
>>> #convert back to dBZ
>>> recovered_dbz= radar_tools.exponential_zr(rate_mmh, r_to_dbz=True)
>>> print(recovered_dbz)
[-10.  -5.   0.   5.  10.  20.  30.  40.  50.  60.]
domutils.radar_tools.read_h5_composite.read_h5_composite(odim_file=None, latlon=False, qced=True, no_data=-9999.0, undetect=-3333.0, latlon_file=None, verbose=0)

Read reflectivity and quality index from OdimH5 composite

Odim H5 files are a bit annoying in that datasets and quality index are just numbered (data1, data2, …) and not named following what they contain. Also, one cannot expect that a dataset named ‘quality3’ will always contain the same quantity. The result is that one has to loop over each item and check that the contents match what we are looking for.

This routine does that as well as converting the Byte data found in H5 files to more tangible quantities such as reflectivity in dBZ and quality index between 0. and 1. Special values are also assigned to no_data and undetect values. UTC time zone is assumed for datestamp in H5 file

Parameters:

  • odim_file (Optional[str]) – /path/to/odim/composite.h5

  • latlon (Optional[bool]) – When true, will output latitudes and longitudes of the odim grid

  • qced (Optional[bool]) – When true (default), Quality Controlled reflectivity (DBZH) will be returned. When false, raw reflectivity (TH) will be returned.

  • no_data (Optional[float]) – Value that will be assigned to missing values

  • undetect (Optional[float]) – Value that will be assigned to valid measurement of no precipitation

  • latlon_file (Optional[str]) – pickle file containing latlon of domain (only used if latlon=True)

Returns:

If no or invalid file present

or a dictionary containing:

’reflectivity’: (ndarray) 2D reflectivity

’total_quality_index’: (ndarray) 2D quality index

’valid_date’: (python datetime object) date of validity

’latitudes’: (ndarray) 2d latitudes of data (conditional on latlon == True)

’longitudes’: (ndarray) 2d longitudes of data (conditional on latlon == True)

Return type:

None

Example

>>> #read odim H5 file
>>> import os, inspect
>>> import domutils.radar_tools as radar_tools
>>> currentdir = os.path.dirname(os.path.abspath(inspect.getfile(inspect.currentframe())))
>>> parentdir = os.path.dirname(currentdir) #directory where this package lives
>>> out_dict = radar_tools.read_h5_composite(parentdir + '/test_data/odimh5_radar_composites/2019/10/31/qcomp_201910311630.h5')
>>> h5_reflectivity        = out_dict['reflectivity']
>>> h5_total_quality_index = out_dict['total_quality_index']
>>> h5_valid_date          = out_dict['valid_date']
>>> print(h5_reflectivity.shape)
(2882, 2032)
>>> print(h5_valid_date)
2019-10-31 16:30:00+00:00
domutils.radar_tools.read_h5_vol.read_h5_vol(odim_file=None, latlon=False, elevations='all', quantities='all', include_quality=False, no_data=-9999.0, undetect=-3333.0)

Read reflectivity and quality index from OdimH5 composite

Odim H5 files are a bit annoying in that datasets and quality index are just numbered (data1, data2, …) and not named following what they contain. Also, one cannot expect that a dataset named ‘quality3’ will always contain the same quantity. The result is that one has to loop over each item and check that the contents match what we are looking for. This routine does that as well as converting the Byte data found in H5 files to more tangible quantities such as reflectivity in dBZ and quality index between 0. and 1. Special values are also assigned to no_data and undetect values. UTC time zone is assumed for datestamp in H5 file

This code is for Volume scans in ODIM format. Reading these files is a bit different from reading mosaics so a different reader is necessary

For more verbose outputs, set logging level in calling handler

Parameters:

  • odim_file (Optional[str]) – /path/to/odim/composite.h5

  • latlon (Optional[bool]) – When true, will output latitudes and longitudes of the returned PPIs

  • elevations (Optional[Any]) – float or list of float for the desired nominal elevation. Set to ‘all’ to return all elevations in a file

  • quantities (include_quality Quality field will be included along side) – desired quantities [‘dbzh’, ‘th’, ‘rhohv’, ‘uphidp’, ‘wradh’, ‘phidp’, ‘zdr’, ‘kdp’, ‘sqih’, ‘vradh’, ‘dr’, etc] set to ‘all’ to return everything

  • quantities –

  • no_data (Optional[float]) – Value that will be assigned to missing values

  • undetect (Optional[float]) – Value that will be assigned to valid measurement of no precipitation

Returns:

If no or invalid file present or desired elevation or quantities not found.

or a dictionary which mimics the structure of odim h5 files:

# dict['elevation1']['quantity1'](2D numpy array, PPI)
#                   ['quantity2'](2D numpy array, PPI)
#                      ...
#                   ['quality_name1'](2D numpy array, PPI)
#                   ['quality_name2'](2D numpy array, PPI)
#                      ...
#                   [nominal_elevation] float
#                   [elevations] (1D numpy array, rows of values)
#                   [azimuths]   (1D numpy array, rows of values)
#                   [ranges]     (1D numpy array, columns of values)
#                   [latitudes]  (2D numpy array)
#                   [longitudes] (2D numpy array)
# dict['elevation2']['quantity1'](2D numpy array)
#                     ...

Return type:

None

Example

>>> #read odim H5 file
>>> import os, inspect
>>> import domutils.radar_tools as radar_tools
>>> currentdir = os.path.dirname(os.path.abspath(inspect.getfile(inspect.currentframe())))
>>> parentdir = os.path.dirname(currentdir) #directory where this package lives
>>> odim_file = parentdir+'/test_data//odimh5_radar_volume_scans/2019071120_24_ODIMH5_PVOL6S_VOL.qc.casbv.h5'
>>> res = radar_tools.read_h5_vol(odim_file=odim_file,
...                               elevations=[0.4],
...                               quantities=['all'],
...                               include_quality=True,
...                               latlon=True)
>>> #check returned radar dictionary
>>> print(res.keys())
dict_keys(['radar_height', 'radar_lat', 'radar_lon', 'date_str', '0.4'])
>>> #
>>> #check returned PPI dictionary
>>> print(res['0.4'].keys())
dict_keys(['dbzh', 'vradh', 'th', 'rhohv', 'uphidp', 'wradh', 'phidp', 'zdr', 'dr', 'kdp', 'sqih', 'quality_beamblockage', 'quality_att', 'quality_broad', 'quality_qi_total', 'nominal_elevation', 'azimuths', 'elevations', 'ranges', 'latitudes', 'longitudes', 'm43_heights'])
domutils.radar_tools.read_fst_composite.read_fst_composite(fst_file=None, valid_date=None, latlon=False, no_data=-9999.0, undetect=-3333.0, verbose=0)

Read reflectivity or precip_rate from CMC standard files

Validity date is obtained via the datev attribute of the entry in the standard file being read. Quality index is set to 1 wherever data is not missing

Parameters:

  • fst_file (Optional[str]) – /path/to/fst/composite.std .stnd .fst or no ‘extention’

  • valid_date (Optional[datetime]) – Datetime object for the time where observations are desired

  • latlon (Optional[bool]) – When true, will output latitudes and longitudes

  • no_data (Optional[float]) – Value that will be assigned to missing values

  • undetect (Optional[float]) – Value that will be assigned to valid measurement of no precipitation

Returns:

If no or invalid file present

or a dictionary containing:

’reflectivity’: (ndarray) 2D reflectivity

’total_quality_index’: (ndarray) 2D quality index

’valid_date’: (python datetime object) date of validity

’latitudes’: (ndarray) 2d latitudes of data (conditional on latlon = True)

’longitudes’: (ndarray) 2d longitudes of data (conditional on latlon = True)

Return type:

None

Example

>>> #read fst file
>>> import os, inspect
>>> import domutils.radar_tools as radar_tools
>>> currentdir = os.path.dirname(os.path.abspath(inspect.getfile(inspect.currentframe())))
>>> parentdir = os.path.dirname(currentdir) #directory where this package lives
>>> out_dict = radar_tools.read_fst_composite(parentdir + '/test_data/std_radar_mosaics/2019103116_30ref_4.0km.stnd')
>>> reflectivity        = out_dict['reflectivity']
>>> total_quality_index = out_dict['total_quality_index']
>>> valid_date          = out_dict['valid_date']
>>> print(reflectivity.shape)
(1650, 1500)
>>> print(valid_date)
2019-10-31 16:30:00+00:00
domutils.radar_tools.median_filter.apply_inds(data, inds)

Apply inds computed from speckle_inds

This is just a shortcut to a 1 line piece of code:

data.ravel()[indices] = filtered bla

domutils.radar_tools.median_filter.get_inds(data, window=3)

Indices for speckle filtering through median filter

Applying filter is just a matter of “data.ravel()[indices] = filtered” which is performed by the function apply_inds() below.

Data points at the edge are repeated for filtering near the domain border

Returning indices allows to apply the same median filtering on the quality index as on the precipitation data itself.

Parameters:

  • data – source data, should be 2D ndarray

  • window – Size of boxcar window default is 3x3 Should be an odd integer

Returns:

indices for applying speckle filter data[indices] = filtered

Example

>>> import numpy as np
>>> import domutils.radar_tools as radar_tools
>>> data = np.array([[1,2,3,0],
...                  [8,9,4,0],
...                  [7,6,5,0]])
>>> inds = radar_tools.median_filter.get_inds(data)
>>> print(inds)
[[ 1  2  1  7]
 [ 8 10  2 11]
 [ 8  9 10 11]]
>>> #median filtered data on a 3x3 window
>>> print(data.ravel()[inds])
[[2 3 2 0]
 [7 5 3 0]
 [7 6 5 0]]
domutils.radar_tools.median_filter.remap_data(data, window=3, mode='clip', return_flat_inds=False)

Puts 2D data into 3D array for any type of window (boxcar) operations

In compiled code, median filtering and other types of boxcar filtering is performed with nested loops which are very inefficient in Python.

This function remaps original data of shape M x N to a 3D array of shape M x N x W with W being the size of the filtering window.

At each point i,j, the values along the 3rd dimension are all the data values participating in the window averaging.

Numpy operations can then be leveraged for fast computation of these filters.

Parameters:

  • data – source data, should be 2D ndarray

  • window – Size of boxcar window default is 3x3 Should be an odd integer

  • mode – mode argument to np.ravel_multi_index Use ‘clip’ (default) to repeat values found at the edge use ‘wrap’ to wrap around

Returns:

a M x N x W remapping of the original data

if return_flat_inds == True, we return the following:

remapped_data: a M x N x W remapping of the original data row_inds: a M x N array containg i index of output array col_inds: a M x N array containg j index of output array flat_inds: a M x N x W array containg flat indices toinput data

Return type:

remapped_data

Example

>>> import numpy as np
>>> import domutils.radar_tools as radar_tools
>>> data = np.array([[1,2,3,0],
...                  [8,9,4,0],
...                  [7,6,5,0]])
>>> remapped_data = radar_tools.median_filter.remap_data(data, window=3, mode='clip')
>>> #
>>> #data points in window for point i=0, j=0
>>> print(remapped_data[0,0,:])
[1. 1. 2. 1. 1. 2. 8. 8. 9.]
>>> #
>>> print(remapped_data[1,0,:])
[1. 1. 2. 8. 8. 9. 7. 7. 6.]
>>> #
>>> #data points in window for point i=1, j=1
>>> print(remapped_data[1,1,:])
[1. 2. 3. 8. 9. 4. 7. 6. 5.]
>>> #
>>> #sum of data points in boxcar
>>> box_sum = np.sum(remapped_data, axis=2)
>>> print(box_sum)
[[33. 33. 23. 10.]
 [49. 45. 29. 12.]
 [65. 57. 35. 14.]]
>>> #
>>> #boxcar average
>>> box_avg = np.mean(remapped_data, axis=2)
>>> print(box_avg)
[[3.66666667 3.66666667 2.55555556 1.11111111]
 [5.44444444 5.         3.22222222 1.33333333]
 [7.22222222 6.33333333 3.88888889 1.55555556]]