How do I index my data according to the years

[DOptimusPrime]

I have this dataset. I want to calculate the wind speed for each year. The data is from 2000 to 2013. I can calculate the wind speed.
Is there any a more efficient way that I can index my dataset. The code I have right now I have to rewrite the indexing part for each indvidivual year manual. Is there a way better way I can do the indexing and break the data into individual years?

opens all the netcdf4 files containing eu10 and ev10

#present climate data
present_eu10 = xr.open_mfdataset('/chinook2/wgallus/climate_output/EU10.nc')
present_ev10 = xr.open_mfdataset('/chinook2/wgallus/climate_output/EV10.nc')

#future climate data
future_eu10 = xr.open_mfdataset('/chinook2/wgallus/climate_output_future/PGW_EU10.nc')
future_ev10= xr.open_mfdataset('/chinook2/wgallus/climate_output_future/PGW_EV10.nc')

Indexing the data according to the different years.

year2000eu = present_eu10['EU10'].isel(Time = slice(0, 2208))
year2000ev = present_ev10['EV10'].isel(Time = slice(0, 2208))
future2000eu = future_eu10['EU10'].isel(Time = slice(0, 2208))
future2000ev = future_ev10['EV10'].isel(Time = slice(0, 2208))

ds1_2001 = present_eu10['EU10'].isel(Time = slice(2208, 10968))
ds2_2001 = present_ev10['EV10'].isel(Time = slice(2208, 10968))
ds3_2001 = future_eu10['EU10'].isel(Time = slice(2208, 10968))
ds4_2001 = future_ev10['EV10'].isel(Time = slice(2208, 10968))

ds1_2002 = present_eu10['EU10'].isel(Time = slice(10968, 19728))
ds2_2002 = present_ev10['EV10'].isel(Time = slice(10968, 19728))
ds3_2002 = future_eu10['EU10'].isel(Time = slice( 10968, 19728))
ds4_2002 = future_ev10['EV10'].isel(Time = slice( 10968, 19728))

ds1_2003 = present_eu10['EU10'].isel(Time = slice(19728, 28488))
ds2_2003 = present_ev10['EV10'].isel(Time = slice(19728, 28488))
ds3_2003 = future_eu10['EU10'].isel(Time = slice(19728,28488))
ds4_2003 = future_ev10['EV10'].isel(Time = slice(19728,28488))

ds1_2004 = present_eu10['EU10'].isel(Time = slice(28488, 37272))
ds2_2004 = present_ev10['EV10'].isel(Time = slice(28488, 37272))
ds3_2004 = future_eu10['EU10'].isel(Time = slice(28488, 37272))
ds4_2004 = future_ev10['EV10'].isel(Time = slice(28488, 37272))

ds1_2005 = present_eu10['EU10'].sel(Time = "2005")
ds2_2005 = present_ev10['EV10'].sel(Time = "2005")
ds3_2005 = future_eu10['EU10'].sel(Time = "2005")
ds4_2005 = future_ev10['EV10'].sel(Time = "2005")

ds1_2006 = present_eu10['EU10'].sel(Time = "2006")
ds2_2006 = present_ev10['EV10'].sel(Time = "2006")
ds3_2006 = future_eu10['EU10'].sel(Time = "2006")
ds4_2006 = future_ev10['EV10'].sel(Time = "2006")

ds1_2007 = present_eu10['EU10'].sel(Time = "2007")
ds2_2007 = present_ev10['EV10'].sel(Time = "2007")
ds3_2007 = future_eu10['EU10'].sel(Time = "2007")
ds4_2007 = future_ev10['EV10'].sel(Time = "2007")

ds1_2008 = present_eu10['EU10'].sel(Time = "2008")
ds2_2008 = present_ev10['EV10'].sel(Time = "2008")
ds3_2008 = future_eu10['EU10'].sel(Time = "2008")
ds4_2008 = future_ev10['EV10'].sel(Time = "2008")

ds1_2009 = present_eu10['EU10'].sel(Time = "2009")
ds2_2009 = present_ev10['EV10'].sel(Time = "2009")
ds3_2009 = future_eu10['EU10'].sel(Time = "2009")
ds4_2009 = future_ev10['EV10'].sel(Time = "2009")

ds1_2010 = present_eu10['EU10'].sel(Time = "2010")
ds2_2010 = present_ev10['EV10'].sel(Time = "2010")
ds3_2010 = future_eu10['EU10'].sel(Time = "2010")
ds4_2010 = future_ev10['EV10'].sel(Time = "2010")

ds1_2011 = present_eu10['EU10'].sel(Time = "2011")
ds2_2011 = present_ev10['EV10'].sel(Time = "2011")
ds3_2011 = future_eu10['EU10'].sel(Time = "2011")
ds4_2011 = future_ev10['EV10'].sel(Time = "2011")

ds1_2012 = present_eu10['EU10'].sel(Time = "2012")
ds2_2012 = present_ev10['EV10'].sel(Time = "2012")
ds3_2012 = future_eu10['EU10'].sel(Time = "2012")
ds4_2012 = future_ev10['EV10'].sel(Time = "2012")

ds1_2013 = present_eu10['EU10'].sel(Time = "2013")
ds2_2013 = present_ev10['EV10'].sel(Time = "2013")
ds3_2013 = future_eu10['EU10'].sel(Time = "2013")
ds4_2013 = future_ev10['EV10'].sel(Time = "2013")

[keewis]

this sounds like a groupby operation:

present_eu10.groupby("Time.year")

Does that help? If not, can you give a bit more context on the datasets and what you're doing with the extracted groups?

[DOptimusPrime]

Thank you for the answer. What I am trying to do is this. I have this data and I want to calculate the wind speed and get the maximum wind speed for each year. Here is a picture of what is in the folder. I opened all the files containing eu10 and ev10. To calculate the speed is : (eu10 2 + ev102) ** 0.5 . "groupby" allows me to group them by year. But then I cannot calculate the speed. I get a TypeError
I hope the folliowing will give you more context

opens all the netcdf4 files containing eu10 and ev10

#present climate data
present_eu10 = xr.open_mfdataset('/chinook2/wgallus/climate_output/EU10.nc').groupby("Time.year")
present_ev10 = xr.open_mfdataset('/chinook2/wgallus/climate_output/EV10.nc').groupby("Time.year")
#future climate data
future_eu10 = xr.open_mfdataset('/chinook2/wgallus/climate_output_future/PGW_EU10.nc').groupby("Time.year")
future_ev10= xr.open_mfdataset('/chinook2/wgallus/climate_output_future/PGW_EV10.nc').groupby("Time.year")

speed = (present_eu10 **2 + present_ev10 ** 2) ** 0.5
TypeError Traceback (most recent call last)
Cell In[29], line 1
----> 1 speed = (present_eu10 **2 + present_ev10 ** 2) ** 0.5

TypeError: GroupBy objects only support binary ops when the other argument is a Dataset or DataArray

This is what in my folder. I open only the files containing eu10 and ev10
for data in (glob.glob('/chinook2/wgallus/climate_output_future/*.nc')):
print(data)

/chinook2/wgallus/climate_output_future/wrf2d_d01_PGW_EU10_200010-200012.nc
/chinook2/wgallus/climate_output_future/wrf2d_d01_PGW_EV10_200010-200012.nc
/chinook2/wgallus/climate_output_future/wrf2d_d01_PGW_PREC_ACC_NC_200010-200012.nc
/chinook2/wgallus/climate_output_future/wrf3d_d01_PGW_QRAIN_200010.nc
/chinook2/wgallus/climate_output_future/wrf3d_d01_PGW_QRAIN_200011.nc
/chinook2/wgallus/climate_output_future/wrf3d_d01_PGW_QRAIN_200012.nc
/chinook2/wgallus/climate_output_future/wrf2d_d01_PGW_T2_200010-200012.nc
/chinook2/wgallus/climate_output_future/wrf2d_d01_PGW_EU10_200101-200103.nc
/chinook2/wgallus/climate_output_future/wrf2d_d01_PGW_EV10_200101-200103.nc
/chinook2/wgallus/climate_output_future/wrf2d_d01_PGW_PREC_ACC_NC_200101-200103.nc
/chinook2/wgallus/climate_output_future/wrf3d_d01_PGW_QRAIN_200101.nc
/chinook2/wgallus/climate_output_future/wrf3d_d01_PGW_QRAIN_200102.nc
/chinook2/wgallus/climate_output_future/wrf3d_d01_PGW_QRAIN_200103.nc
/chinook2/wgallus/climate_output_future/wrf2d_d01_PGW_T2_200101-200103.nc
/chinook2/wgallus/climate_output_future/wrf2d_d01_PGW_EU10_200104-200106.nc
/chinook2/wgallus/climate_output_future/wrf2d_d01_PGW_EV10_200104-200106.nc
/chinook2/wgallus/climate_output_future/wrf2d_d01_PGW_PREC_ACC_NC_200104-200106.nc
/chinook2/wgallus/climate_output_future/wrf3d_d01_PGW_QRAIN_200104.nc
/chinook2/wgallus/climate_output_future/wrf3d_d01_PGW_QRAIN_200105.nc
/chinook2/wgallus/climate_output_future/wrf3d_d01_PGW_QRAIN_200106.nc
/chinook2/wgallus/climate_output_future/wrf2d_d01_PGW_T2_200104-200106.nc
/chinook2/wgallus/climate_output_future/wrf2d_d01_PGW_EU10_200107-200109.nc
/chinook2/wgallus/climate_output_future/wrf2d_d01_PGW_EV10_200107-200109.nc

[DanielAdriaansen]

Can you provide the output from:
print(present_eu10)

and

print(present_eu10['Time'])

[DanielAdriaansen]

Hi,

I thought that datetime had a timedelta object of "years" but I was wrong:
https://docs.python.org/3/library/datetime.html#timedelta-objects

You'll somehow need to increment time forward the duration you want, since if you use days=365, I don't think you'll account for leap years.

Perhaps a better approach is to create a list of years you'd like ahead of time like:
years = [str(x) for x in range(2000,2014,1)]

And then do a for loop over each year in that list instead of the while loop I suggested like:

years = [str(x) for x in range(2000,2014,1)]
for current_year in years:

  # Calculate the wind speed
  speed = (present_eu10.sel(Time=current_year) **2 + present_ev10.sel(Time=current_year) ** 2) ** 0.5

Or, the suggestion below from @keewis is also very good. Ultimately there are several approaches to accomplish what you want to do I believe.

[DOptimusPrime]

Thank you for the suggestion both of them gave me errors.
With this one in particular is it saying that I am trying to perform an operation on something that is not an Array or integer
speed = (present_eu10.sel(Time=current_year) **2 + present_ev10.sel(Time=current_year) ** 2) ** 0.5

[DOptimusPrime]

I don't know why but I think sel() and isel() both return Array but in this case that is the error I get

[keewis]

yeah, this is pretty weird. Can you tell us the output of

print(present_eu10["EU10"].data)

?

All arrays that xarray can wrap should support __pow__ with integers.

It would also be useful to know what your environment is. Could you tell us the output of xr.show_versions()?

[DOptimusPrime]

print(present_eu10["EU10"].data)
dask.array<concatenate, shape=(113952, 1015, 1359), dtype=float32, chunksize=(2208, 1015, 1359), chunktype=numpy.ndarray>

This is what I use
https://jupyterhub.agron.iastate.edu/

[keewis]

If I understand correctly, what you want to do is

1. compute the magnitude of the wind speed from wind vectors
2. get the maximum wind speed per year
3. compute how many times per year the magnitude exceeds a certain threshold

Breaking the task up this way makes the code very easy to read and understand:

wind_speed = np.sqrt(present_eu10["EU10"] ** 2 + present_ev["EV10"] ** 2)
maximum_per_year = wind_speed.groupby("Time.year").max()
n_exceeded = (wind_speed > threshold).groupby("Time.year").sum()

if threshold is something like DataArray([10, 15, 20], dims="thresholds"), this will allow you to compute all thresholds in a single operation. If you're instead interested in binning the wind speeds: that is possible as well, but unfortunately somewhat tricky at the moment (we're working on that, though).

[DOptimusPrime]

How do I retrieve the data that is:
get the maximum wind speed per year
compute how many times per year the magnitude exceeds a certain threshold

[keewis]

you mean, because the data is still a dask array, so the computation didn't actually happen, yet? Just call .compute() on it. What I would do is collect everything in a Dataset and compute everything together:

wind_speed = np.sqrt(present_eu10["EU10"] ** 2 + present_ev["EV10"] ** 2)
maximum_per_year = wind_speed.groupby("Time.year").max()
n_exceeded = (wind_speed > threshold).groupby("Time.year").sum()

result = xr.Dataset({"maximum": maximum_per_year, "n_exceeded": n_exceeded}, coords={"threshold": threshold})

# to look at the values:
computed = result.compute()
computed

[DOptimusPrime]

I have been running it for about an hour an half now and still no output. I will paste what I get when it is ready.

[DOptimusPrime]

[keewis]

well, that can happen if your data is big (you can see that by displaying result and looking at the array's data repr). If you want to know more about dask, see the dask documentation and this section from the xarray docs.

[DOptimusPrime]