Instructions for download and extraction of data

2123df7e · Lars Kaleschke · eb3ea04c · 2123df7e
Commit 2123df7e authored 2 years ago by Lars Kaleschke
--- a/SMOS_ESSD_Ant.ipynb
+++ b/SMOS_ESSD_Ant.ipynb
@@ -24,7 +24,21 @@
    "\n",
    " * Helicopter-based EM Bird (HEM) data¶\n",
    " * Surface and Under-Ice Trawl (SUIT)\n",
-    " * Upward-Looking Sonars (ULS)\n"
+    " * Upward-Looking Sonars (ULS)\n",
+    " \n",
+    "The data sources are given below.\n",
+    "\n",
+    "## Module versions\n",
+    "\n",
+    "The code runs with the following versions:\n",
+    "\n",
+    "    python                    3.9.9           h62f1059_0_cpython    conda-forge\n",
+    "    xarray                    0.19.0             pyhd8ed1ab_0    conda-forge\n",
+    "    matplotlib                3.4.2            py39hf3d152e_0    conda-forge\n",
+    "    basemap                   1.2.2            py39h2103e0b_3    conda-forge\n",
+    "    pandas                    1.3.1            py39hde0f152_0    conda-forge\n",
+    "    pyproj                    3.2.1            py39ha9a7ae0_2    conda-forge\n",
+    "    "
   ]
  },
  {
@@ -229,7 +243,7 @@
    "\n",
    "Behrendt, A., Dierking, W., Fahrbach, E., and Witte, H.: Sea ice draft in the Weddell Sea, measured by upward looking sonars, Earth Syst. Sci. Data, 5, 209–226, https://doi.org/10.5194/essd-5-209-2013, 2013.\n",
    "\n",
-    "Behrendt, A., W. Dierking, and H. Witte(2015), Thermodynamic sea ice growthin the central Weddell Sea, observed inupward-looking sonar data,J. Geophys.Res. Oceans,120, 2270–2286,doi:10.1002/2014JC010408\n",
+    "Behrendt, A., W. Dierking, and H. Witte (2015), Thermodynamic sea ice growth in the central Weddell Sea, observed in upward-looking sonar data, J. Geophys. Res. Oceans,120, 2270–2286,doi:10.1002/2014JC010408\n",
    "\n",
    "The sea ice thickness in the Atlantic sector of the Southern Ocean, PhD thesis Axel Behrendt, 2013\n",
    "http://wwwdb.dbod.de/login?url=https://d-nb.info/1072077566/34\n",
@@ -274,9 +288,9 @@
   ],
   "source": [
    "# Download data from Pangaea https://doi.org/10.5194/essd-5-209-2013 \n",
-    "# and store in the directory ULS_dir\n",
+    "# and store in the directory ULS_dir.\n",
    "# Download ZIP file containing all datasets as tab-delimited text\n",
-    "\n",
+    "# unzip Behrendt_2013.zip\n",
    "ULS_dir=\"repo_data/ULS_Ant_Behrendt_2013/\"\n",
    "!ls $ULS_dir"
   ]
 %% Cell type:markdown id:32432818 tags:
 # SMOS-derived Antarctic thin sea-ice thickness from 2010 to 2020: data description and validation in the Weddell Sea
 This is the main notebook that generates figures and statistics used in the ESSD paper with the same title.
 ## SMOS data
 Tian-Kunze, Xiangshan; Kaleschke, Lars (2021): SMOS-derived sea ice thickness in the Antarctic from 2010 to 2020. PANGAEA, https://doi.org/10.1594/PANGAEA.934732
 As a first step the notebook entitled **SMOS sea ice thickness monthly mean and climatology** has to be used to generate monthly means and the climatology from the daily product.
 ## Validation data
 Three main datasets are used
 * Helicopter-based EM Bird (HEM) data¶
 * Surface and Under-Ice Trawl (SUIT)
 * Upward-Looking Sonars (ULS)
+The data sources are given below.
+## Module versions
+The code runs with the following versions:
+    python                    3.9.9           h62f1059_0_cpython    conda-forge
+    xarray                    0.19.0             pyhd8ed1ab_0    conda-forge
+    matplotlib                3.4.2            py39hf3d152e_0    conda-forge
+    basemap                   1.2.2            py39h2103e0b_3    conda-forge
+    pandas                    1.3.1            py39hde0f152_0    conda-forge
+    pyproj                    3.2.1            py39ha9a7ae0_2    conda-forge
 %% Cell type:code id:718fd223 tags:
 ``` python
 # Directory for daily product (data source)
 smos_datadir ='/isibhv/projects/siral/product/smos/SMOS_v620_L3C_Sea_v3/south/'
 # Directory for program output
 smos_monthly ='repo_data/SMOS_monthly/'
 ```
 %% Cell type:markdown id:7102d92a tags:
 ## Helicopter-based EM Bird (HEM) data
 HEM data have been submitted to Pangaea but the reference is not yet available. Therefore, the HEM data are included in this repository.
 %% Cell type:code id:5a5c52c4 tags:
 ``` python
 # This directory contains the HEM data, profiles and gridded.
 grdvaldir='repo_data/HEM'
 valdir='repo_data/HEM'
 #!ls $grdvaldir
 !ls $valdir
 ```
 %% Output
    Ant_HEM_2013-06-19_1343_sit.nc	HEM_Ant29_20130619T134315_20130619T134258.nc
    Ant_HEM_2013-06-20_1405_sit.nc	HEM_Ant29_20130620T140530_20130620T132327.nc
    Ant_HEM_2013-06-21_1145_sit.nc	HEM_Ant29_20130621T114536_20130621T114519.nc
    Ant_HEM_2013-07-07_1221_sit.nc	HEM_Ant29_20130707T122107_20130707T122051.nc
 %% Cell type:code id:c27b6af0 tags:
 ``` python
 #!gzip repo_data/HEM/*.nc
 !gunzip repo_data/HEM/*.gz
 ```
 %% Cell type:markdown id:50a679c0 tags:
 ## SMOS grid definitions
 %% Cell type:code id:712e8b75 tags:
 ``` python
 %matplotlib notebook
 from pylab import *
 from mpl_toolkits.basemap import Basemap
 from matplotlib.colors import LinearSegmentedColormap,ListedColormap,Normalize
 import cmocean
 import pandas as pd
 import xarray as xr
 import glob
 import os.path
 import pyproj
 pol=pyproj.Proj("+init=EPSG:3412") # NSIDC Polar Stereographic South
 def mapxy(x,y):
    lon,lat=pol(x*1000,y*1000, inverse=True)
    return lat,lon
 def mapll(lat,lon):
    x,y=pol(lon,lat, inverse=False)
    return x,y
 def def_grids():
    grids={}
    grids.update({'Arc': [1,3750.0,5850.0,-3850.0,-5350.0]})
    grids.update({'Ant': [-1,3950.0,4350.0,-3950.0,-3950.0]})
    return grids
 hemis='Ant'
 grids=def_grids()
 cs=12.5
 sgn,cds= grids[hemis][0],grids[hemis][1:5]
 X0,Y0,X2,Y2=cds[0],cds[1],cds[2],cds[3]
 x0,y0,x2,y2=cds[0],cds[1],cds[2],cds[3]
 ```
 %% Output
    /isibhv/projects/SMOS_ESL/isibhv_miniconda3/envs/jupyterlab/lib/python3.9/site-packages/pyproj/crs/crs.py:131: FutureWarning: '+init=<authority>:<code>' syntax is deprecated. '<authority>:<code>' is the preferred initialization method. When making the change, be mindful of axis order changes: https://pyproj4.github.io/pyproj/stable/gotchas.html#axis-order-changes-in-proj-6
      in_crs_string = _prepare_from_proj_string(in_crs_string)
 %% Cell type:markdown id:3f57b058 tags:
 ### SMOS climatological monthly mean and interannual standard deviation
 %% Cell type:code id:acca6c4f tags:
 ``` python
 smos_months=['04','05','06','07','08','09','10']
 sy,sx=(664, 632)
 SIT=zeros((sy,sx,7))
 fdir='/isibhv/projects/siral/product/smos/SMOS_v620_L3C_Sea_v3/south/monthly/'
 #fdir='/isibhv/projects/siral/product/smos/SMOS_v724_L3C_Sea_v3.3/south/monthly/'
 for i,month in enumerate(smos_months):
    smos_fn=fdir+'climatology_SMOS_Icethickness_v3.2_south_'+month+'.nc'
    #smos_fn=fdir+'monthly_SMOS_Icethickness_v3.2_south_2010'+month+'.nc'
    print(smos_fn)
    sit=xr.open_dataarray(smos_fn)
    SIT[:,:,i]=sit
 SIT=ma.masked_invalid(SIT)
 sitm=SIT.mean(axis=2)
 # Interannual variability (std)
 smos_months=['04','05','06','07','08','09','10']
 sy,sx=(664, 632)
 SITs=zeros((sy,sx,7))
 fdir='/isibhv/projects/siral/product/smos/SMOS_v620_L3C_Sea_v3/south/monthly/'
 #fdir='/isibhv/projects/siral/product/smos/SMOS_v724_L3C_Sea_v3.3/south/monthly/'
 for i,month in enumerate(smos_months):
    smos_fn=fdir+'climatology_SMOS_Icethickness_v3.2_south_'+month+'.nc'
    #smos_fn=fdir+'monthly_SMOS_Icethickness_v3.2_south_2010'+month+'.nc'
    flistpat=fdir+'monthly_SMOS_Icethickness_v3.2_south_20??'+month+'.nc'
    flist=!ls $flistpat
    nry=len(flist) # number years
    SITsy=zeros((sy,sx,nry))
    for j,smos_fn in enumerate(flist):
        #print(j,smos_fn)
        sit=xr.open_dataarray(smos_fn)
        SITsy[:,:,j]=sit
    SITsy=ma.masked_invalid(SITsy)
    sits=SITsy.std(axis=2)
    SITs[:,:,i]=sits
 ```
 %% Output
    /isibhv/projects/siral/product/smos/SMOS_v620_L3C_Sea_v3/south/monthly/climatology_SMOS_Icethickness_v3.2_south_04.nc
    /isibhv/projects/siral/product/smos/SMOS_v620_L3C_Sea_v3/south/monthly/climatology_SMOS_Icethickness_v3.2_south_05.nc
    /isibhv/projects/siral/product/smos/SMOS_v620_L3C_Sea_v3/south/monthly/climatology_SMOS_Icethickness_v3.2_south_06.nc
    /isibhv/projects/siral/product/smos/SMOS_v620_L3C_Sea_v3/south/monthly/climatology_SMOS_Icethickness_v3.2_south_07.nc
    /isibhv/projects/siral/product/smos/SMOS_v620_L3C_Sea_v3/south/monthly/climatology_SMOS_Icethickness_v3.2_south_08.nc
    /isibhv/projects/siral/product/smos/SMOS_v620_L3C_Sea_v3/south/monthly/climatology_SMOS_Icethickness_v3.2_south_09.nc
    /isibhv/projects/siral/product/smos/SMOS_v620_L3C_Sea_v3/south/monthly/climatology_SMOS_Icethickness_v3.2_south_10.nc
 %% Cell type:markdown id:37887b6b tags:
 ## ULS Data
 ### References
 Behrendt, A., Dierking, W., Fahrbach, E., and Witte, H.: Sea ice draft in the Weddell Sea, measured by upward looking sonars, Earth Syst. Sci. Data, 5, 209–226, https://doi.org/10.5194/essd-5-209-2013, 2013.
-Behrendt, A., W. Dierking, and H. Witte(2015), Thermodynamic sea ice growthin the central Weddell Sea, observed inupward-looking sonar data,J. Geophys.Res. Oceans,120, 2270–2286,doi:10.1002/2014JC010408
+Behrendt, A., W. Dierking, and H. Witte (2015), Thermodynamic sea ice growth in the central Weddell Sea, observed in upward-looking sonar data, J. Geophys. Res. Oceans,120, 2270–2286,doi:10.1002/2014JC010408
 The sea ice thickness in the Atlantic sector of the Southern Ocean, PhD thesis Axel Behrendt, 2013
 http://wwwdb.dbod.de/login?url=https://d-nb.info/1072077566/34
 Harms, S., Fahrbach, E., & Strass, V. H. (2001). Sea ice transports in the Weddell Sea. Journal of Geophysical Research: Oceans, 106(C5), 9057-9073.
 Behrendt, Axel; Dierking, Wolfgang; Fahrbach, Eberhard; Witte, Hannelore (2013): Sea ice draft measured by upward looking sonars in the Weddell Sea (Antarctica). PANGAEA, https://doi.org/10.1594/PANGAEA.785565, Supplement to: Behrendt, A et al. (2013): Sea ice draft in the Weddell Sea, measured by upward looking sonars. Earth System Science Data, 5(1), 209-226, https://doi.org/10.5194/essd-5-209-2013
 %% Cell type:code id:e6b5bd4c tags:
 ``` python
 # Download data from Pangaea https://doi.org/10.5194/essd-5-209-2013
-# and store in the directory ULS_dir
+# and store in the directory ULS_dir.
 # Download ZIP file containing all datasets as tab-delimited text
+# unzip Behrendt_2013.zip
 ULS_dir="repo_data/ULS_Ant_Behrendt_2013/"
 !ls $ULS_dir
 ```
 %% Output
    AWI206-4_sea-ice_draft.tab  AWI229-6_sea-ice_draft.tab
    AWI206-6_sea-ice_draft.tab  AWI229-8_sea-ice_draft.tab
    AWI207-2_sea-ice_draft.tab  AWI230-2_sea-ice_draft.tab
    AWI207-4_sea-ice_draft.tab  AWI230-3_sea-ice_draft.tab
    AWI207-6_sea-ice_draft.tab  AWI231-1_sea-ice_draft.tab
    AWI207-7_sea-ice_draft.tab  AWI231-2_sea-ice_draft.tab
    AWI208-3_sea-ice_draft.tab  AWI231-3_sea-ice_draft.tab
    AWI208-5_sea-ice_draft.tab  AWI231-4_sea-ice_draft.tab
    AWI209-3_sea-ice_draft.tab  AWI231-6_sea-ice_draft.tab
    AWI210-2_sea-ice_draft.tab  AWI231-7_sea-ice_draft.tab
    AWI212-2_sea-ice_draft.tab  AWI232-1_sea-ice_draft.tab
    AWI217_sea-ice_draft.tab    AWI232-4_sea-ice_draft.tab
    AWI227-3_sea-ice_draft.tab  AWI232-5_sea-ice_draft.tab
    AWI227-4_sea-ice_draft.tab  AWI232-6_sea-ice_draft.tab
    AWI227-6_sea-ice_draft.tab  AWI232-8_sea-ice_draft.tab
    AWI229-1_sea-ice_draft.tab  AWI232-9_sea-ice_draft.tab
    AWI229-2_sea-ice_draft.tab  AWI233-2_sea-ice_draft.tab
    AWI229-3_sea-ice_draft.tab  AWI233-6_sea-ice_draft.tab
    AWI229-4_sea-ice_draft.tab  summary.txt
    AWI229-5_sea-ice_draft.tab
 %% Cell type:markdown id:2a8102d4 tags:
 ### Functions to read in ULS data
 %% Cell type:code id:9334555e tags:
 ``` python
 from pylab import *
 import pandas as pd
 import xarray as xr
 import glob
 import os.path
 def ts_fillnan(ts):
    dr=pd.date_range(ts.index[0],ts.index[-1])
    ind=pd.Series(data=ones(len(dr)),index=dr)
    ind[:]=nan
    ind[ts.index]=ts[ts.index]
    return ind
 def headerlenght(fn):
    fid=open(fn)
    for l in range(40):
        line=fid.readline()
        if "*/" in line:
            return(l+1)
 def latlon_fromheader(fn):
    fid=open(fn)
    for l in range(40):
        line=fid.readline()
        if "Coverage:" in line:
            lat,lon=line.replace('\n','').split('LATITUDE: ')[1].split('* LONGITUDE:')
            lat,lon=float(lat),float(lon)
            return(lat,lon)
 fl=!ls $ULS_dir/*.tab
 L=[]
 for fn in fl:
    awi_id=os.path.basename(fn)[0:6]
    L.append(awi_id)
 awi_ids=unique(array(L))
 print('Unique ULS data positions', awi_ids)
 ```
 %% Output
    Unique ULS data positions ['AWI206' 'AWI207' 'AWI208' 'AWI209' 'AWI210' 'AWI212' 'AWI217' 'AWI227'
     'AWI229' 'AWI230' 'AWI231' 'AWI232' 'AWI233']
 %% Cell type:markdown id:05037629 tags:
 ### Read ULS data and calculate monthly means
 %% Cell type:code id:58aef07b tags:
 ``` python
 ULS_Dict_daily={}
 ULS_Dict_monthly={}
 ULS_latlon={}
 for uls_id in awi_ids:
    fl=glob.glob(ULS_dir+uls_id+'*.tab')
    fl.sort()
    ULS_dict={}
    for fn in fl:
        print(fn)
        ULS_name=os.path.basename(fn)[0:8]
        skiprows=headerlenght(fn)
        ULS_latlon[ULS_name[0:6]]=latlon_fromheader(fn)
        uls=pd.read_csv(fn,sep='\t',skiprows=skiprows)
        #print(ULS_name,uls.keys())
        for k in uls.keys():
            if 'Draft' in k:
                if 'model' in k:
                    draft_param=k
                else:
                    if 'zero line' in k:
                        draft_param=k
        print('--> Selected parameter:',draft_param)
        draft=array(uls[draft_param].values).astype(float)
        index=pd.to_datetime(uls['Date/Time'])
        index.name='Date'
        ULS_dict[ULS_name]=pd.Series(data=draft,index=index)
    ULS_1d={}
    for ULS_name in ULS_dict:
        print('Calculating daily averages', ULS_name)
        ULS_1d[ULS_name]=ULS_dict[ULS_name].resample('1d').mean()
    L=[]
    for uls in ULS_1d:
        L.append(ULS_1d[uls])
    ts=pd.concat(L)
    ts1m=ts.resample('m').mean()
    ULS_Dict_monthly[uls_id]=ts1m
    print('Finished ',uls_id)
 xa1=xr.Dataset(ULS_Dict_monthly)# Monthly mean
 ```
 %% Output
    repo_data/ULS_Ant_Behrendt_2013/AWI206-4_sea-ice_draft.tab
    --> Selected parameter: Draft [m] (zero line correction, Calculated)
    repo_data/ULS_Ant_Behrendt_2013/AWI206-6_sea-ice_draft.tab
    --> Selected parameter: Draft [m] (zero line correction, Calculated)
    Calculating daily averages AWI206-4
    Calculating daily averages AWI206-6
    Finished  AWI206
    repo_data/ULS_Ant_Behrendt_2013/AWI207-2_sea-ice_draft.tab
    --> Selected parameter: Draft [m] (model correction, Calculated)
    repo_data/ULS_Ant_Behrendt_2013/AWI207-4_sea-ice_draft.tab
    --> Selected parameter: Draft [m] (model correction, Calculated)
    repo_data/ULS_Ant_Behrendt_2013/AWI207-6_sea-ice_draft.tab
    --> Selected parameter: Draft [m] (model correction, Calculated)
    repo_data/ULS_Ant_Behrendt_2013/AWI207-7_sea-ice_draft.tab
    --> Selected parameter: Draft [m] (model correction, Calculated)
    Calculating daily averages AWI207-2
    Calculating daily averages AWI207-4
    Calculating daily averages AWI207-6
    Calculating daily averages AWI207-7
    Finished  AWI207
    repo_data/ULS_Ant_Behrendt_2013/AWI208-3_sea-ice_draft.tab
    --> Selected parameter: Draft [m] (model correction, Calculated)
    repo_data/ULS_Ant_Behrendt_2013/AWI208-5_sea-ice_draft.tab
    --> Selected parameter: Draft [m] (model correction, Calculated)
    Calculating daily averages AWI208-3
    Calculating daily averages AWI208-5
    Finished  AWI208
    repo_data/ULS_Ant_Behrendt_2013/AWI209-3_sea-ice_draft.tab
    --> Selected parameter: Draft [m] (model correction, Calculated)
    Calculating daily averages AWI209-3
    Finished  AWI209
    repo_data/ULS_Ant_Behrendt_2013/AWI210-2_sea-ice_draft.tab
    --> Selected parameter: Draft [m] (model correction, Calculated)
    Calculating daily averages AWI210-2
    Finished  AWI210
    repo_data/ULS_Ant_Behrendt_2013/AWI212-2_sea-ice_draft.tab
    --> Selected parameter: Draft [m] (model correction, Calculated)
    Calculating daily averages AWI212-2
    Finished  AWI212
    repo_data/ULS_Ant_Behrendt_2013/AWI217_sea-ice_draft.tab
    --> Selected parameter: Draft [m] (model correction, Calculated)
    Calculating daily averages AWI217_s
    Finished  AWI217
    repo_data/ULS_Ant_Behrendt_2013/AWI227-3_sea-ice_draft.tab
    --> Selected parameter: Draft [m] (zero line correction, Calculated)
    repo_data/ULS_Ant_Behrendt_2013/AWI227-4_sea-ice_draft.tab
    --> Selected parameter: Draft [m] (model correction, Calculated)
    repo_data/ULS_Ant_Behrendt_2013/AWI227-6_sea-ice_draft.tab
    --> Selected parameter: Draft [m] (model correction, Calculated)
    Calculating daily averages AWI227-3
    Calculating daily averages AWI227-4
    Calculating daily averages AWI227-6
    Finished  AWI227
    repo_data/ULS_Ant_Behrendt_2013/AWI229-1_sea-ice_draft.tab
    --> Selected parameter: Draft [m] (zero line correction, Calculated)
    repo_data/ULS_Ant_Behrendt_2013/AWI229-2_sea-ice_draft.tab
    --> Selected parameter: Draft [m] (zero line correction, Calculated)
    repo_data/ULS_Ant_Behrendt_2013/AWI229-3_sea-ice_draft.tab
    --> Selected parameter: Draft [m] (zero line correction, Calculated)
    repo_data/ULS_Ant_Behrendt_2013/AWI229-4_sea-ice_draft.tab
    --> Selected parameter: Draft [m] (zero line correction, Calculated)
    repo_data/ULS_Ant_Behrendt_2013/AWI229-5_sea-ice_draft.tab
    --> Selected parameter: Draft [m] (zero line correction, Calculated)
    repo_data/ULS_Ant_Behrendt_2013/AWI229-6_sea-ice_draft.tab
    --> Selected parameter: Draft [m] (zero line correction, Calculated)
    repo_data/ULS_Ant_Behrendt_2013/AWI229-8_sea-ice_draft.tab
    --> Selected parameter: Draft [m] (zero line correction, Calculated)
    Calculating daily averages AWI229-1
    Calculating daily averages AWI229-2
    Calculating daily averages AWI229-3
    Calculating daily averages AWI229-4
    Calculating daily averages AWI229-5
    Calculating daily averages AWI229-6
    Calculating daily averages AWI229-8
    Finished  AWI229
    repo_data/ULS_Ant_Behrendt_2013/AWI230-2_sea-ice_draft.tab
    --> Selected parameter: Draft [m] (zero line correction, Calculated)
    repo_data/ULS_Ant_Behrendt_2013/AWI230-3_sea-ice_draft.tab
    --> Selected parameter: Draft [m] (model correction, Calculated)
    Calculating daily averages AWI230-2
    Calculating daily averages AWI230-3
    Finished  AWI230
    repo_data/ULS_Ant_Behrendt_2013/AWI231-1_sea-ice_draft.tab
    --> Selected parameter: Draft [m] (model correction, Calculated)
    repo_data/ULS_Ant_Behrendt_2013/AWI231-2_sea-ice_draft.tab
    --> Selected parameter: Draft [m] (model correction, Calculated)
    repo_data/ULS_Ant_Behrendt_2013/AWI231-3_sea-ice_draft.tab
    --> Selected parameter: Draft [m] (model correction, Calculated)
    repo_data/ULS_Ant_Behrendt_2013/AWI231-4_sea-ice_draft.tab
    --> Selected parameter: Draft [m] (model correction, Calculated)
    repo_data/ULS_Ant_Behrendt_2013/AWI231-6_sea-ice_draft.tab
    --> Selected parameter: Draft [m] (model correction, Calculated)
    repo_data/ULS_Ant_Behrendt_2013/AWI231-7_sea-ice_draft.tab
    --> Selected parameter: Draft [m] (model correction, Calculated)
    Calculating daily averages AWI231-1
    Calculating daily averages AWI231-2
    Calculating daily averages AWI231-3
    Calculating daily averages AWI231-4
    Calculating daily averages AWI231-6
    Calculating daily averages AWI231-7
    Finished  AWI231
    repo_data/ULS_Ant_Behrendt_2013/AWI232-1_sea-ice_draft.tab
    --> Selected parameter: Draft [m] (model correction, Calculated)
    repo_data/ULS_Ant_Behrendt_2013/AWI232-4_sea-ice_draft.tab
    --> Selected parameter: Draft [m] (model correction, Calculated)
    repo_data/ULS_Ant_Behrendt_2013/AWI232-5_sea-ice_draft.tab
    --> Selected parameter: Draft [m] (model correction, Calculated)
    repo_data/ULS_Ant_Behrendt_2013/AWI232-6_sea-ice_draft.tab
    --> Selected parameter: Draft [m] (model correction, Calculated)
    repo_data/ULS_Ant_Behrendt_2013/AWI232-8_sea-ice_draft.tab
    --> Selected parameter: Draft [m] (model correction, Calculated)
    repo_data/ULS_Ant_Behrendt_2013/AWI232-9_sea-ice_draft.tab
    --> Selected parameter: Draft [m] (model correction, Calculated)
    Calculating daily averages AWI232-1
    Calculating daily averages AWI232-4
    Calculating daily averages AWI232-5
    Calculating daily averages AWI232-6
    Calculating daily averages AWI232-8
    Calculating daily averages AWI232-9
    Finished  AWI232
    repo_data/ULS_Ant_Behrendt_2013/AWI233-2_sea-ice_draft.tab
    --> Selected parameter: Draft [m] (model correction, Calculated)
    repo_data/ULS_Ant_Behrendt_2013/AWI233-6_sea-ice_draft.tab
    --> Selected parameter: Draft [m] (model correction, Calculated)
    Calculating daily averages AWI233-2
    Calculating daily averages AWI233-6
    Finished  AWI233
 %% Cell type:markdown id:95268a3d tags:
 ### Calculate seasonal cycle - climatological monthly means and interannual standard deviation
 %% Cell type:code id:1572d837 tags:
 ``` python
 ULS_seasonal={}
 ULS_seasonal_std={}
 for k in xa1:
    ts=xa1[k]
    dates=pd.to_datetime(ts.indexes['Date'])
    month_index=dates.to_period(freq='M')
    df=pd.DataFrame(ts)
    df=pd.DataFrame(data=ts.data,index=dates,columns=['Draft'])
    #df.index.name='Date'
    df['year-month']=month_index
    df['month']=month_index.month
    seasonal=df.groupby('month').mean()
    seasonal_std=df.groupby('month').std()
    ULS_seasonal[k]=seasonal
    ULS_seasonal_std[k]=seasonal_std
 xas1=xr.Dataset(ULS_seasonal)
 xass1=xr.Dataset(ULS_seasonal_std)
 xas2=xr.Dataset(ULS_seasonal,attrs=ULS_latlon)
 #xas2.to_netcdf('ULS_draft_monthly_seasonal.nc')
 ```
 %% Cell type:markdown id:2cd221dc tags:
 ### Select suitable ULS data
 The data have to fulfill following criteria:
 * Range of SMOS sea ice thickness 0 - 1.5 m
 * Long time series, at least 3 years of data
 %% Cell type:code id:0e0dc14c tags:
 ``` python
 for k in xa1:
    if xa1[k].max()<1.5:
        xa1[k].plot(label=k)
 ylabel('Draft [m]')
 legend(fontsize=8)
 title('ULS < 1.5 m')
 figure()
 #Selected: AWI227, AWI229, AWI231
 k='AWI227'
 xa1[k].plot(label=k)
 k='AWI229'
 xa1[k].plot(label=k)
 k='AWI231'
 xa1[k].plot(label=k)
 legend()
 title('Selected ULS with 47+ months of data')
 ```
 %% Output
    Text(0.5, 1.0, 'Selected ULS with 47+ months of data')
 %% Cell type:markdown id:98affb73 tags:
 ### Extract SMOS values at ULS positions
 %% Cell type:code id:08558000 tags:
 ``` python
 SMOS_dict={}
 SMOS_dict_s={}
 for k in xas2:
    lat,lon=xas2.attrs[k]
    ts1=xas2[k]
    psxy=mapll(lat,lon)
    posxy=(int((psxy[0]/1000-x2)/cs),int((psxy[1]/1000-y2)/cs))
    xi,yi=posxy
    #SIT_0[yi-1:yi+1,xi-1:xi+1]=-1
    #print(k,lat,lon,xi,yi)
    L=[nan]*3 # Month 1-3 not available
    Ls=[nan]*3
    for mi in range(7):
        L.append(SIT[yi,xi,mi])
        Ls.append(SITs[yi,xi,mi])
    L.append(nan)
    L.append(nan)
    Ls.append(nan)
    Ls.append(nan)
    SMOS_dict[k]=array(L)
    SMOS_dict_s[k]=array(Ls)
 ```
 %% Cell type:code id:32971b21 tags:
 ``` python
 selected=['AWI209','AWI227','AWI229','AWI231']
 selected=['AWI227','AWI229','AWI231']
 symdic={'AWI227':'o','AWI229':'v','AWI231':'^'}
 coldic={'AWI227':'b','AWI229':'r','AWI231':'c'}
 coldic={'AWI227':'k','AWI229':'k','AWI231':'k'}
 XL=[] # List for mean
 YL=[]
 XLS=[] # Std
 YLS=[]
 Lab=[] # Labels
 xx=arange(12)+1
 ci=1
 for k in selected:
    subplot(3,1,ci)
    lat,lon=xas2.attrs[k]
    mu=float(xas2[k].mean())
    mumax=float(xas2[k].max())
    nrm=sum(isfinite(array(xa1[k])))
    print(k,nrm)
    ts1=xas2[k].values[:,0]
    ts1=0.028+ts1*1.012
    tss1=xass1[k].values[:,0]
    fill_between(xx,ts1-tss1,ts1+tss1,alpha=0.3,color=coldic[k])
    plot(xx,ts1,label='ULS '+k,linestyle=':',color=coldic[k])
    ts2=SMOS_dict[k]
    tss2=SMOS_dict_s[k]
    errorbar(xx,ts2,yerr=tss2,color=coldic[k],label='SMOS @'+k)
    ci+=1
    xlim([3.9,10.5])
    ylim([0,1.3])
    legend(fontsize=8,loc=2)
    xticks([])
    ylabel('Thickness [m]')
    # Store value pairs into list (mean, std)
    ind=isfinite(ts2)
    YL.extend(ts2[ind])
    XL.extend(ts1[ind])
    XLS.extend(tss1[ind])
    YLS.extend(tss2[ind])
    Lab.extend(array(((k+',')*12).split(',')[:-1])[ind])
 xticks(ticks=[4,5,6,7,8,9,10],labels=['April','May','June','July','Aug.','Sept.','Oct.'])
 tight_layout()
 #savefig('Seasonal_thickness_ULS_SMOS.pdf')
 x=array(XL) # ULS
 y=array(YL) # SMOS
 xstd=array(XLS)
 ystd=array(YLS)
 ```
 %% Output
    AWI227 47
    AWI229 129
    AWI231 99
 %% Cell type:code id:0041713f tags:
 ``` python
 from sklearn.metrics import mean_squared_error
 import scipy.stats as stats
 slope, intercept, r_value, p_value, std_err = stats.linregress(x, y)
 print('Standard metric -----------')
 print("Mean VAL, SMOS: %1.2f +- %1.2f ,  %1.2f +- %1.2f"  % (mean(x),std(x),mean(y),std(y)))
 print("N: ",len(x))
 md=np.mean(x-y)
 print("Mean difference VAL-SMOS: %1.2f" % md)
 rmsd=mean_squared_error(x,y,squared=False)
 print("RMSD VAL-SMOS: %1.2f" % rmsd)
 print("R: %1.2f" % r_value)
 plt.figure(figsize=(6,6))
 xn=np.linspace(0,1.4)
 plt.plot(xn,xn,'k:',label='Identity')
 plt.plot(xn,xn*slope+intercept,'k-',label='Regression')
 for i,lab in enumerate(Lab):
    marker=symdic[lab]
    if i/7==int(i/7):
        label=lab
    else:
        label=''
    plt.plot([x[i]],[y[i]],color='k',marker=marker,ms=10,alpha=1.0,label=label)
 plt.errorbar(x,y,xerr=xstd,yerr=ystd,fmt='.',color='k')
 grid()
 xlabel('ULS sea ice thickness [m]')
 ylabel('SMOS sea ice thickness [m]')
 xlim([-0.1,1.3])
 ylim([-0.1,1.3])
 legend()
 tight_layout()
 #savefig('Scatter_SMOS_ULS_seasonal.pdf')
 #savefig('Scatter_SMOS_ULS_seasonal.png')
 ```
 %% Output
    Standard metric -----------
    Mean VAL, SMOS: 0.46 +- 0.32 ,  0.49 +- 0.38
    N:  21
    Mean difference VAL-SMOS: -0.03
    RMSD VAL-SMOS: 0.14
    R: 0.94
 %% Cell type:markdown id:68f8da1c tags:
 ### SUIT
 Castellani, Giulia; Flores, Hauke; Lange, Benjamin Allen; Schaafsma, Fokje L; Ehrlich, Julia; David, Carmen L; van Franeker, Jan Andries; Meijboom, Andre; van Dorssen, Michiel; Vortkamp, Martina; Van de Putte, Anton P (2019): Sea ice draft, under-ice irradiance and radiance, and surface water temperature, salinity and chlorophyl a from Surface and Under Ice Trawl (SUIT) measurements in 2012-2017. Alfred Wegener Institute, Helmholtz Centre for Polar and Marine Research, Bremerhaven, PANGAEA, https://doi.org/10.1594/PANGAEA.902056
 %% Cell type:code id:9506be52 tags:
 ``` python
 fn0='repo_data/PANGAEA_902260_PS81.tab'
 dat=pd.read_csv(fn0,sep='\t',skiprows=35)
 di=pd.to_datetime(dat['Date/Time'])+pd.to_timedelta(dat['Time [s]'].values,unit='s')
 Lat=pd.Series(data=dat['Latitude'].values,index=di)
 Lon=pd.Series(data=dat['Longitude'].values,index=di)
 Thick=pd.Series(data=dat['Ice thick [m]'].values,index=di)
 Distance=pd.Series(data=dat['Distance [m]'].values,index=di)
 from numpy import unique, array
 days=unique(array(dat['Date/Time']))
 for day in days:
    print(day)
 dat.head()
 ```
 %% Output
    2013-08-29
    2013-09-09
    2013-09-10
    2013-09-11
    2013-09-12
    2013-09-16
    2013-09-28
    2013-09-29
    2013-09-30
    2013-10-02
            Event   Date/Time   Latitude  Longitude  Haul  Station      Time [s]  \
    0  PS81/549-1  2013-08-29 -61.251660 -42.066770     2      549  43027.002375
    1  PS81/549-1  2013-08-29 -61.249863 -42.066779     2      549  43028.032945
    2  PS81/549-1  2013-08-29 -61.252591 -42.066787     2      549  43029.054942
    3  PS81/549-1  2013-08-29 -61.251435 -42.066796     2      549  43030.066205
    4  PS81/549-1  2013-08-29 -61.251739 -42.066805     2      549  43031.082798
       Distance [m]  Sea ice draft [m]  Ice thick [m]  Chl a [mg/m**3]  Temp [°C]  \
    0           0.0          -0.802515       0.875152         0.122137  -1.859637
    1           0.5          -0.698829       0.762081         0.126799  -1.859591
    2           1.0          -0.426149       0.464721         0.130338  -1.859484
    3           1.5          -0.711586       0.775993         0.132577  -1.859596
    4           2.0          -0.801220       0.873740         0.134420  -1.859658
             Sal
    0  34.241937
    1  34.186280
    2  34.247329
    3  34.253710
    4  34.194110
 %% Cell type:code id:ea04ace6 tags:
 ``` python
 # Create SUIT profile plots, not used in the paper
 close('all')
 figure(figsize=(9,18))
 spn=1
 thresh=0.03
 L=[]
 for day in days:
    thick0=Thick.loc[day]
    dist0=Distance.loc[day]
    lat0=Lat.loc[day]
    lon0=Lon.loc[day]
    mlat=mean(lat0)
    mlon=mean(lon0)
    if len(thick0)>1000:
        subplot(9,1,spn)
        spn+=1
        plot(dist0,thick0,'.',ms=1,alpha=0.5)
        h_median=median(thick0)
        h_mean=mean(thick0)
        h0_str='$h_{median}$'+str('={:01.2f} m'.format(h_median))
        h1_str='$h_{mean}$'+str('={:01.2f} m'.format(h_mean))
        x0,x1=0,max(dist0)
        plot([x0,x1],[h_median,h_median],'k:',label=h0_str)
        plot([x0,x1],[h_mean,h_mean],'k.-',ms=0,lw=0.5,label=h1_str)
        # Ice concentration / open water fraction
        c0=thick0<=0.05
        c1=thick0>0.05
        ice_conc=100-sum(c0)/sum(c1)*100
        c_str=str('$C$={:01.1f} %'.format(ice_conc))
        plot(dist0[c0],thick0[c0],'r.',ms=1,alpha=0.5,label=c_str)
        legend(loc=1,fontsize=8,title=day)
        ylim([-0.1,5])
        L.append([day,mlat,mlon,h_median,h_mean,ice_conc])
        ylabel('Thickness [m]',fontsize=8)
 xlabel('SUIT profile length [m]')
 #savefig('figs/SUIT_profiles.pdf',dpi=200)
 L=array(L)
 D={'Date':L[:,0],'Latitude':L[:,1].astype(float),'Longitude':L[:,2].astype(float),
   'H median':L[:,3].astype(float),'H mean':L[:,4].astype(float),'Ice concentration':L[:,5].astype(float)}
 SUIT_xa=xr.Dataset(D)
 #SUIT_xa.to_netcdf('data/SUIT_thickness_validation.nc')
 #close('all')
 ```
 %% Output
 %% Cell type:markdown id:65151b57 tags:
 ### Generate overview map
 %% Cell type:code id:61064cd1 tags:
 ``` python
 import warnings
 warnings.filterwarnings('ignore')
 i=4
 # 0 April 1 May 2 June 3 July 4 August 5 Septmber 6 October
 #sit0=SIT[:,:,i]
 sit0=sitm
 #sit0=SIT[:,:] # September 2010
 if True:
    sit0=ma.masked_invalid(sit0)
    #Extract subimage
    x0,y0,x2,y2=1200,3700.0,-3000.0,700.0
    sit_region=sit0[int((y2-Y2)/cs):int((y0-Y2)/cs),int((x2-X2)/cs):int((x0-X2)/cs)]
    Cur=array(mapxy(x0,y0))
    Cll=array(mapxy(x2,y2))
    palette = cm.gist_earth#bone
    palette = cm.viridis
    palette = cm.plasma
    palette.set_over("#ffffff")
    palette.set_bad("#d0d0d0")
    palette.set_under("#0000f0")
    close()
    fig=plt.figure()
    m = Basemap(projection='stere',resolution='l',llcrnrlon=Cll[1],llcrnrlat=Cll[0],\
                urcrnrlon=Cur[1],urcrnrlat=Cur[0],lat_ts=-70,lon_0=0,lat_0=-90,rsphere=(6378273.,6356889.))#,round=True)
    m.imshow(sit_region,interpolation='nearest',cmap=palette,vmin=0,vmax=1.0)
    m.drawmeridians(arange(0, 360, 30),labels=[0,0,1,0],linewidth=0.5)
    m.drawparallels(arange(-90, 90, 5),labels=[0,0,0,1],linewidth=0.5)
    plt.colorbar(shrink=0.5,extend='max',ticks=[0,0.5,1.0],pad=0.01,orientation='vertical').set_label(label='Sea-ice thickness [m]',size=10)#.set_label_position('left')
    for k in ['AWI227', 'AWI229', 'AWI231']:
        lat,lon=xas2.attrs[k]
        #mu=float(xas2[k][i+4].mean())
        mu=float(xas2[k].mean())
        xsc,ysc=m(lon,lat)
        ec='k'
        mumax=float(xas2[k].max())
        label=''
        ec='w'
        label='ULS '+k
        msize=50
        scatter(xsc,ysc,s=msize,color=palette(mu),vmin=0,vmax=1,edgecolors=ec,label=label)
        va='bottom'
        text(xsc,ysc,k[3:],va=va,fontsize=10)
    # SUIT
    for day in days:
        thick0=mean(Thick.loc[day])
        dist0=Distance.loc[day]
        lat0=Lat.loc[day]
        lon0=Lon.loc[day]
        mlat=mean(lat0)
        mlon=mean(lon0)
        xsc,ysc=m(mlon,mlat)
        ec='g'
        label=''
        scatter(xsc,ysc,s=30,marker='s',color=palette(thick0),vmin=0,vmax=1,edgecolors=ec,label=label)
    label='SUIT'
    scatter(xsc,ysc,s=30,marker='s',color=palette(thick0),vmin=0,vmax=1,edgecolors=ec,label=label)
    # HEM
    fl=!ls $valdir/HEM_Ant29*.nc
    for fn in fl:
        xhem=xr.open_mfdataset(fn)
        thickness=xhem.THICKNESS.values
        ind=isfinite(thickness)
        thickness=thickness[ind]
        xlat=xhem.LATITUDE.values[ind][::100]
        xlon=xhem.LONGITUDE.values[ind][::100]
        x,y=m(xlon,xlat)
        #plot(x,y,'c.',ms=0.1)
        mxlat=mean(xlat)
        mxlon=mean(xlon)
        mthickness=mean(thickness)
        xsc,ysc=m(mxlon,mxlat)
        ec='c'
        label=''
        scatter(xsc,ysc,marker='d',s=40,color=palette(mthickness),vmin=0,vmax=1,edgecolors=ec,label=label)
    label='HEM'
    scatter(xsc,ysc,marker='d',s=40,color=palette(mthickness),vmin=0,vmax=1,edgecolors=ec,label=label)
 legend(loc='lower right',title='Validation sites',fontsize=6)
 tight_layout()
 #savefig('ULS_map_SMOS_mean.pdf')
 print('finished')
 ```
 %% Output
    finished
 %% Cell type:markdown id:33fa8270 tags:
 ## SUIT single days
 %% Cell type:code id:4972e4b9 tags:
 ``` python
 Outdir0=smos_datadir
 year2=2020 # This year
 year1=2010 # Start of SMOS data
 Ny=year2-year1
 years=list(range(year1,year2+1))
 dr=pd.date_range(str(year1)+'-01-01',str(year2)+'-12-31',freq='1d')
 #dr=dr.to_period('D')
 Nd=len(dr) # Total days in period
 nans=zeros(Nd)
 nans[:]=nan
 ts2=pd.Series(data=nans,index=dr)
 ts=pd.Series()
 df=pd.DataFrame()
 date_list=list(dr.astype(str))
 fl=[]
 L_date=[]
 for datestr in date_list:
    yyyy,mm,dd,hh=datestr[0:4],datestr[5:7],datestr[8:10],datestr[11:13]
    Outdir=Outdir0+'/'+yyyy+'/'+mm+'/'
    fn=Outdir+'SMOS_Icethickness_v3.2_south_'+yyyy+mm+dd+'.nc'
    fl1=glob.glob(Outdir+'*'+yyyy+mm+dd+'*.nc')
    if fl1==[]:
        fn=''
    else:
        fn=fl1[0]
    if os.path.exists(fn):
        L_date.append(pd.to_datetime( yyyy+mm+dd))
        fl.append(fn)
 ind=array(L_date)
 ind2=pd.Series(data=ones(len(ind)),index=ind)
 #ind2=ind2.to_period('12h')
 smos_index=ts2.add(ind2,fill_value=0)
 file_ind=pd.Series(array(fl),index=ind)
 figure(figsize=(15,5))
 isfinite(smos_index).astype(float).plot.area(stacked=False)
 #savefig('SMOS_L3_TB_NH_available.png',dpi=150)
 #savefig('SMOS_L3_SH_available.png',dpi=150)
 ```
 %% Output
    <AxesSubplot:>
 %% Cell type:code id:960ee8b3 tags:
 ``` python
 file_ind=pd.Series(array(fl),index=ind)
 fn0=file_ind[day]
 xa=xr.open_mfdataset(fn0)
 land=array(xa.land)[0,:,:]
 sit=array(xa.sea_ice_thickness)[0,:,:]
 sy,sx=sit.shape[0],sit.shape[1]
 lat=array(xa.latitude)
 lon=array(xa.longitude)
 tc='k'
 #sit=sit0
 ```
 %% Cell type:code id:d8f63bf7 tags:
 ``` python
 days=array(['2013-08-29', '2013-09-09', '2013-09-10', '2013-09-11',
       '2013-09-12', '2013-09-16', '2013-09-28', '2013-09-29',
       '2013-09-30', '2013-10-02'], dtype=object)
 days2=array(['2013-08-29', '2013-09-09', '2013-09-10', '2013-09-11',
       '2013-09-12', '2013-09-16', '2013-09-29',
       '2013-09-30', ], dtype=object)
 #days=days2
 dayL=[]
 fl0=[]
 for day in days:
    fn0=file_ind[day]
    fl0.append(fn0)
 xa=xr.open_mfdataset(fl0)
 sit=xa.sea_ice_thickness.mean(dim='time')
 situ=xa.ice_thickness_uncertainty.mean(dim='time')
 satr=xa.saturation_ratio.mean(dim='time')
 ```
 %% Cell type:code id:b05dfedd tags:
 ``` python
 SITL=[]
 vals=[]
 #days_selected=days
 days=array(['2013-08-29', '2013-09-09', '2013-09-10', '2013-09-11',
       '2013-09-12', '2013-09-16', '2013-09-28', '2013-09-29',
       '2013-09-30', '2013-10-02'], dtype=object)
 days2=array(['2013-08-29', '2013-09-09', '2013-09-10', '2013-09-11',
       '2013-09-12', '2013-09-16', '2013-09-29',
       '2013-09-30', ], dtype=object)
 #days=days2
 dayL=[]
 for day in days:
    fn0=file_ind[day]
    xa=xr.open_mfdataset(fn0)
    sit=ma.masked_invalid((xa.sea_ice_thickness)[0,:,:])
    situ=ma.masked_invalid((xa.ice_thickness_uncertainty)[0,:,:])
    satr=ma.masked_invalid((xa.saturation_ratio)[0,:,:])
    #SITL.append(sit)
    lat0=Lat.loc[day].median()
    lon0=Lon.loc[day].median()
    thick0=Thick.loc[day]#.median()
    mthick=thick0.median()
    #mthick=thick0.mean()
    sthick=thick0.std()
    q25thick=quantile(thick0,0.25)
    q75thick=quantile(thick0,0.75)
    c0=thick0<=0.05
    c1=thick0>0.05
    ice_conc=100-sum(c0)/sum(c1)*100
    #if day=='2013-10-02':
    if len(Thick.loc[day])>1000:
        iy,ix=unravel_index(argmin(abs(lat[:]-lat0)+abs(lon[:]-lon0)),lat.shape) # get closest index
        vals.append([mthick,sthick,q25thick,q75thick,ice_conc,(sit[iy,ix]),situ[iy,ix],satr[iy,ix]])
        dayL.append(day)
 v=array(vals)
 ```
 %% Cell type:code id:9a11a24b tags:
 ``` python
 from sklearn.metrics import mean_squared_error
 #x,y,xe,ye=v[:,0],v[:,3],v[:,1],v[:,4]
 x,y,xe0,xe1,ye1=v[:,0],v[:,5],v[:,2],v[:,3],v[:,6]
 dmax=v[:,5]/v[:,7]*100
 ye=array([ye1,dmax-y])
 xe=array([xe0,xe1])
 slope, intercept, r_value, p_value, std_err = stats.linregress(x, y)
 print("slope: %f    intercept: %f" % (slope, intercept))
 print("R: %f" % r_value)
 print("R-squared: %f" % r_value**2)
 print("Std err: %f" % std_err)
 md=np.mean(x-y)
 print("Mean difference VAL-SMOS: %f" % md)
 print('Standard metric -----------')
 print("Mean VAL, SMOS: %1.2f +- %1.2f ,  %1.2f +- %1.2f"  % (mean(x),std(x),mean(y),std(y)))
 print("N: ",len(x))
 md=np.mean(x-y)
 print("Mean difference VAL-SMOS: %1.2f" % md)
 rmsd=mean_squared_error(x,y,squared=False)
 print("RMSD VAL-SMOS: %1.2f" % rmsd)
 print("R: %1.2f" % r_value)
 #plt.plot(x,y,'o')
 plt.figure(figsize=(6,6))
 xn=np.linspace(0,2)
 plt.plot(xn,xn,'k:',label='Identity')
 #plt.plot(xn,xn*slope+intercept,'k-',label='Ordinary regression')
 plt.plot(x,y,'ko',ms=5,alpha=1.0,label='SMOS pixel based')
 plt.errorbar(x,y,ms=5,xerr=xe,yerr=ye,fmt='ko',alpha=1.0)
 #plt.errorbar(x,y,xerr=xe,yerr=ye,fmt='ko',alpha=0.2,label='SMOS pixel based')
 plt.plot(x,dmax,'r^',ms=5,alpha=1.0,label='$d_{max}$')
 #plt.plot(xn,xn*wslope+wintercept,'k:',label='Weighted regression')
 #plt.plot(xn,xn,'k:')
 #plt.grid()
 plt.xlabel('Validation thickness [m]')
 plt.ylabel('SMOS thickness [m]')
 #plt.title('SMOS SIT version 3.2')
 #plt.title('SMOS SIT version 2')
 plt.legend(loc=2)
 plt.xlim([-0.1,1.5])
 plt.ylim([-0.1,1.5])
 #plt.savefig('SMOS_v2_scatter_fits.png',dpi=200)
 for i,day in enumerate(dayL):
    text(x[i]-0.05,y[i]-0.05,day[5:],color='r',fontsize=8)
 plt.tight_layout()
 grid()
 #plt.savefig('figs/SMOS_SUIT_scatter.pdf')
 ```
 %% Output
    slope: 0.486081    intercept: 0.055360
    R: 0.714217
    R-squared: 0.510105
    Std err: 0.180045
    Mean difference VAL-SMOS: 0.246493
    Standard metric -----------
    Mean VAL, SMOS: 0.59 +- 0.24 ,  0.34 +- 0.16
    N:  9
    Mean difference VAL-SMOS: 0.25
    RMSD VAL-SMOS: 0.30
    R: 0.71
 %% Cell type:markdown id:fcda86b3 tags:
 ## HEM single days
 %% Cell type:code id:a4b3bcca tags:
 ``` python
 XA=xr.open_mfdataset('repo_data/HEM/Ant_*_*_sit.nc')
 XA
 ```
 %% Output
    <xarray.Dataset>
    Dimensions:            (x: 632, y: 664, time: 4)
    Coordinates:
      * x                  (x) float64 -3.95e+03 -3.937e+03 ... 3.937e+03 3.95e+03
      * y                  (y) float64 -3.95e+03 -3.937e+03 ... 4.337e+03 4.35e+03
      * time               (time) datetime64[ns] 2013-06-19T13:43:00 ... 2013-07-...
    Data variables:
        sea_ice_thickness  (time, y, x) float32 dask.array<chunksize=(1, 664, 632), meta=np.ndarray>
 %% Cell type:code id:ba8503a1 tags:
 ``` python
 zval=[]
 zsmos=[]
 zsmosu=[] # Uncertainty
 zsmosr=[] # Saturation
 zvalm=[] # Mean for days
 zsmosm=[]
 zvals=[] # std for days
 zsmoss=[]
 for ti in XA.time:
    SITval=ma.masked_invalid(XA.sel(time=ti).to_array())[0,:,:]
    datestr=str(np.array(ti))[0:10]
    yyyy,mm,dd,hh=datestr[0:4],datestr[5:7],datestr[8:10],datestr[11:13]
    Outdir=Outdir0+'/'+yyyy+'/'+mm+'/'
    fn=Outdir+'SMOS_Icethickness_v3.2_south_'+yyyy+mm+dd+'.nc'
    print(fn)
    xa=xr.open_mfdataset(fn)
    SIT0=np.array(xa.sea_ice_thickness)[0,:,:]
    SITU=np.array(xa.ice_thickness_uncertainty)[0,:,:]
    SATR=np.array(xa.saturation_ratio)[0,:,:]
    ind=np.logical_not(SITval.mask)
    zval.extend(SITval[ind].flatten().data)
    zsmos.extend(SIT0[ind].flatten().data)
    zsmosu.extend(SITU[ind].flatten().data)
    zsmosr.extend(SATR[ind].flatten().data)
    zvalm.append(mean(SITval[ind].flatten().data))
    zsmosm.append(mean(SIT0[ind].flatten().data))
    zvals.append(std(SITval[ind].flatten().data))
    zsmoss.append(std(SIT0[ind].flatten().data))
    print(len(zval))
 x=np.array(zval)
 y=np.array(zsmos)
 ye=np.array(zsmosu)
 dmax=y/np.array(zsmosr)*100
 xm=np.array(zvalm)
 ym=np.array(zsmosm)
 xms=np.array(zvals)
 yms=np.array(zsmoss)
 ```
 %% Output
    /isibhv/projects/siral/product/smos/SMOS_v620_L3C_Sea_v3/south//2013/06/SMOS_Icethickness_v3.2_south_20130619.nc
    11
    /isibhv/projects/siral/product/smos/SMOS_v620_L3C_Sea_v3/south//2013/06/SMOS_Icethickness_v3.2_south_20130620.nc
    25
    /isibhv/projects/siral/product/smos/SMOS_v620_L3C_Sea_v3/south//2013/06/SMOS_Icethickness_v3.2_south_20130621.nc
    41
    /isibhv/projects/siral/product/smos/SMOS_v620_L3C_Sea_v3/south//2013/07/SMOS_Icethickness_v3.2_south_20130707.nc
    53
 %% Cell type:code id:c42f298b tags:
 ``` python
 ye2=array([ye,dmax-y])
 ye3=array([0,dmax-y])
 plt.figure(figsize=(6,6))
 xn=np.linspace(0,2)
 plt.plot(xn,xn,'k:',label='Identity')
 plt.plot(x,y,'ko',ms=3,alpha=0.5,label='SMOS pixel based')
 plt.errorbar(xm,ym,xerr=xms,yerr=yms,ms=10,fmt='ko',label='Flight mean')
 plt.plot(x,dmax,'r^',ms=5,alpha=1.0,label='$d_{max}$')
 plt.errorbar(x,y,yerr=ye2,ms=3,fmt='ko',alpha=0.3)#,label='SMOS pixel based')
 plt.xlabel('Validation thickness [m]')
 plt.ylabel('SMOS thickness [m]')
 plt.legend(loc=2)
 plt.xlim([0.0,2.0])
 plt.ylim([0.0,2.0])
 plt.grid()
 #savefig('Ant_Val_HEM_Suit_Scatter.png',dpi=300)
 #savefig('figs/Ant_Val_HEM_Scatter.pdf')
 ```
 %% Output