ArcticSnow / TopoPyScale

TopoPyScale: a Python library to perform simplistic climate downscaling at the hillslope scale
https://topopyscale.readthedocs.io
MIT License
39 stars 9 forks source link

topo_scale.py downscaling fails on last itersation of SW/LW section #48

Closed joelfiddes closed 1 year ago

joelfiddes commented 1 year ago

config:

split:
    IO: False       # Flag to split downscaling in time or not
    time: 2         # number of years to split timeline in
    space: None     # NOT IMPLEMENTED

fails at this line (l.327):

horizon = horizon_da.sel(x=row.x, y=row.y, azimuth=np.rad2deg(ds_solar.azimuth.isel(point_id=31)), method='nearest')

with:

IndexError: Index is not smaller than dimension 31 >= 31

Explanation:

The object ds_solar.azimuth is:

Out[2]: 
<xarray.DataArray 'azimuth' (point_id: 31, time: 50376)>
dask.array<open_dataset-6baaeed24838583c892469b5b4861c12azimuth, shape=(31, 50376), dtype=float32, chunksize=(31, 50376), chunktype=numpy.ndarray>
Coordinates:
  * point_id  (point_id) int64 1 2 3 4 5 6 7 8 9 ... 23 24 25 26 27 28 29 30 31
  * time      (time) datetime64[ns] 2016-09-01 ... 2022-05-31T23:00:00

point_id = 1:31

the iterator here is the issue (l.266):

for i, row in df_centroids.iterrows():

makes values 1:31 SHOULD be 0:30

this works: ds_solar.azimuth.isel(point_id=0)

this causes fail: ds_solar.azimuth.isel(point_id=31)

More info:

        print('Downscaling LW, SW for point: {} out of {}'.format(pt_id+1,
                                                                  df_centroids.point_id.max()+1))

this print statement assumes a python index as adds one. for this example print statemnts give iterrations 2:32, should be 1:31.

I guess this is connected to the last changes of making downscaling split in time?

Solution:

add -1 here:

horizon = horizon_da.sel(x=row.x, y=row.y, azimuth=np.rad2deg(ds_solar.azimuth.isel(point_id=pt_id-1)), method='nearest')

All other incidences of pt_id are correct with current indexing.

removed +1 here:

        print('Downscaling LW, SW for point: {} out of {}'.format(pt_id,
                                                                  df_centroids.point_id.max()))

to get correct index in print statement

ArcticSnow commented 1 year ago

This is odd the cluster index (point_id) 0 is not present in the ds_solar dataset. It should be there as this means one cluster is then missing.