<jupyter_start><jupyter_text>Accessing non-SNOTEL met site data from the SnowEx period  imports<jupyter_code>%%capture
!pip install git+https://github.com/friedrichknuth/driveanon.git
    
import driveanon
import pandas as pd
import matplotlib.pyplot as plt
import datetime<jupyter_output><empty_output><jupyter_text>load in datasets, format, and append to single dataframe<jupyter_code>#for accessing google drive data, blob ids and downloading
blob_ids = ['1Ioh1huZZpIJLRtoyGvWZfvJkaSQh0-QS', #to get site locations
            '1YT9bNveID55Lo5j4upMGrA7vEl5FWXp_', #Mesa West .csv
            '1v2zUMtjBSohk0OX6WCirZYkddysLiZY7', #GMSP .csv
            '1cwIsJQ5kCVcHpw6pRgywf6TKsmeg2RZC' #Mesa East .csv
            ] 
            #'', 
            #'' ] #Mesa Middle .csv
    
for blob_id in blob_ids:
    f = driveanon.save(blob_id)
#load in csv files and trim to relevant variables
files = ['MW-V7_clean.csv','GMSP-V2.csv','ME-V3.csv']
cols = ['station_id','Datetime_UTC','latitude','longitude','AirTC_10ft_Avg_C','IRtarget_nadir2_C_Avg','SnowDepth_m',
        'TCstringC_s50cm_Avg','TCstringC_s10cm_Avg','TCstringC_0cm_Avg',
        'TCstringC_10cm_Avg','TCstringC_50cm_Avg','TCstringC_100cm_Avg',
        'TCstringC_200cm_Avg'];

#date range for trimming
start_date = datetime.datetime(2020,1,1,0,0)
end_date = datetime.datetime(2020,3,31,23,59)

DF = pd.DataFrame(columns=cols)
for file in files:
    station = file[:-4]
    
    #hard coded this in cause I ran out of patience
    if station[:2] == 'GM':
        lat = 39.05084
        lon = -108.06144
    elif station[:2] == 'MM':
        lat = 39.03954
        lon = -107.94174
    elif station[:2] == 'ME':
        lat = 39.10358
        lon = -107.88383
    elif station[:2] == 'MW':
        lat = 39.03388
        lon = -108.21399
        
    df = pd.read_csv('./'+file)
    
    #converts any NAN text to NaN and numeric columns (fixes error raised)
    for (columnName, columnData) in df.iteritems():
        df[columnName] = pd.to_numeric(df[columnName],errors='coerce')
    
    #create and add site column
    df.insert(0,'station_id',station)
    
    #add 
    df.insert(1, 'latitude', lat)
    df.insert(2, 'longitude', lon)
    
    #create and add datetime column
    full_dates = pd.to_datetime(dict(year=df['Year'], month=df['Month'], day=df['Day'], hour = df['Hour'], minute = df['Minute']))
    df.insert(3, 'Datetime_UTC', full_dates)
    
    #trim to study period (January - March 2020)
    mask = (df['Datetime_UTC'] > start_date) & (df['Datetime_UTC'] <= end_date)
    df = df.loc[mask]
    
    #trim to necessary columns
    df = df[cols]
    
    #append data frame
    DF = DF.append(df, ignore_index=True)
    
DF.shape
DF.head()
DF.nunique()
DF['station_id'].value_counts()<jupyter_output><empty_output><jupyter_text>Use DF dataframe to add sites to plot (3 unique sites), and add to ERA comparison plots<jupyter_code>## old method ##
#colors = {'MW-V7_clean':1, 'GMSP-V2.csv':2, 'ME-V3':3} 
#plt.scatter(DF['Datetime_UTC'],DF['AirTC_10ft_Avg_C'],c=DF['station_id'].map(colors),s=5,marker='x',cmap='brg')
variable = 'TCstringC_200cm_Avg'

MW = DF[DF['station_id'] == 'MW-V7_clean']
GMSP = DF[DF['station_id'] == 'GMSP-V2']
ME = DF[DF['station_id'] == 'ME-V3']

plt.figure(figsize=(15,5))
plt.plot(MW['Datetime_UTC'], MW[variable], color='b',linestyle = ':',label='Mesa West ' + variable)
plt.plot(GMSP['Datetime_UTC'], GMSP[variable], color='g',linestyle = ':',label='GMSP-Skyway ' + variable)
plt.plot(ME['Datetime_UTC'], ME[variable], color='r',linestyle = ':',label='Mesa East ' + variable)
plt.ylabel('Temperature (C)')

#to trim by date
d1 = datetime.datetime(year=2020,month=1,day=31)
d2 = datetime.datetime(year=2020,month=2,day=29)
delta = d2 - d1

plt.xlim(d1,d2)
plt.grid(color='k', linestyle='-', linewidth=1)
plt.legend()
plt.savefig('./images/mesa_sites_T2m.jpg')
plt.show()
#read in PitSurfaceTemperatures.csv (see QueryPractice to show how to save these as a .csv)
f = 'PitSurfaceTemps.csv'
pitMax_df = pd.read_csv('./'+f)
pitMax_df.head()
#add datetime column
pit_datestr = pitMax_df['date'] + ' ' + pitMax_df['time']
pit_dates =  pd.to_datetime(pitMax_df['date'] + ' ' + pitMax_df['time'])#, format='%d%b%Y:%H:%M:%S.%f')
pitMax_df.insert(4, 'Datetime_UTC', pit_dates.dt.tz_localize(None)+datetime.timedelta(hours=6))
#remove timezone from datetime

#pitMax_df['Datetime_UTC'] = pitMax_df['Datetime_UTC'].replace(tzinfo=None)
pitMax_df.head()<jupyter_output><empty_output><jupyter_text>Get ERA5 2m air temperature data<jupyter_code>!pip install -q cdsapi
import cdsapi
c = cdsapi.Client()
c.retrieve(
    'reanalysis-era5-land',
    {
        'format': 'netcdf',
        'variable': [
            '2m_temperature', 'snow_albedo', 'snow_density',
            'snow_depth', 'snow_depth_water_equivalent', 'temperature_of_snow_layer',
            'total_precipitation',
        ],
        'year': '2020',
        'month': [
            '01', '02', '03',
        ],
        'day': [
            '01', '02', '03',
            '04', '05', '06',
            '07', '08', '09',
            '10', '11', '12',
            '13', '14', '15',
            '16', '17', '18',
            '19', '20', '21',
            '22', '23', '24',
            '25', '26', '27',
            '28', '29', '30',
            '31',
        ],
        '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': [
            39.2, -108.4, 38.8,
            -107.7,
        ],
    },
    'data/era5_2020-02-01_2020-02-29.nc')
import xarray as xr
import numpy as np
ds = xr.open_dataset('data/era5_2020-02-01_2020-02-29.nc')
print(ds)
ds['t2m'].values = ds['t2m'].values - 273.15
ds['t2m'].attrs['units'] = 'C'
ds['t2m'].attrs['long_name'] = 'Temperature (2 m above surface)'
ds['t2m'].mean(dim=('latitude', 'longitude'))
ds['t2m'].time
T = [t - np.timedelta64(0, 'h') for t in ds['time'].values] #convert to correct time
# SAME PLOT AS ABOVE BUT WITH SNOW PIT SURFACE TEMPS
variable = 'TCstringC_200cm_Avg'

MW = DF[DF['station_id'] == 'MW-V7_clean']
GMSP = DF[DF['station_id'] == 'GMSP-V2']
ME = DF[DF['station_id'] == 'ME-V3']

plt.figure(figsize=(15,5))
plt.plot(MW['Datetime_UTC'], MW[variable], color='b',linestyle = '--',label='Mesa West ' + variable)
plt.plot(GMSP['Datetime_UTC'], GMSP[variable], color='g',linestyle = '--',label='GMSP-Skyway ' + variable)
plt.plot(ME['Datetime_UTC'], ME[variable], color='r',linestyle = '--',label='Mesa East ' + variable)
#ERA5
plt.plot(T,ds['t2m'].mean(dim=('latitude', 'longitude')),color='k',linestyle = '-',label='ERA5 2m temperature')
plt.ylabel('Temperature (C)')
plt.scatter(pitMax_df['Datetime_UTC'],pitMax_df['temperature_C'],c='k',marker='x',label='pit snow surface obs.')
plt.xticks(rotation=45)
plt.title('Pits, Met sites, and ERA5 on the Mesa')

#to trim by date
d1 = datetime.datetime(year=2020,month=1,day=28)
d2 = datetime.datetime(year=2020,month=2,day=12)
delta = d2 - d1

plt.xlim(d1,d2)
plt.grid(color='k', linestyle='-', linewidth=1)
plt.legend()
plt.savefig('./images/mesa_sites_T2m_ERA5.jpg')
plt.show()<jupyter_output><empty_output><jupyter_text>Plotting sites on a map (from geospatial cookbook tutorial)<jupyter_code>import os
import requests
import rasterio as rio
import pandas as pd
from rasterio import plot, mask
import geopandas as gpd
from shapely.geometry import Point
import contextily as ctx
gm_bounds = [-108.34115668,   38.82320553, -107.72839859,   39.19563035]
#List of all products hosted by OpenTopography GlobalDEM API
demtype_list = ["SRTMGL3", "SRTMGL1", "SRTMGL1_E", "AW3D30", "AW3D30_E", "SRTM15Plus"]
#1 arcsec
demtype = "SRTMGL1_E"
base_url="https://portal.opentopography.org/API/globaldem?demtype={}&west={}&south={}&east={}&north={}&outputFormat=GTiff"
base_url.format(demtype, *gm_bounds)
#Paste API key below (from user account)
api_key = '<KEY>'
if api_key:
    demtype_list.extend(['NASADEM', 'COP30', 'COP90'])
    demtype = 'COP30'
def get_OT_GlobalDEM(demtype, bounds, out_fn=None, api_key=None):
    if out_fn is None:
        out_fn = '{}.tif'.format(demtype)
    
    if not os.path.exists(out_fn):
        #Prepare API request url
        #Bounds should be [minlon, minlat, maxlon, maxlat]
        url = base_url.format(demtype, *bounds)
        if api_key is not None:
            url = f'{url}&API_Key={api_key}'
        print(url)
        #Get
        response = requests.get(url)
        #Should check for 200
        #Write to disk
        open(out_fn, 'wb').write(response.content)
out_fn = f"GM_{demtype}.tif"
out_fn
get_OT_GlobalDEM(demtype, gm_bounds, out_fn, api_key)
!ls -lh $out_fn
#convert met site data to a geodataframe for plotting with DEM basemap
met_sites = pd.read_csv('./MesaMetSites.csv')
#met_sites['geometry'] = [Point(row['longitude'], row['latitude']) for row in met_sites]
met_sites_gdf = gpd.GeoDataFrame(met_sites, geometry=gpd.points_from_xy(met_sites.longitude, met_sites.latitude))
met_sites_gdf.head()
met_sites_gdf.crs = 4236
met_sites_gdf.crs
with rio.open(out_fn) as src:
    rio.plot.show(src)
import hvplot.pandas
from geoviews import tile_sources as gvts
#met_sites_gdf.hvplot(hover_cols=['index','name'])
map_tiles = gvts.EsriImagery
map_tiles * met_sites_gdf.to_crs('EPSG:3857').hvplot(hover_cols=['index','siteName'])<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1707827542904981, 'fraction_numerical': 0.05473001255755546, 'mean_word_length': 3.8985135827780626, 'pattern_counts': {'":': 0, '<': 23, '<?xml version=': 0, '>': 23, 'https://': 2, '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': '5', 'dir': 'jupyter-structured-clean-dedup', 'id': '34221', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 4068, 'n_tokens_neox': 3698, 'n_words': 721, 'path': 'contributors/Jeremy_Johnston/MetSite_ClipExtract.ipynb', 'repo_name': 'GiuliaMazzotti/hot-pow', 'size': '9583'}