<jupyter_start><jupyter_text>Light calculationsHere we use pysolar python modules to perform the light calculations<jupyter_code>latitude_deg=37.77
longitude_deg=-122.41
latitudes=np.arange(0,90,1)

d = datetime.datetime.now(datetime.timezone(offset=datetime.timedelta(hours=-8)))
thirty_minutes = datetime.timedelta(hours = 2)
for _ in range(48):
    timestamp = d.ctime()
    altitude_deg = solar.get_altitude(latitude_deg, longitude_deg, d)
    azimuth_deg = solar.get_azimuth(latitude_deg, longitude_deg, d)
    power = radiation.get_radiation_direct(d, altitude_deg)
    
    test3=pysolar.util.global_irradiance_clear(latitude_deg, longitude_deg, d)
    
    
    #if (altitude_deg > 0):
   # print(timestamp, "PST: {:.2f} {:.2f} ".format(power,test3))
    d = d + thirty_minutes     

grid_labels = ['gr'] 

member_ids = ['r1i1p1f1', 'r1i1p1f2', 'r1i1p1f3', 'r2i1p1f1', 'r2i1p1f2', 'r2i1p1f3', 'r3i1p1f1', 'r3i1p1f2', 
    'r4i1p1f1', 'r4i1p1f2', 'r4i1p1f3', 'r5i1p1f1', 'r5i1p1f2', 'r5i1p1f3', 'r6i1p1f1', 'r6i1p1f2', 'r6i1p1f3', 'r7i1p1f1', 'r7i1p1f2',
    'r7i1p1f3', 'r8i1p1f1', 'r8i1p1f2', 'r8i1p1f3', 'r9i1p1f1', 'r9i1p1f2', 'r9i1p1f3', 'r102i1p1f1', 'r102i1p1f2', 'r102i1p1f3']
member_ids = ['r1i1p1f1']

experiment_ids = ['1pctCO2'] #'abrupt-4xCO2',
institution_ids = [ 'NOAA-GFDL'] 
source_ids = ['GFDL-ESM4']
variable_ids = ["sithick","siconc","sisnthick","sisnconc"]
table_ids=['SImon','SImon','SImon','SImon']
#variable_ids = ["sithick"]<jupyter_output><empty_output><jupyter_text>Query dataThe dataset returned from the query need to be checked for time dimension. If the time array uses a relative date (0001 instead of 20001) then the values need to be added 2000 prior to converting to datetime64 format.https://unidata.github.io/cftime/api.htmlhttp://xarray.pydata.org/en/stable/time-series.html<jupyter_code>def plot_timeseries(ds_dict,ds_area):

    plt.figure(figsize=(12,8))

    for name, ds in ds_dict.items():
    
      #  total_area = ds_area.areacella.sum(dim=['lat', 'lon'])
     #   print(ds_area.areacella)
        if name=='sithick':
            ta_timeseries = ds.sithick.mean(dim=['lat','lon']) #* ds_area.areacella).sum(dim=['lon', 'lat']) / total_area
        if name=='sisnthick':
            ta_timeseries = ds.sisnthick.mean(dim=['lat','lon']) #* ds_area.areacella).sum(dim=['lon', 'lat']) / total_area
        if name=='siconc':
            ta_timeseries = ds.siconc.mean(dim=['lat','lon']) #* ds_area.areacella).sum(dim=['lon', 'lat']) / total_area
        if name=='sisnconc':
            ta_timeseries = ds.sisnconc.mean(dim=['lat','lon']) #* ds_area.areacella).sum(dim=['lon', 'lat']) / total_area
        if name=='tas':
            ta_timeseries = ds.tas.mean(dim=['lat','lon']) #* ds_area.areacella).sum(dim=['lon', 'lat']) / total_area
         
        plt.ylabel(r'global-mean {}'.format(name))
        ta_timeseries.plot()
        ta_timeseries.rolling(time=12).mean().plot(color='r')

        plt.xlabel('time')
     #   plt.xlim([1850,2100]);
        plt.show()
  #  plot_map(ds)
    #plt.savefig('../figures/ssp585_global_warming.png',dpi=100,bbox_inches='tight')
def plot_map(ds):
    plt.figure(figsize=(14,6))
    ax = plt.axes(projection=ccrs.PlateCarree())
    ax.set_global()
    ds.sithick[0].plot.pcolormesh(ax=ax, transform=ccrs.PlateCarree(), x='xc', y='yc', add_colorbar=False)
    ax.coastlines()
    ax.set_ylim([0,90]);
    plt.show()
def perform_query(query_string):
    df_sub = df.query(query_string)
    mapper = fs.get_mapper(df_sub.zstore.values[-1])
    ds = xr.open_zarr(mapper, consolidated=True)
    time_object=ds['time'].values[0]
    
    # Convert if necesssary
    if (time_object.year==1):
      
        times=ds['time'].values
        times_plus_2000=[]
        for t in times:
            times_plus_2000.append(cftime.DatetimeNoLeap(t.year+2000,t.month,t.day,t.hour))
        ds['time'].values=times_plus_2000
        ds = xr.decode_cf(ds)
    print("=> Dates extracted range from {} to {}\n".format(ds['time'].values[0],ds['time'].values[-1])) 
    return ds
from dask.distributed import Client, progress
client = Client(processes=False) #threads_per_worker=2, n_workers=16)
#client.cluster.scale(2)
client()
def calculate_light(lon,lat,when,snow_thickness,ice_thickness):
    altitude_deg = solar.get_altitude(lat, lon, when)
    surface_light = radiation.get_radiation_direct(when, altitude_deg)
    
    return calculate_light_under_sea_ice(snow_thickness,ice_thickness,surface_light)
    
def calculateArea(lat0,lat1,lon0,lon1,areaIce):

    earthRadius = 6371000
    rad = np.pi / 180.0

    """    -180 <= lon0 < lon1 <= 180
            -90 <= lat0 < lat1 <= 90
            areaIce is in percent
    """

    area = earthRadius**2 * (np.sin(lat1*rad)-np.sin(lat0*rad)) * (lon1 - lon0) * rad
    return area * (areaIce)
def calculate_light_under_sea_ice(snow_thickness,ice_thickness,surface_light,debug=False):
    # Now calculate the amount of light below the sea ice and snow (Perovich 1996,
    # but see  Jin 2006 Annals of Glaciology)
    attenuationSnow = 20  #unit : m-1
    attenuationIceTop10cm = 5
    attenuationIceBelowSurface = 1
    missing_value=1e20
    
    if 0 < snow_thickness  < missing_value:
        albedo = 0.9
        if debug is True:
            print("Albedo for snow and ice covered water: {}".format(albedo))
        surface_light = (1.0 - albedo) * surface_light
    if snow_thickness == 0 and 0 < ice_thickness < missing_value:
        albedo = 0.5
        if debug is True:
            print("Albedo for ice covered water: {}".format(albedo))
        surface_light = (1.0 - albedo) * surface_light
    if (snow_thickness == 0 and ice_thickness == 0) or (snow_thickness==missing_value and ice_thickness==missing_value):
        albedo = 0.06
        if debug is True:
            print("Albedo for open water: {}".format(albedo))
        surface_light = (1.0 - albedo) * surface_light

    Eb = surface_light
    if debug is True:
        print("\nSurface light {}".format(Eb))
    if 0 < snow_thickness < missing_value:
        Eb = surface_light * np.exp(attenuationSnow * (-snow_thickness))
        if debug is True:
            print("Eb with snow (%s m) : {}".format(snow_thickness, Eb))

    if  0.1 <= ice_thickness < missing_value:

        Eb = Eb * np.exp(attenuationIceTop10cm * (-0.1))
        if debug is True:
            print("Eb with ice top (%s m) : {}".format(ice_thickness, Eb))
        Eb = Eb * np.exp(attenuationIceBelowSurface * (-(ice_thickness - 0.1)))
        if debug is True:
            print("Eb with ice below top (%s m) : {}".format(ice_thickness - 0.1, Eb))
    else:
        Eb = Eb * np.exp(attenuationIceTop10cm * (-ice_thickness))
        if debug is True:
            print("Eb with ice top (%s m) : {}".format(ice_thickness, Eb))

    #print "Eb", Eb, "snow", snowthickness, "ice", icethickness, "albedo", albedo

   # if snowthickness==missing_value or icethickness==missing_value or albedo==missing_value:
   #     return missing_value
    return Eb
def run_tasks_for_dataset(dset_dict):
    #client.cluster.scale(10)  # ask for ten 4-thread workers
    
    snow_thickness=dset_dict['sisnthick'].sisnthick
    ice_thickness=dset_dict['sithick'].sithick
   
    dd=snow_thickness['time'].values[0]
    when=datetime.datetime(dd.year,dd.month,dd.day,tzinfo=datetime.timezone.utc)

    z=np.zeros((len(snow_thickness.lon.values),len(snow_thickness.lat.values)))
    n, m = np.shape(z)
    
    zr = [delayed(calculate_light)(snow_thickness.lon[i],snow_thickness.lat[j],when,snow_thickness[0,j,i],ice_thickness[0,j,i]) for i in range(n) for j in range(m)]
    final=dask.compute(zr)
    
    
    zr=np.array(final).reshape((n,m))
    print(np.shape(zr))
 #   ds = xr.DataArray(zr, dims=['lon', 'lat'],
  #                         coords={'lon': x, #snow_thickness.lon,
  #                                 'lat': y}) #snow_thickness.lat})
#    ds.to_netcdf('test.nc')

    
    plt.figure(figsize=(12,12))
    plt.imshow(zr)
    plt.colorbar()<jupyter_output><empty_output><jupyter_text>Useful linkshttps://climate-cms.org/2019/11/12/Calendars-and-monthly-data.htmlhttp://xarray.pydata.org/en/stable/examples/monthly-means.htmlhttps://rabernat.github.io/research_computing_2018/xarray-tips-and-tricks.htmlhttp://meteo.unican.es/work/xarray_seminar/xArray_seminar.html<jupyter_code># Dicitionary to hold the queried variables
dset_dict = {}
first=True

for experiment_id in experiment_ids:
    for grid_label in grid_labels:
        for source_id in source_ids:
            for member_id in member_ids:
                for variable_id,table_id in zip(variable_ids,table_ids):
                    
                    # Historical
                    query_string="source_id=='{}'and table_id=='{}' and member_id=='{}' and grid_label=='{}' and experiment_id=='esm-hist' and variable_id=='{}'".format(source_id,table_id,member_id,grid_label,variable_id)
                    print("Running historical query on data: \n ==> {}\n".format(query_string))
                    ds_hist=perform_query(query_string)
                    
                    # Future projection depending on choice in experiment_id
                    query_string="source_id=='{}'and table_id=='{}' and member_id=='{}' and grid_label=='{}' and experiment_id=='{}' and variable_id=='{}'".format(source_id,table_id,member_id,grid_label,experiment_id,variable_id)
                    print("Running projections query on data: \n ==> {}\n".format(query_string))
                    ds_proj=perform_query(query_string)
                   
                    if first:
                        df_area = df.query("variable_id == 'areacella' and source_id =='{}'".format(source_id))
                        ds_area = xr.open_zarr(fs.get_mapper(df_area.zstore.values[0]), consolidated=True)
                        first=False
                    
                    # Concatentate the historical and projections datasets
                    ds = xr.concat([ds_hist, ds_proj],dim='time')
                    
                    # Remove the duplicate overlapping times (e.g. 2001-2014)
                    _, index = np.unique(ds['time'], return_index=True)
                    ds=ds.isel(time=index)
                    
                    # Save the dataset for variable_id in the dictionary
                    dset_dict[variable_id]=ds
                   
    
run_tasks_for_dataset(dset_dict)
                
# Plot the global average timeseries for each variable
plot_timeseries(dset_dict,ds_area)<jupyter_output>Running historical query on data: 
 ==> source_id=='GFDL-ESM4'and table_id=='SImon' and member_id=='r1i1p1f1' and grid_label=='gr' and experiment_id=='esm-hist' and variable_id=='sithick'
{'fraction_non_alphanumeric': 0.14082278481012658, 'fraction_numerical': 0.03145941921072226, 'mean_word_length': 3.1168582375478926, 'pattern_counts': {'":': 0, '<': 26, '<?xml version=': 0, '>': 17, 'https://': 3, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '3', 'dir': 'jupyter-structured-clean-dedup', 'id': '566639', 'license': 'MIT', 'n_tokens_mistral': 4116, 'n_tokens_neox': 3813, 'n_words': 673, 'path': '.ipynb_checkpoints/seaice-light-cmip6_v4-checkpoint.ipynb', 'repo_name': 'trondkr/okokyst_toolbox', 'size': '10744'}