Downlowd for TCWV notebook

[CarboneLMD]

Notebook description

Hi @malmans2,

I need to download all of this datas please :

# Set parameters
year_start = 2003
year_stop = 2017
area = [30, -180, -30, 180]



# Define requests
requests = {}

time_request = {
    "year": [str(year) for year in range(year_start, year_stop + 1)],
    "month": [f"{month:02d}" for month in range(1, 13)],
    "day": [f"{day:02d}" for day in range(1, 31 + 1)],
}

requests["satellite-total-column-water-vapour-land-ocean"] = time_request | {
    "product": "near_infrared_hoaps_combined",
    "horizontal_aggregation": ["0_5_x_0_5"],
    "temporal_aggregation": "daily",
    "variable": "all",
}
requests["satellite-cloud-properties"] = time_request | {
    "product_family": "clara_a3",
    "origin": "eumetsat",
    "variable": ["cloud_fraction"],
    "climate_data_record_type": "thematic_climate_data_record",
    "time_aggregation": "daily_mean",
    "area": area,
}
requests["reanalysis-era5-single-levels"] = time_request | {
    "variable": ["skin_temperature"],
    "product_type": ["reanalysis"],
    "time": ["00:00", "01:00", "02:00",
        "03:00", "04:00", "05:00",
        "06:00", "07:00", "08:00",
        "09:00", "10:00", "11:00",
        "12:00", "13:00", "14:00",
        "15:00", "16:00", "17:00",
        "18:00", "19:00", "20:00",
        "21:00", "22:00", "23:00"
    ],
    "area": area,
}
requests["satellite-earth-radiation-budget"] = time_request | {
    "product_family": "clara_a3",
    "origin": "eumetsat",
    "variable": ["outgoing_longwave_radiation"],
    "climate_data_record_type": "thematic_climate_data_record",
    "time_aggregation": "daily_mean",
    "area": area,
}

Is it possible to aggregate "reanalysis-era5-single-levels" into daily_mean ?

Thank you !

Notebook link or upload

No response

Anything else we need to know?

No response

Environment

Details

[malmans2]

Can I use this dataset for daily means of ERA5? https://cds.climate.copernicus.eu/datasets/derived-era5-single-levels-daily-statistics?tab=download

[CarboneLMD]

Yes, i guess. I will see if i find the same results again

[malmans2]

I'm running the script to download data. I'll let you know when it's done.
The chunking you are using in the code you shared is too big. Please don't run the code as it is as I'm caching different chunks.

[malmans2]

Everything is cached. Here is the code to retrieve the data:

from c3s_eqc_automatic_quality_control import download, utils

# Set parameters
year_start = 2003
year_stop = 2017

# Region
lon_slice = slice(-180, 180)
lat_slice = slice(-30, 30)

# Define requests
requests = {}

time_request = {
    "year": [str(year) for year in range(year_start, year_stop + 1)],
    "month": [f"{month:02d}" for month in range(1, 13)],
    "day": [f"{day:02d}" for day in range(1, 32)],
}

requests["satellite-total-column-water-vapour-land-ocean"] = time_request | {
    "product": "near_infrared_hoaps_combined",
    "horizontal_aggregation": ["0_5_x_0_5"],
    "temporal_aggregation": "daily",
    "variable": "all",
}
requests["satellite-cloud-properties"] = time_request | {
    "product_family": "clara_a3",
    "origin": "eumetsat",
    "variable": ["cloud_fraction"],
    "climate_data_record_type": "thematic_climate_data_record",
    "time_aggregation": "daily_mean",
}
requests["derived-era5-single-levels-daily-statistics"] = time_request | {
    "product_type": "reanalysis",
    "variable": ["skin_temperature"],
    "daily_statistic": "daily_mean",
    "time_zone": "utc+00:00",
    "frequency": "1_hourly",
}
requests["satellite-earth-radiation-budget"] = time_request | {
    "product_family": "clara_a3",
    "origin": "eumetsat",
    "variable": ["outgoing_longwave_radiation"],
    "climate_data_record_type": "thematic_climate_data_record",
    "time_aggregation": "daily_mean",
}

# Download
datasets = {}
for collection_id, request in requests.items():
    datasets[collection_id] = download.download_and_transform(
        collection_id,
        request,
        chunks={"year": 1, "month": 1},
        transform_func = utils.regionalise,
        transform_func_kwargs = {
            "lon_slice": lon_slice,
            "lat_slice": lat_slice,
        }
    )

[CarboneLMD]

Hi @malmans2,

Thank's ! it works well.

First of all, i know you have a lot of work to do actually, so you can respond to me later. In fact, i already have all the figures for this notebook, so i can work on it.

I have an other question about this notebook. In fact, all the figures are made with a dataframe. But this dataframe takes a lot of time to be done (~ 20 minutes). So i'm wondering if you know how to make it faster. Here is the beginning of the notebook from the dataframe :

import cartopy.crs as ccrs
import matplotlib.pyplot as plt
import pandas as pd
import xarray as xr
from scipy.stats import linregress
from xhistogram.xarray import histogram

import math
import numpy as np 
from c3s_eqc_automatic_quality_control import diagnostics, download, plot, utils

import os
os.environ["CDSAPI_RC"] = os.path.expanduser("~/carbone_sullivan/.cdsapirc")  

# Set parameters
year_start = 2003
year_stop = 2017

# Region
lon_slice = slice(-180, 180)
lat_slice = slice(-30, 30)

# Define requests
requests = {}

time_request = {
    "year": [str(year) for year in range(year_start, year_stop + 1)],
    "month": [f"{month:02d}" for month in range(1, 13)],
    "day": [f"{day:02d}" for day in range(1, 32)],
}

requests["satellite-total-column-water-vapour-land-ocean"] = time_request | {
    "product": "near_infrared_hoaps_combined",
    "horizontal_aggregation": ["0_5_x_0_5"],
    "temporal_aggregation": "daily",
    "variable": "all",
}
requests["satellite-cloud-properties"] = time_request | {
    "product_family": "clara_a3",
    "origin": "eumetsat",
    "variable": ["cloud_fraction"],
    "climate_data_record_type": "thematic_climate_data_record",
    "time_aggregation": "daily_mean",
}
requests["derived-era5-single-levels-daily-statistics"] = time_request | {
    "product_type": "reanalysis",
    "variable": ["skin_temperature"],
    "daily_statistic": "daily_mean",
    "time_zone": "utc+00:00",
    "frequency": "1_hourly",
}
requests["satellite-earth-radiation-budget"] = time_request | {
    "product_family": "clara_a3",
    "origin": "eumetsat",
    "variable": ["outgoing_longwave_radiation"],
    "climate_data_record_type": "thematic_climate_data_record",
    "time_aggregation": "daily_mean",
}

# Download
datasets = {}
for collection_id, request in requests.items():
    datasets[collection_id] = download.download_and_transform(
        collection_id,
        request,
        chunks={"year": 1, "month": 1},
        transform_func = utils.regionalise,
        transform_func_kwargs = {
            "lon_slice": lon_slice,
            "lat_slice": lat_slice,
        }
    )
 

For the next part below, i transform all the datasets in order to have the same geospatial resolution ( longitudes and latitudes) to merge them in one dataset (ds_daily). I know that i can make it with download fonction.

#coarsen datasets to 0.5° grid 
ds_clouds = datasets["satellite-cloud-properties"][["cfc"]].coarsen(latitude=2, longitude=2).mean()

ds_flux = datasets["satellite-earth-radiation-budget"][["LW_flux"]].coarsen(latitude=2, longitude=2).mean()

ds_tcwv   = datasets["satellite-total-column-water-vapour-land-ocean"][["tcwv", "surface_type_flag"]]

ds_sst = datasets["derived-era5-single-levels-daily-statistics"]
new_lat = datasets["satellite-earth-radiation-budget"].latitude.values
new_lon = datasets["satellite-earth-radiation-budget"].longitude.values
ds_sst = ds_sst.interp(latitude=new_lat, longitude = new_lon)
ds_sst = ds_sst.coarsen(latitude=2, longitude=2).mean()

# Merge datasets
ds_daily = xr.merge(
    [
        ds_clouds,
        ds_flux,
        ds_sst,
        ds_tcwv
    ]
)

Then, i compute clear sky greenhouse effect and i add it in ds_daily. afterwards, i create a clear sky mask in order to filter the dataset ds_daily (keep only clear sky zone) and i compute the monthly mean (ds_monthly) :

# Compute clear sky greenhouse effect and add to the merged dataset
ds_daily["GCS"] = (5.67 * 10 ** (-8) * (ds_daily["skt"] ** 4) - ds_daily["LW_flux"])
ds_daily["GCS"].attrs = {"long_name": "clear sky greenhouse effect", "units": "W m−2"}

# Selection of clear sky regions (at most 10\% of cloud fraction within each grid cell) and quality of retrieved SST and TCVW (SST > 0°K and TCWV > 0 mm)
mask = (ds_daily["cfc"] < 10) & (ds_daily["skt"] > 0) & (ds_daily["tcwv"] > 0) ## & (ds_daily["GCS"] > 120)

ds_daily = ds_daily.where(mask)

# Compute montly mean from daily & masked data
ds_monthly = ds_daily.resample(time="M").mean()

Finally, i transform the monthly dataset into three dataframes based on the surface type (land, ocean and both) :

# variables and surfaces
variables = ["tcwv", "skt", "LW_flux", "GCS"]
surface = ["land", "ocean", "land_ocean"]

df_monthly = {}

for idx, surf in enumerate(surface):
    series = {}

    # selection according to surface_type_flag 
    if idx < 2:
        ds_sel = ds_monthly.where(ds_tcwv["surface_type_flag"] == idx)
    else:
        ds_sel = ds_monthly

    for var in variables:
        da = ds_sel[var]

        df_var = da.to_dataframe(name=var).reset_index()  

        cols = list(df_var.columns)
        index_cols = [c for c in cols if c not in ("realization","latitude", "longitude", var)]
        
        series[var] = df_var.set_index(index_cols)[var]

    # Building the final DataFrame  
    df_monthly[f"df_monthly_{surf}"] = pd.DataFrame(series).dropna().sort_index()

This part makes lot of time.

Thank you for your great help

[malmans2]

Let's proceed step by step. It looks like we need to cache all data on the same grid.
However, I'm not fully following what you are doing with "derived-era5-single-levels-daily-statistics"

It looks like the "derived-era5-single-levels-daily-statistics" grid overlaps the target grid. Can we just select the same grid points of the target grid?

[malmans2]

@CarboneLMD I started a template for this use case.

I got to the point where I open all datasets, regrid them when needed, and merge them into a single dataset.
Then I compute monthly means.

I'm not sure I understand what you're trying to do when you build the dataframes. First, you probably don't need to use pandas at all. Everything can be done more efficiently with xarray. More importantly, I don't understand this line:

ds_sel = ds_monthly.where(ds_tcwv["surface_type_flag"] == idx)

Isn't ds_tcwv a daily variable? Why are you applying a daily mask to a monthly dataset?

Here is the template: https://github.com/bopen/c3s-eqc-toolbox-template/blob/main/notebooks/wp5/tcwv.ipynb

[CarboneLMD]

Hi @malman2 ,

Thank you !

Yes, we can select the same grid points of the target grid.

ds_tcwv["surface_type_flag"] is just a surface type (land, ocean ...), so daily or monthly, it's the same thing i guess. May be, i should rename it ?

[malmans2]

As far as I can tell it's not the same. surface_type_flag has the time dimension and changes over time, therefore the monthly average is different compared to the daily values. If you mask a monthly dataset with a daily dataset, then you are broadcasting the monthly to daily.

How would you like to perform the masking?

[CarboneLMD]

HI @malmans2,

Happy new year !

Sorry for my late reply, i was busy finishing other notebooks. I understand why the surface types are not the same between daily and monthly. I think they are take into account a sea_ice, so the ocean surface changes following the season. As we work over the tropics, we are not affected by this. But, in order to be consistent, we can use the monthly average of the surface types.

I have another request for you :

Can you add to the download files the sea level pressure standardized data (lower array) available at this link https://www.cpc.ncep.noaa.gov/data/indices/soi for the years from 2003 to 2017 please ?

Thank you very much

[malmans2]

Here is the code to retrieve those tables:

import pandas as pd

year_stop = 2025

kwargs = {
    "filepath_or_buffer": "https://www.cpc.ncep.noaa.gov/data/indices/soi",
    "sep": r"\s+",
    "nrows": year_stop - 1950,
    "index_col": "YEAR",
}
df_slp_anomaly = pd.read_csv(skiprows=3, **kwargs)
df_slp_standardized = pd.read_csv(skiprows=87, **kwargs)