gfs.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import ee\n",
    "import geemap\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import os\n",
    "#from paths import *\n",
    "import requests\n",
    "import pandas as pd\n",
    "import xarray as xr\n",
    "from os import listdir\n",
    "from datetime import datetime, timedelta, date\n",
    "import contextlib"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Initialize the Earth Engine module\n",
    "ee.Initialize()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#########################################################################\n",
    "############################ USER INPUTS ################################\n",
    "#########################################################################\n",
    "# PATHS\n",
    "# path to temporary folder to store tif files from gee\n",
    "TIFpath = '/nfs/depot/cce_u1/hill/dfh/op_snowmodel/get_met_data/GEE_Downloads_or_test/'\n",
    "TIFpath2 = '/nfs/depot/cce_u1/hill/dfh/op_snowmodel/get_met_data/GEE_Downloads_or_gfs/'\n",
    "# path to where you want your output met .dat fime\n",
    "OUTpath = '/nfs/depot/cce_u1/hill/dfh/op_snowmodel/or_snowmodel/met/mm_or_test.dat'\n",
    "\n",
    "# DOMAIN\n",
    "# choose the modeling domain\n",
    "domain = 'OR'\n",
    "\n",
    "# TIME\n",
    "# choose if want to set 'manual' or 'auto' date \n",
    "date_flag = 'auto'\n",
    "# If you choose 'manual' set your dates below  \n",
    "# This will start on the 'begin' date at 0:00 and the last iteration will \n",
    "# be on the day before the 'end' date below.\n",
    "st_dt = '2023-02-10'\n",
    "ed_dt = '2023-02-15'\n",
    "#########################################################################"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Date setup function\n",
    "def set_dates(st_dt,ed_dt,date_flag):\n",
    "    if date_flag == 'auto':\n",
    "        # ###automatically select date based on today's date \n",
    "        hoy = date.today()\n",
    "        #in line below change from 2 to 1? 11/23/2022\n",
    "        antes = timedelta(days = 1)\n",
    "        #end date 3 days before today's date\n",
    "        fecha = hoy - antes\n",
    "        eddt = fecha.strftime(\"%Y-%m-%d\") \n",
    "        #whole water year\n",
    "        #change the 3 to 2? 11/23/2022\n",
    "        if (hoy.month == 10) & (hoy.day == 2):\n",
    "            eddt = fecha.strftime(\"%Y-%m-%d\") \n",
    "            stdt = str(hoy.year - 1)+'-10-01'\n",
    "        #start dates\n",
    "        elif fecha.month <10:\n",
    "            stdt = str(fecha.year - 1)+'-10-01'\n",
    "            stdt='2023-02-17' #temp addition (Feb 16 2023). Delete when done testing\n",
    "        else:\n",
    "            stdt = str(fecha.year)+'-10-01'\n",
    "            stdt='2023-02-17' #temp addition (Feb 16 2023). Delete when done testing\n",
    "    elif date_flag == 'manual':\n",
    "        stdt = st_dt\n",
    "        # add one day to end date because GEE ends on date before last date\n",
    "        eddt = (datetime.strptime(ed_dt, \"%Y-%m-%d\")+timedelta(days = 1)).strftime(\"%Y-%m-%d\") \n",
    "    return stdt, eddt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Download CFSv2 met data function\n",
    "def get_cfsv2(domain, TIFpath, stdt, eddt):\n",
    "    # in GEE the last iteration is on the day before the 'end' date below\n",
    "    # we adjust this here since it is not intuative\n",
    "    #eddt = (datetime.strptime(eddt, '%Y-%m-%d')+timedelta(days = 1)).strftime('%Y-%m-%d')\n",
    "    \n",
    "    #create directory with initiation date for ensemble if it doesn't exist\n",
    "    get_ipython().system('mkdir -p $TIFpath')\n",
    "\n",
    "    #path to CSO domains\n",
    "    domains_resp = requests.get(\"https://raw.githubusercontent.com/snowmodel-tools/preprocess_python/master/CSO_domains.json\")\n",
    "    domains = domains_resp.json()\n",
    "\n",
    "    minLat = domains[domain]['Bbox']['latmin']\n",
    "    #// Input the minimum long, lower left corner\n",
    "    minLong = domains[domain]['Bbox']['lonmin']\n",
    "    #// Input the max lat, upper right corner\n",
    "    maxLat = domains[domain]['Bbox']['latmax']\n",
    "    #// Input the max Long, upper right corner\n",
    "    maxLong = domains[domain]['Bbox']['lonmax']\n",
    "\n",
    "    #/ These are the min and max corners of your reanalysis in Lat, Long (create a slightly larger box)\n",
    "    #// Input the minimum lat, lower left corner\n",
    "    minLatMET = (minLat - 0.25);\n",
    "    #// print(minLat2);\n",
    "    #// Input the minimum long, lower left corner\n",
    "    minLongMET = (minLong - 0.5);\n",
    "    #// Input the max lat, upper right corner\n",
    "    maxLatMET = (maxLat + 0.25);\n",
    "    #// Input the max Long, upper right corner\n",
    "    maxLongMET = (maxLong + 0.5);\n",
    "\n",
    "    '''// Define the final output projection using EPSG codes'''\n",
    "    epsg_code = domains[domain]['mod_proj']\n",
    "\n",
    "    my_domain = ee.Geometry.Rectangle(**{'coords':[minLong,minLat,maxLong,maxLat],'proj': 'EPSG:4326','geodesic':True,});\n",
    "    my_domain_met = ee.Geometry.Rectangle([minLongMET,minLatMET,maxLongMET,maxLatMET])\n",
    "    \n",
    "    # download reanalysis data\n",
    "    cfsv2 = ee.ImageCollection('NOAA/CFSV2/FOR6H')         .filterBounds(my_domain_met)         .filter(ee.Filter.date(stdt,eddt))\n",
    "    data = cfsv2.select('Temperature_height_above_ground',         'Geopotential_height_surface',         'u-component_of_wind_height_above_ground',         'v-component_of_wind_height_above_ground',         'Pressure_surface',         'Specific_humidity_height_above_ground',         'Precipitation_rate_surface_6_Hour_Average')\n",
    "    with contextlib.redirect_stdout(None):\n",
    "        geemap.ee_export_image_collection(data, out_dir=TIFpath,region=my_domain_met,scale=22200,crs=epsg_code)\n",
    "    #print(data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Download GFS met data function\n",
    "def get_gfs(domain, TIFpath2, stdt2, eddt2):\n",
    "    \n",
    "    #create directory with initiation date for ensemble if it doesn't exist\n",
    "    get_ipython().system('mkdir -p $TIFpath2')\n",
    "\n",
    "    #path to CSO domains\n",
    "    domains_resp = requests.get(\"https://raw.githubusercontent.com/snowmodel-tools/preprocess_python/master/CSO_domains.json\")\n",
    "    domains = domains_resp.json()\n",
    "\n",
    "    minLat = domains[domain]['Bbox']['latmin']\n",
    "    #// Input the minimum long, lower left corner\n",
    "    minLong = domains[domain]['Bbox']['lonmin']\n",
    "    #// Input the max lat, upper right corner\n",
    "    maxLat = domains[domain]['Bbox']['latmax']\n",
    "    #// Input the max Long, upper right corner\n",
    "    maxLong = domains[domain]['Bbox']['lonmax']\n",
    "\n",
    "    #/ These are the min and max corners of your reanalysis in Lat, Long (create a slightly larger box)\n",
    "    #// Input the minimum lat, lower left corner\n",
    "    minLatMET = (minLat - 0.25);\n",
    "    #// print(minLat2);\n",
    "    #// Input the minimum long, lower left corner\n",
    "    minLongMET = (minLong - 0.5);\n",
    "    #// Input the max lat, upper right corner\n",
    "    maxLatMET = (maxLat + 0.25);\n",
    "    #// Input the max Long, upper right corner\n",
    "    maxLongMET = (maxLong + 0.5);\n",
    "\n",
    "    '''// Define the final output projection using EPSG codes'''\n",
    "    epsg_code = domains[domain]['mod_proj']\n",
    "\n",
    "    my_domain = ee.Geometry.Rectangle(**{'coords':[minLong,minLat,maxLong,maxLat],'proj': 'EPSG:4326','geodesic':True,});\n",
    "    my_domain_met = ee.Geometry.Rectangle([minLongMET,minLatMET,maxLongMET,maxLatMET])\n",
    "    \n",
    "    # download reanalysis data\n",
    "    print(stdt2)\n",
    "    print(eddt2)\n",
    "    gfs = ee.ImageCollection('NOAA/GFS0P25').filterBounds(my_domain_met).filter(ee.Filter.date(stdt2,eddt2)).filter(ee.Filter.Or(ee.Filter.eq('forecast_hours',6), ee.Filter.eq('forecast_hours',12),ee.Filter.eq('forecast_hours',18),ee.Filter.eq('forecast_hours',24), ee.Filter.eq('forecast_hours',30),ee.Filter.eq('forecast_hours',36),ee.Filter.eq('forecast_hours',42),ee.Filter.eq('forecast_hours',48),ee.Filter.eq('forecast_hours',54),ee.Filter.eq('forecast_hours',60),ee.Filter.eq('forecast_hours',66),ee.Filter.eq('forecast_hours',72)))\n",
    "    data2 = gfs.select('temperature_2m_above_ground', 'u_component_of_wind_10m_above_ground', 'v_component_of_wind_10m_above_ground',     'relative_humidity_2m_above_ground',         'total_precipitation_surface')\n",
    "    with contextlib.redirect_stdout(None):\n",
    "        geemap.ee_export_image_collection(data2, out_dir=TIFpath2,region=my_domain_met,scale=22200,crs=epsg_code)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# function to check for missing dates (for the CFS data)\n",
    "def missing_slice_check(stdt, eddt, TIFpath):\n",
    "    # create a 6-hourly timeseries with no missing values from the start to end date\n",
    "    timesin = pd.date_range(start=stdt, end=eddt, freq='6H')[:-1]\n",
    "    for time in timesin:\n",
    "        nam = time.strftime('%Y%m%d%H')\n",
    "        \n",
    "    # compile list of all tif files downloaded from gee\n",
    "    gee_times =[]\n",
    "\n",
    "    for file in listdir(TIFpath):\n",
    "        if file.endswith(\"tif\"):\n",
    "            datetmp = datetime.strptime(file[:-4], '%Y%m%d%H')\n",
    "            gee_times.append(datetmp)\n",
    "    gee_times = sorted(gee_times)\n",
    "\n",
    "    # check for to see if all time slices downloaded from GEE\n",
    "    if len(timesin) != len(gee_times):\n",
    "    #### on 4/16 Nina edited code to print all missing timeslices\n",
    "        print('gee is missing timeslice(s):\\n',timesin[~timesin.isin(gee_times)].values)\n",
    "        \n",
    "        # if 4 or more consecutive timeslices are missing - quit the function\n",
    "        duration = []\n",
    "        for i in range(len(gee_times)-1):\n",
    "            time_delta = gee_times[i+1] - gee_times[i]\n",
    "            duration.append(time_delta.total_seconds()/60/60)\n",
    "        if max(duration) >= 48:    \n",
    "            print('at least two full days of met data are missing - quitting function')\n",
    "\n",
    "        # if there are less than 4 missing consecutive time slices \n",
    "        # repeat the met data from the last valid time slice \n",
    "        else:\t\n",
    "            missing_idx = np.where(~timesin.isin(gee_times))\t\n",
    "            missing_dt = timesin[missing_idx]\t\n",
    "            if len(missing_dt)==1:\t\n",
    "                if missing_idx == 0:\t\n",
    "                    pre_dt=TIFpath+timesin[np.squeeze(missing_idx)+1].strftime('%Y%m%d%H')+'.tif'\t\n",
    "                    mis_dt = TIFpath+timesin[np.squeeze(missing_idx)].strftime('%Y%m%d%H')+'.tif' \t\n",
    "                    get_ipython().system('cp $pre_dt $mis_dt')\t\n",
    "                    print('replaced', timesin[np.squeeze(missing_idx)].strftime('%Y%m%d%H'),' with ', timesin[np.squeeze(missing_idx)-1].strftime('%Y%m%d%H'))\t\n",
    "                else:\t\n",
    "                    pre_dt=TIFpath+timesin[np.squeeze(missing_idx)-1].strftime('%Y%m%d%H')+'.tif'\t\n",
    "                    mis_dt = TIFpath+timesin[np.squeeze(missing_idx)].strftime('%Y%m%d%H')+'.tif' \t\n",
    "                    get_ipython().system('cp $pre_dt $mis_dt')\t\n",
    "                    print('replaced', timesin[np.squeeze(missing_idx)].strftime('%Y%m%d%H'),' with ', timesin[np.squeeze(missing_idx)-1].strftime('%Y%m%d%H'))\t\n",
    "            else:\t\n",
    "                for j in range(len(missing_dt)):\t\n",
    "                    if np.squeeze(missing_idx)[j] == 0:\t\n",
    "                        print('choose earlier start date so missing time slices can be filled in')\t\n",
    "                    else:\t\n",
    "                        pre_dt=TIFpath+timesin[np.squeeze(missing_idx)[j]-1].strftime('%Y%m%d%H')+'.tif'\t\n",
    "                        mis_dt = TIFpath+timesin[np.squeeze(missing_idx)[j]].strftime('%Y%m%d%H')+'.tif' \t\n",
    "                        get_ipython().system('cp $pre_dt $mis_dt')\t\n",
    "                        print('replaced', timesin[np.squeeze(missing_idx)[j]].strftime('%Y%m%d%H'),' with ', timesin[np.squeeze(missing_idx)[j]-1].strftime('%Y%m%d%H'))\t\n",
    "                        "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Format gee files for SnowModel function\n",
    "# this function will take care of the CFS data...there will be a second function that will try to append the GFS data...Ha. Good luck.\n",
    "def MET2SM(TIFpath, OUTpath, stdt, eddt):\n",
    "    \n",
    "    # create a 6-hourly timeseries with no missing values from the start to end date\n",
    "    timesin = pd.date_range(start=stdt, end=eddt, freq='6H')[:-1]\n",
    "    \n",
    "    #load first tif to get dimensions\n",
    "    ar = xr.open_rasterio(TIFpath+timesin[0].strftime('%Y%m%d%H')+'.tif')\n",
    "    #print(TIFpath+timesin[0].strftime('%Y%m%d%H')+'.tif')\n",
    "    \n",
    "    # empty arrays for each met variable\n",
    "    T = np.empty((len(timesin),ar.shape[1],ar.shape[2]))\n",
    "    Z = np.empty((len(timesin),ar.shape[1],ar.shape[2]))\n",
    "    U = np.empty((len(timesin),ar.shape[1],ar.shape[2]))\n",
    "    V = np.empty((len(timesin),ar.shape[1],ar.shape[2]))\n",
    "    P = np.empty((len(timesin),ar.shape[1],ar.shape[2]))\n",
    "    H = np.empty((len(timesin),ar.shape[1],ar.shape[2]))\n",
    "    PR = np.empty((len(timesin),ar.shape[1],ar.shape[2]))\n",
    "\n",
    "    # extract met data from tifs \n",
    "    for i in range(len(timesin)):\n",
    "\n",
    "        #load tif file\n",
    "        nam = TIFpath+timesin[i].strftime('%Y%m%d%H')+'.tif'\n",
    "        ar = xr.open_rasterio(nam)\n",
    "        T[i,:,:] = ar[0,:,:]\n",
    "        Z[i,:,:] = ar[1,:,:]\n",
    "        U[i,:,:] = ar[2,:,:]\n",
    "        V[i,:,:] = ar[3,:,:]\n",
    "        P[i,:,:] = ar[4,:,:]\n",
    "        H[i,:,:] = ar[5,:,:]\n",
    "        PR[i,:,:] = ar[6,:,:]\n",
    "\n",
    "    #number of timesteps per dat \n",
    "    pointsperday = 4\n",
    "\n",
    "    #compute number of grid points and time steps from size of 3d matrix\n",
    "    t,y,x=PR.shape\n",
    "    gridpts=x*y\n",
    "    tsteps=t\n",
    "\n",
    "    #create y m d h vectors\n",
    "    year = timesin.year\n",
    "    month = timesin.month\n",
    "    day = timesin.day\n",
    "    hour = timesin.hour\n",
    "\n",
    "    #create ID numbers for the grid points\n",
    "    ID=1000000+np.linspace(1,gridpts,gridpts)\n",
    "\n",
    "    #create matrices of x and y values\n",
    "    X, Y = np.meshgrid(ar.x.values, ar.y.values)\n",
    "    X=X.flatten(order='F')\n",
    "    Y=Y.flatten(order='F')\n",
    "\n",
    "    #elevation is static (doesn't change with time)\n",
    "    elev=Z[1,:,:].flatten(order='F')\n",
    "\n",
    "    #find number of grid points with <0 elevation. Note: this is related to the\n",
    "    #subroutine met_data_check in the preprocess_code.f. that subroutine seems\n",
    "    #to suggest that negative elevations are ok (say, death valley). But, the\n",
    "    #code itself checks for negative elevations and stops execution is any\n",
    "    #negatives are found.\n",
    "    I = np.where(elev>=0)\n",
    "    validgridpts=np.shape(I)[1]\n",
    "\n",
    "    #remove data at points with neg elevations\n",
    "    ID=ID[I]\n",
    "    X=X[I]\n",
    "    Y=Y[I]\n",
    "    elev=elev[I]\n",
    "\n",
    "    #we are now ready to begin our main loop over the time steps.\n",
    "    fid= open(OUTpath,\"w+\")\n",
    "\n",
    "    for j in range(tsteps):\n",
    "        #first we write the number of grid points\n",
    "        fid.write('{0:6d}\\n'.format(validgridpts))\n",
    "\n",
    "        #prep data matrix for this time step. First, grab the jth time slice\n",
    "        Prtmp=PR[j,:,:].flatten(order='F')\n",
    "        Htmp=H[j,:,:].flatten(order='F')\n",
    "        Ptmp=P[j,:,:].flatten(order='F')\n",
    "        Ttmp=T[j,:,:].flatten(order='F')\n",
    "        Utmp=U[j,:,:].flatten(order='F')\n",
    "        Vtmp=V[j,:,:].flatten(order='F')\n",
    "\n",
    "        #remove data at points with neg elevations\n",
    "        Prtmp=Prtmp[I]\n",
    "        Htmp=Htmp[I]\n",
    "        Ptmp=Ptmp[I]\n",
    "        Ttmp=Ttmp[I]\n",
    "        Utmp=Utmp[I]\n",
    "        Vtmp=Vtmp[I]\n",
    "\n",
    "        #convert precip rate to precip DEPTH (mm) during time interval\n",
    "        Prtmp=Prtmp*24*3600/pointsperday\n",
    "\n",
    "        #convert specific hum. to RH from Clausius-Clapeyron. T is still in K\n",
    "        RHtmp=0.263*Ptmp*Htmp*(np.exp(17.67*(Ttmp-273.16)/(Ttmp-29.65)))**(-1)\n",
    "\n",
    "        #compute wind speed\n",
    "        SPDtmp=np.sqrt(Utmp**2+Vtmp**2)\n",
    "\n",
    "        #compute wind direction. 0-360, with 0 being true north! 90 east, etc.\n",
    "        DIRtmp=np.degrees(np.arctan2(Utmp,Vtmp))\n",
    "        K=np.where(DIRtmp>=180)\n",
    "        J=np.where(DIRtmp<180)\n",
    "        DIRtmp[K]=DIRtmp[K]-180\n",
    "        DIRtmp[J]=DIRtmp[J]+180\n",
    "\n",
    "        #put T in C\n",
    "        Ttmp=Ttmp-273.16\n",
    "\n",
    "        for z in range(len(Prtmp)): \n",
    "\n",
    "            fid.write('{:5.0f}\\t'.format(int(year[j]))+'{:5.0f}\\t'.format(int(month[j]))+\n",
    "                      '{:3.0f}\\t'.format(int(day[j]))+'{:6.3f}\\t'.format(hour[j])+\n",
    "                      '{:9.0f}\\t'.format(int(ID[z]))+'{:12.1f}\\t'.format(X[z])+\n",
    "                      '{:12.1f}\\t'.format(Y[z])+'{:8.1f}\\t'.format(elev[z])+\n",
    "                      '{:9.2f}\\t'.format(Ttmp[z])+'{:9.2f}\\t'.format(RHtmp[z])+\n",
    "                      '{:9.2f}\\t'.format(SPDtmp[z])+'{:9.2f}\\t'.format(DIRtmp[z])+\n",
    "                      '{:9.2f}\\n'.format(Prtmp[z]))\n",
    "    fid.close()\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Format gee files for SnowModel function\n",
    "# this function will take care of the GFS data.\n",
    "def MET2SM2(TIFpath, TIFpath2, OUTpath, stdt, eddt):\n",
    "    \n",
    "    eddt_tmp=(datetime.strptime(eddt,'%Y-%m-%d')+timedelta(days=3)).strftime('%Y-%m-%d')\n",
    "    \n",
    "    # create a 6-hourly timeseries with no missing values from the start to end date of the GFS data\n",
    "    timesin = pd.date_range(start=eddt, end=eddt_tmp, freq='6H')[:-1]\n",
    "    \n",
    "    #dump listing of tifs (GFS) into array. Critical there are no extra files in there! \n",
    "    #there is probably a better way, but this seems the easiest for now...\n",
    "    x = sorted(os.listdir(TIFpath2))\n",
    "    \n",
    "    #load first tif to get dimensions\n",
    "    ar = xr.open_rasterio(TIFpath2+x[0])\n",
    "    \n",
    "    # empty arrays for each met variable\n",
    "    T = np.empty((len(timesin),ar.shape[1],ar.shape[2]))\n",
    "    U = np.empty((len(timesin),ar.shape[1],ar.shape[2]))\n",
    "    V = np.empty((len(timesin),ar.shape[1],ar.shape[2]))\n",
    "    H = np.empty((len(timesin),ar.shape[1],ar.shape[2]))\n",
    "    PR = np.empty((len(timesin),ar.shape[1],ar.shape[2]))\n",
    "    \n",
    "    #note...Z data not in GFS, but GFS and CFS are same size. We will have to read Z from one of the CFS files, but \n",
    "    #go ahead and initialize size here.\n",
    "    Z = np.empty((len(timesin),ar.shape[1],ar.shape[2]))\n",
    "\n",
    "    # extract met data from tifs \n",
    "    for i in range(len(timesin)):\n",
    "\n",
    "        #load tif file\n",
    "        nam = TIFpath2+x[i]\n",
    "        ar = xr.open_rasterio(nam)\n",
    "        T[i,:,:] = ar[0,:,:]\n",
    "        U[i,:,:] = ar[1,:,:]\n",
    "        V[i,:,:] = ar[2,:,:]\n",
    "        H[i,:,:] = ar[3,:,:]\n",
    "        PR[i,:,:] = ar[4,:,:]\n",
    "\n",
    "    #number of timesteps per dat \n",
    "    pointsperday = 4\n",
    "\n",
    "    #compute number of grid points and time steps from size of 3d matrix\n",
    "    t,y,x=PR.shape\n",
    "    gridpts=x*y\n",
    "    tsteps=t\n",
    "\n",
    "    #create y m d h vectors\n",
    "    year = timesin.year\n",
    "    month = timesin.month\n",
    "    day = timesin.day\n",
    "    hour = timesin.hour\n",
    "\n",
    "    #create ID numbers for the grid points\n",
    "    ID=1000000+np.linspace(1,gridpts,gridpts)\n",
    "\n",
    "    #create matrices of x and y values\n",
    "    X, Y = np.meshgrid(ar.x.values, ar.y.values)\n",
    "    X=X.flatten(order='F')\n",
    "    Y=Y.flatten(order='F')\n",
    "    \n",
    "    #ok...so now let us deal with the Z data from one of the CFS files.\n",
    "    # create a 6-hourly timeseries with no missing values from the start to end date\n",
    "    timesin_tmp = pd.date_range(start=stdt, end=eddt, freq='6H')[:-1]\n",
    "    \n",
    "    #load first tif to get dimensions\n",
    "    ar2 = xr.open_rasterio(TIFpath+timesin_tmp[0].strftime('%Y%m%d%H')+'.tif')\n",
    "    print(TIFpath+timesin_tmp[0].strftime('%Y%m%d%H')+'.tif')\n",
    "    Z[0,:,:] = ar2[1,:,:]\n",
    "\n",
    "    #elevation is static (doesn't change with time)\n",
    "    elev=Z[0,:,:].flatten(order='F')\n",
    "\n",
    "    #find number of grid points with <0 elevation. Note: this is related to the\n",
    "    #subroutine met_data_check in the preprocess_code.f. that subroutine seems\n",
    "    #to suggest that negative elevations are ok (say, death valley). But, the\n",
    "    #code itself checks for negative elevations and stops execution is any\n",
    "    #negatives are found.\n",
    "    I = np.where(elev>=0)\n",
    "    validgridpts=np.shape(I)[1]\n",
    "\n",
    "    #remove data at points with neg elevations\n",
    "    ID=ID[I]\n",
    "    X=X[I]\n",
    "    Y=Y[I]\n",
    "    elev=elev[I]\n",
    "\n",
    "    #we are now ready to begin our main loop over the time steps. Note use of 'a' option\n",
    "    #since we are appending to the snowmodel met file generated for the CFS data.\n",
    "    fid= open(OUTpath,\"a\")\n",
    "\n",
    "    for j in range(tsteps):\n",
    "        #first we write the number of grid points\n",
    "        fid.write('{0:6d}\\n'.format(validgridpts))\n",
    "\n",
    "        #prep data matrix for this time step. First, grab the jth time slice\n",
    "        #note that precip is super confusing. The GEE data catalog documentation is pretty\n",
    "        #poor. A better reference: https://benny.istan.to/blog/20210416-how-to-get-daily-rainfall-forecast-data-from-gfs-part-2\n",
    "        #from that...it SEEMS, that if I just take the precip grids every 6 hours (which I have done),\n",
    "        #then the grids will be precip depth over the previous 6 hour interval...which is precisely what\n",
    "        #we want. It may take further investigation to confirm this.\n",
    "        Prtmp=PR[j,:,:].flatten(order='F')\n",
    "        Htmp=H[j,:,:].flatten(order='F')\n",
    "        Ttmp=T[j,:,:].flatten(order='F')\n",
    "        Utmp=U[j,:,:].flatten(order='F')\n",
    "        Vtmp=V[j,:,:].flatten(order='F')\n",
    "\n",
    "        #remove data at points with neg elevations\n",
    "        Prtmp=Prtmp[I] #gfs precip is already depth...no need to convert\n",
    "        Htmp=Htmp[I] #gfs humidity is relative already...no need to convert\n",
    "        Ttmp=Ttmp[I] #gfs temp is C already, not K, no need to convert...\n",
    "        Utmp=Utmp[I]\n",
    "        Vtmp=Vtmp[I]\n",
    "\n",
    "        #compute wind speed\n",
    "        SPDtmp=np.sqrt(Utmp**2+Vtmp**2)\n",
    "\n",
    "        #compute wind direction. 0-360, with 0 being true north! 90 east, etc.\n",
    "        DIRtmp=np.degrees(np.arctan2(Utmp,Vtmp))\n",
    "        K=np.where(DIRtmp>=180)\n",
    "        J=np.where(DIRtmp<180)\n",
    "        DIRtmp[K]=DIRtmp[K]-180\n",
    "        DIRtmp[J]=DIRtmp[J]+180\n",
    "\n",
    "        for z in range(len(Prtmp)): \n",
    "\n",
    "            fid.write('{:5.0f}\\t'.format(int(year[j]))+'{:5.0f}\\t'.format(int(month[j]))+\n",
    "                      '{:3.0f}\\t'.format(int(day[j]))+'{:6.3f}\\t'.format(hour[j])+\n",
    "                      '{:9.0f}\\t'.format(int(ID[z]))+'{:12.1f}\\t'.format(X[z])+\n",
    "                      '{:12.1f}\\t'.format(Y[z])+'{:8.1f}\\t'.format(elev[z])+\n",
    "                      '{:9.2f}\\t'.format(Ttmp[z])+'{:9.2f}\\t'.format(Htmp[z])+\n",
    "                      '{:9.2f}\\t'.format(SPDtmp[z])+'{:9.2f}\\t'.format(DIRtmp[z])+\n",
    "                      '{:9.2f}\\n'.format(Prtmp[z]))\n",
    "    fid.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# set time parameters for CFS2\n",
    "stdt, eddt = set_dates(st_dt,ed_dt,date_flag)\n",
    "print(stdt)\n",
    "print(eddt)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# set time parameters for GFS. This will grab three additional days of forecast. The timing is tricky\n",
    "# the cfs call grabs data that will end at hour 18 on the day before eddt. So, we wish to pick up with hour 0\n",
    "# on eddt itself. BUT, the GFS forecast does not issue a complete forecast at hour 0. So, we need to\n",
    "# go back to hour 18 of the day before eddt and grab GFS starting then. Whew. In addition to this, GFS forecasts\n",
    "# are every hour, and we only want things every six hours, to match CFS.\n",
    "stdt_tmp=(datetime.strptime(eddt,'%Y-%m-%d')+timedelta(days=-1)).strftime('%Y-%m-%d')\n",
    "eddt_tmp=(datetime.strptime(eddt,'%Y-%m-%d')+timedelta(days=0)).strftime('%Y-%m-%d')\n",
    "# append string to help narrow down \n",
    "stdt2=stdt_tmp + 'T18:00:00Z'\n",
    "eddt2=eddt_tmp + 'T00:00:00Z' #this will end up grabbing 384 hours of forecast (too much...trim later with filter)\n",
    "print(stdt2)\n",
    "print(eddt2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# download GEE data for CFS\n",
    "get_cfsv2(domain, TIFpath, stdt, eddt)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# download GEE data for GFS\n",
    "get_gfs(domain, TIFpath2, stdt2, eddt2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# fill in missing time slices or throw error if missing >4 slices. Currently, this only operates on the CFS reanalysis and not the GFS forecast.\n",
    "missing_slice_check(stdt, eddt, TIFpath)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# build SnowModel met file\n",
    "MET2SM(TIFpath, OUTpath, stdt, eddt)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# build SnowModel met file\n",
    "MET2SM2(TIFpath, TIFpath2, OUTpath, stdt, eddt)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# delete directory with tif files \n",
    "get_ipython().system('rm -rf $TIFpath')\n",
    "get_ipython().system('rm -rf $TIFpath2')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [conda env:ee]",
   "language": "python",
   "name": "conda-env-ee-py"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}