fetch_data.py

import os
import copernicusmarine
import numpy as np
from datetime import datetime, timedelta
import netCDF4 as nc
from scipy.interpolate import RegularGridInterpolator, NearestNDInterpolator
from scipy.ndimage import gaussian_filter
import matplotlib.pyplot as plt
from pathlib import Path
import shutil

BASE_DIR = Path(__file__).resolve().parent
REQUIRED_DIR = BASE_DIR / "required_data"
DAILY_DIR = BASE_DIR / "daily_data"
BACKUP_DIR = BASE_DIR / "daily_backup"

DAILY_DIR.mkdir(exist_ok=True)
BACKUP_DIR.mkdir(exist_ok=True)

def clear_directory(directory: Path) -> None:
    for item in directory.iterdir():
        if item.is_dir():
            shutil.rmtree(item)
        else:
            item.unlink()

existing_daily_files = list(DAILY_DIR.glob("*"))
if existing_daily_files:
    clear_directory(BACKUP_DIR)

    for src in existing_daily_files:
        dst = BACKUP_DIR / src.name
        if src.is_dir():
            shutil.copytree(src, dst)
        else:
            shutil.copy2(src, dst)

masks_red=np.load(REQUIRED_DIR / "masks_red.npy")
masks=np.load(REQUIRED_DIR / "masks.npy")
print(masks_red.shape)
latitude_red=np.load(REQUIRED_DIR / "latitude_red.npy")
longitude_red=np.load(REQUIRED_DIR / "longitude_red.npy")

yesterday = datetime.today() - timedelta(days=1)
formatted_date = yesterday.strftime('%Y-%m-%d')

depths=[0,10,30,60,100]

c=np.concatenate([x.reshape(latitude_red.shape[0],longitude_red.shape[0],1) for x in np.meshgrid(latitude_red,longitude_red,indexing="ij")],axis=2)

c=c.reshape(-1,2)


data={
    "Chla":["chl","cmems_mod_med_bgc-pft_anfc_4.2km_P1D-m"],
    "N1p":["po4","cmems_mod_med_bgc-nut_anfc_4.2km_P1D-m"],
    "N3n":["no3","cmems_mod_med_bgc-nut_anfc_4.2km_P1D-m"],
    "S":["so","cmems_mod_med_phy-sal_anfc_4.2km_P1D-m"],
    "T":["thetao","cmems_mod_med_phy-tem_anfc_4.2km_P1D-m"]

}

depths_nc_path = DAILY_DIR / "depths.nc"
if not depths_nc_path.exists():
    copernicusmarine.subset(
            dataset_id="cmems_mod_med_phy_anfc_4.2km_static",
            minimum_longitude=11,
            maximum_longitude=18,
            minimum_latitude=42,
            maximum_latitude=47,
            minimum_depth=0,
            maximum_depth=120,
            start_datetime=formatted_date,
            end_datetime=formatted_date,
            output_filename=str(depths_nc_path),
            force_download=True
        )

d_dataset = nc.Dataset(str(depths_nc_path)).variables["depth"][:]

vars_list =["Chla","N1p","N3n","S","T"]

for i, depth in enumerate(depths):
    for var in vars_list:
        depth_tmp = d_dataset[np.argmin(d_dataset - depth)]

        nc_filename = DAILY_DIR / f"{var}_{depth}.nc"
        npy_filename = DAILY_DIR / f"{var}_{depth}.npy"
        png_filename = DAILY_DIR / f"miao_{var}_{depth}.png"

        if not nc_filename.exists():
            copernicusmarine.subset(
                dataset_id=data[var][1],
                variables=[data[var][0]],
                minimum_longitude=11,
                maximum_longitude=18,
                minimum_latitude=42,
                maximum_latitude=47,
                start_datetime=formatted_date,
                end_datetime=formatted_date,
                minimum_depth=depth_tmp,
                maximum_depth=depth_tmp,
                output_filename=str(nc_filename),
                force_download=True,
            )

        a = nc.Dataset(str(nc_filename))
        longitude = a.variables["longitude"][:].data
        latitude = a.variables["latitude"][:].data
        tmp = a.variables[data[var][0]][0, 0].data
        a.close()

        tmp[tmp > 1e10] = np.nan

        rg = RegularGridInterpolator(
            [latitude, longitude],
            tmp,
            method="nearest",
            fill_value=np.nan,
        )
        miao = rg(c).reshape(latitude_red.shape[0], longitude_red.shape[0])

        valid_mask = ~np.isnan(miao)
        valid_pts = c[valid_mask.ravel()]
        valid_vals = miao[valid_mask]
        nn = NearestNDInterpolator(valid_pts, valid_vals)
        miao = nn(c).reshape(latitude_red.shape[0], longitude_red.shape[0])

        miao = gaussian_filter(miao, sigma=1.5, mode='nearest')
        miao = gaussian_filter(miao, sigma=1.5, mode='nearest')

        plt.imshow(np.ma.masked_array(miao, masks_red[i]))
        plt.savefig(png_filename)
        plt.close()

        miao_masked_flat = miao.reshape(-1)[np.logical_not(masks_red[i]).reshape(-1)]
        np.save(npy_filename, miao_masked_flat)

expected_files = []
expected_files.append(depths_nc_path)
for depth in depths:
    for var in vars_list:
        expected_files.append(DAILY_DIR / f"{var}_{depth}.nc")
        expected_files.append(DAILY_DIR / f"{var}_{depth}.npy")
        expected_files.append(DAILY_DIR / f"miao_{var}_{depth}.png")

for path in expected_files:
    if not path.exists():
        backup_path = BACKUP_DIR / path.name
        if backup_path.exists():
            shutil.copy2(backup_path, path)
            print(f"[RECOVER] Restored {path.name} from backup.")
        else:
            print(f"[WARNING] Missing {path.name} in both daily_data and daily_backup.")