Examples

Reading and plotting ASCAT data

This example program reads and plots ASCAT SSM and SWI data with different masking options. It can be found in the /bin folder of the pytesmo package under the name plot_ASCAT_data.py.

In[1]:

import pytesmo.io.sat.ascat as ascat
import os
import matplotlib.pyplot as plt

In[2]:

#I've downloaded my ASCAT data to a folder on my D drive
path_to_ascat_ssm_data = os.path.join('D:\\','small_projects','cpa_2013_07_userformat_reader',
                                      'data','ASCAT_SSM_25km_ts_WARP5.5_R0.1','data')
path_to_ascat_swi_data = os.path.join('D:\\','small_projects','cpa_2013_07_userformat_reader',
                                      'data','ASCAT_SWI_25km_ts_WARP5.5_R0.1','data')
#path to grid definition file, default name TUW_W54_01_lonlat-ld-land.txt
path_to_grid_definition = os.path.join('D:\\','small_projects','cpa_2013_07_userformat_reader',
                                      'data','auxiliary_data','grid_info')
#path to advisory flags from FTP Server
path_to_adv_flags = os.path.join('D:\\','small_projects','cpa_2013_07_userformat_reader',
                                      'data','auxiliary_data','advisory_flags')

In[3]:

#init the ASCAT_SSM reader with the paths
ascat_SSM_reader = ascat.Ascat_SSM(path_to_ascat_ssm_data,path_to_grid_definition,
                                   advisory_flags_path = path_to_adv_flags)

In[4]:

lon, lat = 16, 48

#reads ssm data nearest to this lon,lat coordinates
ssm_data_raw = ascat_SSM_reader.read_ssm(lon,lat)

#plot the data using pandas builtin plot functionality
ssm_data_raw.plot()
plt.show()
_images/plot_ascat_data_3_0.png

In[5]:

#read the same data but mask observations where the SSF shows frozen
#and where frozen and snow probabilty are greater than 20%
ssm_data_masked = ascat_SSM_reader.read_ssm(lon,lat,mask_ssf=True,mask_frozen_prob=20,mask_snow_prob=20)

#plot the data using pandas builtin plot functionality
#this time using a subplot for each variable in the DataFrame
ssm_data_masked.plot(subplots=True)
plt.show()
_images/plot_ascat_data_4_0.png

In[6]:

#plot raw and masked SSM data in one plot to compare them

ssm_data_raw.data['SSM'].plot(label='raw SSM data')
ssm_data_masked.data['SSM'].plot(label='masked SSM data')
plt.legend()
plt.show()
_images/plot_ascat_data_5_0.png

In[7]:

ascat_SWI_reader = ascat.Ascat_SWI(path_to_ascat_swi_data,path_to_grid_definition,
                                   advisory_flags_path = path_to_adv_flags)


#reads swi data nearest to this lon,lat coordinates
#without any additional keywords all unmasked T values and
#Quality flags will be read
swi_data_raw = ascat_SWI_reader.read_swi(lon,lat)

#plot the data using pandas builtin plot functionality
swi_data_raw.plot()
plt.show()
_images/plot_ascat_data_6_0.png

In[8]:

#read the same data but this time only SWI with a T value
#of 20 is returned
swi_data_T_20 = ascat_SWI_reader.read_swi(lon,lat,T=20)

#plot the data using pandas builtin plot functionality
#this time using a subplot for each variable in the DataFrame
swi_data_T_20.plot(subplots=True)
plt.show()
_images/plot_ascat_data_7_0.png

In[9]:

#you can also mask manually if you prefer
swi_data_T_20.data = swi_data_T_20.data[swi_data_T_20.data['frozen_prob'] < 10]

swi_data_T_20.plot(subplots=True)
plt.show()
_images/plot_ascat_data_8_0.png

Reading and plotting data from the ISMN

This example program chooses a random Network and Station and plots the first variable,depht,sensor combination. To see how to get data for a variable from all stations see the next example.

It can be found in the /bin folder of the pytesmo package under the name plot_ISMN_data.py.

In[1]:

import pytesmo.io.ismn.interface as ismn
import os
import matplotlib.pyplot as plt
import random

In[2]:

#path unzipped file downloaded from the ISMN web portal
#on windows the first string has to be your drive letter
#like 'C:\\'
path_to_ismn_data = os.path.join('D:\\','small_projects','cpa_2013_07_ISMN_userformat_reader',
                                      'header_values_parser_test')

In[3]:

#initialize interface, this can take up to a few minutes the first
#time, since all metadata has to be collected
ISMN_reader = ismn.ISMN_Interface(path_to_ismn_data)

#plot available station on a map
ISMN_reader.plot_station_locations()
_images/plot_ISMN_2_0.png

In[4]:

#select random network and station to plot
networks = ISMN_reader.list_networks()
print "Available Networks:"
print networks
Available Networks:
['OZNET']

In[5]:

network = random.choice(networks)
stations = ISMN_reader.list_stations(network = network)
print "Available Stations in Network %s"%network
print stations
Available Stations in Network OZNET
['Alabama' 'Balranald-Bolton_Park' 'Banandra' 'Benwerrin' 'Bundure'
 'Canberra_Airport' 'Cheverelis' 'Cooma_Airfield' 'Cootamundra_Aerodrome'
 'Cox' 'Crawford' 'Dry_Lake' 'Eulo' 'Evergreen' 'Ginninderra_K4'
 'Ginninderra_K5' 'Griffith_Aerodrome' 'Hay_AWS' 'Keenan' 'Kyeamba_Downs'
 'Kyeamba_Mouth' 'Kyeamba_Station' 'Rochedale' 'S_Coleambally' 'Samarra'
 'Silver_Springs' 'Spring_Bank' 'Strathvale' 'Uri_Park' 'Waitara'
 'Weeroona' 'West_Wyalong_Airfield' 'Widgiewa' 'Wollumbi' 'Wynella'
 'Yamma_Road' 'Yammacoona' 'Yanco_Research_Station']

In[6]:

station = random.choice(stations)
station_obj = ISMN_reader.get_station(station)
print "Available Variables at Station %s"%station
#get the variables that this station measures
variables = station_obj.get_variables()
print variables
Available Variables at Station Evergreen
['precipitation' 'soil moisture' 'soil temperature']

In[7]:

#to make sure the selected variable is not measured
#by different sensors at the same depths
#we also select the first depth and the first sensor
#even if there is only one
depths_from,depths_to = station_obj.get_depths(variables[0])

sensors = station_obj.get_sensors(variables[0],depths_from[0],depths_to[0])

#read the data of the variable, depth, sensor combination
time_series = station_obj.read_variable(variables[0],depth_from=depths_from[0],depth_to=depths_to[0],sensor=sensors[0])

#print information about the selected time series
print "Selected time series is:"
print time_series
Selected time series is:
OZNET Evergreen -0.50 m - -0.50 m precipitation measured with TB4-0.2-mm-tipping-bucket-raingauge

In[8]:

#plot the data
time_series.plot()
#with pandas 0.12 time_series.plot() also works
plt.legend()
plt.show()
_images/plot_ISMN_7_0.png

In[9]:

#we also want to see soil moisture
sm_depht_from,sm_depht_to = station_obj.get_depths('soil moisture')
print sm_depht_from,sm_depht_to
[ 0.   0.   0.3  0.6] [ 0.05  0.3   0.6   0.9 ]

In[10]:

#read sm data measured in first layer 0-0.05m
sm = station_obj.read_variable('soil moisture',depth_from=0,depth_to=0.05)
sm.plot()
plt.show()
_images/plot_ISMN_9_0.png

Comparing ASCAT and insitu data from the ISMN

This example program loops through all insitu stations that measure soil moisture with a depth between 0 and 0.1m it then finds the nearest ASCAT grid point and reads the ASCAT data. After temporal matching and scaling using linear CDF matching it computes several metrics, like the correlation coefficients(Pearson’s, Spearman’s and Kendall’s), Bias, RMSD as well as the Nash–Sutcliffe model efficiency coefficient.

It is stopped after 2 stations to not take to long to run and produce a lot of plots

It can be found in the /bin folder of the pytesmo package under the name compare_ISMN_ASCAT.py.

In[1]:

import pytesmo.io.ismn.interface as ismn
import pytesmo.io.sat.ascat as ascat
import pytesmo.temporal_matching as temp_match
import pytesmo.scaling as scaling
import pytesmo.metrics as metrics

import os
import matplotlib.pyplot as plt

In[2]:

#set paths for ASCAT SSM
path_to_ascat_ssm_data = os.path.join('D:\\','small_projects',
                                      'cpa_2013_07_userformat_reader','data','ASCAT_SSM_25km_ts_WARP5.5_R0.1','data')

path_to_grid_definition = os.path.join('D:\\','small_projects',
                                      'cpa_2013_07_userformat_reader','data','auxiliary_data','grid_info')

path_to_adv_flags = os.path.join('D:\\','small_projects',
                                      'cpa_2013_07_userformat_reader','data','auxiliary_data','advisory_flags')

In[3]:

#set path to ISMN data
path_to_ismn_data =os.path.join('D:\\','small_projects','cpa_2013_07_ISMN_userformat_reader','header_values_parser_test')

In[4]:

#Initialize readers
ascat_SSM_reader = ascat.Ascat_SSM(path_to_ascat_ssm_data,path_to_grid_definition,
                                   advisory_flags_path = path_to_adv_flags)

ISMN_reader = ismn.ISMN_Interface(path_to_ismn_data)

In[5]:

i = 0

label_ascat='SSM'
label_insitu='insitu_sm'

In[6]:

#this loops through all stations that measure soil moisture
for station in ISMN_reader.stations_that_measure('soil moisture'):

    #this loops through all time series of this station that measure soil moisture
    #between 0 and 0.1 meters
    for ISMN_time_series in station.data_for_variable('soil moisture',min_depth=0,max_depth=0.1):

        ascat_time_series = ascat_SSM_reader.read_ssm(ISMN_time_series.longitude,
                                                      ISMN_time_series.latitude,
                                                      mask_ssf=True,
                                                      mask_frozen_prob = 5,
                                                      mask_snow_prob = 5)


        #drop nan values before doing any matching
        ascat_time_series.data = ascat_time_series.data.dropna()
        ISMN_time_series.data = ISMN_time_series.data.dropna()

        #rename the soil moisture column in ISMN_time_series.data to insitu_sm
        #to clearly differentiate the time series when they are plotted together
        ISMN_time_series.data.rename(columns={'soil moisture':label_insitu},inplace=True)

        #get ISMN data that was observerd within +- 1 hour(1/24. day) of the ASCAT observation
        #do not include those indexes where no observation was found
        matched_ISMN_data = temp_match.df_match(ascat_time_series.data,ISMN_time_series.data,
                                                window=1/24.,dropna=True)
        #matched ISMN data is now a dataframe with the same datetime index
        #as ascat_time_series.data and the nearest insitu observation

        #temporal matching also includes distance information
        #but we are not interested in it right now so let's drop it
        matched_ISMN_data = matched_ISMN_data.drop(['distance'],axis=1)

        #this joins the SSM column of the ASCAT data to the matched ISMN data
        matched_data = matched_ISMN_data.join(ascat_time_series.data[label_ascat])

        #continue only with relevant columns
        matched_data = matched_data[[label_ascat,label_insitu]]

        #the plot shows that ISMN and ASCAT are observed in different units
        matched_data.plot(figsize=(15,5),secondary_y=[label_ascat])
        plt.show()

        #this takes the matched_data DataFrame and adds a column
        scaled_data = scaling.add_scaled(matched_data, method='lin_cdf_match',
                                         label_in=label_ascat,label_scale=label_insitu)
        #the label of the scaled data is construced as label_in+'_scaled_'+method
        scaled_ascat_label = label_ascat+'_scaled_'+'lin_cdf_match'
        #now the scaled ascat data and insitu_sm are in the same space
        scaled_data.plot(figsize=(15,5),secondary_y=[label_ascat])
        plt.show()

        plt.scatter(scaled_data[scaled_ascat_label].values,scaled_data[label_insitu].values)
        plt.xlabel(scaled_ascat_label)
        plt.ylabel(label_insitu)
        plt.show()

        #calculate correlation coefficients, RMSD, bias, Nash Sutcliffe
        x, y = scaled_data[scaled_ascat_label].values, scaled_data[label_insitu].values

        print "ISMN time series:",ISMN_time_series
        print "compared to"
        print ascat_time_series
        print "Results:"
        print "Pearson's (R,p_value)", metrics.pearsonr(x, y)
        print "Spearman's (rho,p_value)", metrics.spearmanr(x, y)
        print "Kendalls's (tau,p_value)", metrics.kendalltau(x, y)
        print "RMSD", metrics.rmsd(x, y)
        print "Bias", metrics.bias(x, y)
        print "Nash Sutcliffe", metrics.nash_sutcliffe(x, y)


    i += 1

    #only show the first 2 stations, otherwise this program would run a long time
    #and produce a lot of plots
    if i >= 2:
        break
_images/compare_ASCAT_ISMN_5_0.png _images/compare_ASCAT_ISMN_5_1.png _images/compare_ASCAT_ISMN_5_2.png
ISMN time series: OZNET Alabama 0.00 m - 0.05 m soil moisture measured with Stevens-Hydra-Probe
compared to
ASCAT time series gpi:1884359 lat:-35.342 lon:147.541
Results:
Pearson's (R,p_value) (0.59736953256517777, 1.4810058830429653e-60)
Spearman's (rho,p_value) (0.63684906343988457, 4.8971200217989799e-71)
Kendalls's (tau,p_value) (0.45994629380576146, 4.6771942474849024e-65)
RMSD 0.0807313501609
Bias 0.00258302466701
Nash Sutcliffe 0.221824420266
_images/compare_ASCAT_ISMN_5_4.png _images/compare_ASCAT_ISMN_5_5.png _images/compare_ASCAT_ISMN_5_6.png
ISMN time series: OZNET Balranald-Bolton_Park 0.00 m - 0.08 m soil moisture measured with CS615
compared to
ASCAT time series gpi:1821003 lat:-33.990 lon:146.381
Results:
Pearson's (R,p_value) (0.65811087356086551, 9.1620935528699124e-126)
Spearman's (rho,p_value) (0.65874491635978671, 4.3707663858540222e-126)
Kendalls's (tau,p_value) (0.48451720923430946, 4.6613967263363183e-117)
RMSD 0.0283269899964
Bias -0.000181669876467
Nash Sutcliffe 0.314284186192