From 740a8ffad6ba034ec79a54b09a86841f0d940171 Mon Sep 17 00:00:00 2001 From: rzinke Date: Wed, 11 Sep 2024 01:43:28 +0000 Subject: [PATCH] GNSS station checks --- .../Secular_Requirement_Validation.ipynb | 130 +++++++++++++----- 1 file changed, 94 insertions(+), 36 deletions(-) diff --git a/methods/secular/Secular_Requirement_Validation.ipynb b/methods/secular/Secular_Requirement_Validation.ipynb index 24e059c..1bf1fa1 100644 --- a/methods/secular/Secular_Requirement_Validation.ipynb +++ b/methods/secular/Secular_Requirement_Validation.ipynb @@ -496,6 +496,15 @@ "With this formulation, we can obtain InSAR and GNSS velocity estimates and their formal uncertainties (including in areas where zero deformation is expected). The default InSAR (and GNSS) velocity fit in MintPy is to estimate a mean linear velocity $(v)$ in in the equation. To estimate a secular deformation rate, one can also estimate additional parameters for annual and semi-annual terms to account for seasonal periodic fluctuations in position. Below, we define a set of basis functions that MintPy will invert for, for both InSAR and GNSS data." ] }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "
\n", + "If the time-series span less than one year, then annual and semi-annual fuctions will not be reliably fit for noisy data. Therefore, the period functions will not be included in these analyses if that is the case.\n", + "
" + ] + }, { "cell_type": "code", "execution_count": null, @@ -511,7 +520,18 @@ " 'polyline': [],\n", " 'exp': {},\n", " 'log': {}\n", - "}" + "}\n", + "\n", + "# Determine time span to check whether fitting of periodic functions is appropriate\n", + "ts_metadata = readfile.read_attribute(timeseries_filename)\n", + "start_date = dt.strptime(ts_metadata['START_DATE'], '%Y%m%d')\n", + "end_date = dt.strptime(ts_metadata['END_DATE'], '%Y%m%d')\n", + "ts_timespan = (end_date - start_date).days\n", + "\n", + "if ts_timespan < 365:\n", + " print('Time span is {:d} days'.format(ts_timespan))\n", + " ts_functions['periodic'] = []\n", + " print('No periodic functions will be fit')\n" ] }, { @@ -664,6 +684,13 @@ "### 4.4. Get GNSS Position Time Series " ] }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We download the GNSS station data, and evaluate the station data quality based on data completeness and scatter of residuals." + ] + }, { "cell_type": "code", "execution_count": null, @@ -672,56 +699,89 @@ }, "outputs": [], "source": [ - "use_stn = [] #stations to keep\n", - "bad_stn = [] #stations to toss\n", - "use_lats = [] \n", - "use_lons = []\n", + "# Empty dicts and lists to store GNSS data\n", + "gnss_stns = {}\n", + "bad_stns = {}\n", + "\n", + "for site_name in site_names:\n", + " gnss_stn = gnss.get_gnss_class(gnss_source)(site = site_name)\n", + " gnss_stn.open(print_msg=False)\n", "\n", - "for counter, stn in enumerate(site_names):\n", - " gnss_obj = gnss.get_gnss_class(gnss_source)(site = stn, data_dir = os.path.join(mintpy_dir,'GNSS'))\n", - " gnss_obj.open(print_msg=False)\n", - " \n", " # count number of dates in time range\n", - " dates = gnss_obj.dates\n", + " dates = gnss_stn.dates\n", " range_days = (end_date_gnss - start_date_gnss).days\n", " gnss_count = np.histogram(dates, bins=[start_date_gnss, end_date_gnss])\n", " gnss_count = int(gnss_count[0][0])\n", "\n", - " # for this quick screening check of data quality, we use constant incidence and azimuth angles \n", " # get standard deviation of residuals to linear fit\n", - " disp_los = ut.enu2los(gnss_obj.dis_e, gnss_obj.dis_n, gnss_obj.dis_u, inc_angle, az_angle)\n", + " disp_los = ut.enu2los(gnss_stn.dis_e, gnss_stn.dis_n, gnss_stn.dis_u, inc_angle, az_angle)\n", " disp_detrended = signal.detrend(disp_los)\n", " stn_stdv = np.std(disp_detrended)\n", - " \n", + "\n", " # select GNSS stations based on data completeness and scatter of residuals\n", " disp_detrended = signal.detrend(disp_los)\n", " if range_days*gnss_completeness_threshold <= gnss_count:\n", " if stn_stdv > gnss_residual_stdev_threshold:\n", - " bad_stn.append(stn)\n", + " bad_stns[site_name] = gnss_stn\n", " else:\n", - " use_stn.append(stn)\n", - " use_lats.append(site_lats[counter])\n", - " use_lons.append(site_lons[counter])\n", + " gnss_stns[site_name] = gnss_stn\n", " else:\n", - " bad_stn.append(stn)\n", - "\n", - "site_names = use_stn\n", - "site_lats = use_lats\n", - "site_lons = use_lons\n", + " bad_stns[site_name] = gnss_stn\n", "\n", "# [optional] manually remove additional stations\n", "gnss_to_remove=[]\n", "\n", - "for i, gnss_site in enumerate(gnss_to_remove):\n", - " if gnss_site in site_names:\n", - " site_names.remove(gnss_site)\n", - " if gnss_site not in bad_stn:\n", - " bad_stn.append(gnss_site)\n", + "for site_name in gnss_to_remove:\n", + " bad_stns[site_name] = gnss_stns[site_name]\n", + " del gnss_stns[site_name]\n", + "\n", + "# Final list of site names\n", + "site_names = [*gnss_stns]\n", + "bad_site_names = [*bad_stns]\n", "\n", "print(\"\\nFinal list of {} stations used in analysis:\".format(len(site_names)))\n", "print(site_names)\n", - "print(\"List of {} stations removed from analysis\".format(len(bad_stn)))\n", - "print(bad_stn)" + "print(\"List of {} stations removed from analysis\".format(len(bad_site_names)))\n", + "print(bad_site_names)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Plot stations in map view\n", + "vmin, vmax = -25, 25\n", + "cmap = plt.get_cmap('RdBu_r')\n", + "\n", + "fig, ax = plt.subplots(figsize=[18, 5.5])\n", + "cax = ax.imshow(insar_velocities, cmap=cmap, vmin=vmin, vmax=vmax, interpolation='nearest', extent=(W, E, S, N))\n", + "cbar = fig.colorbar(cax, ax=ax)\n", + "cbar.set_label('LOS velocity [mm/year]')\n", + "\n", + "# Plot valid stations\n", + "for i, gnss_stn in enumerate(gnss_stns.values()):\n", + " # Retrieve station information\n", + " gnss_stn.get_site_lat_lon()\n", + "\n", + " # Plot station\n", + " ax.scatter(gnss_stn.site_lon, gnss_stn.site_lat, s=8**2, c='k',\n", + " label=('valid' if i == 0 else None))\n", + " ax.annotate(gnss_stn.site, (gnss_stn.site_lon, gnss_stn.site_lat))\n", + "\n", + "# Plot rejected stations\n", + "for i, gnss_stn in enumerate(bad_stns.values()):\n", + " # Retrieve station information\n", + " gnss_stn.get_site_lat_lon()\n", + "\n", + " # Plot station\n", + " ax.scatter(gnss_stn.site_lon, gnss_stn.site_lat, s=8**2,\n", + " marker='X', facecolor='none', edgecolor='k',\n", + " label=('rejected' if i == 0 else None))\n", + " ax.annotate(gnss_stn.site, (gnss_stn.site_lon, gnss_stn.site_lat))\n", + "\n", + "ax.legend()\n" ] }, { @@ -749,14 +809,12 @@ "outputs": [], "source": [ "# Empty dicts and lists to store GNSS data\n", - "gnss_stns = {}\n", "gnss_velocities = []\n", "gnss_velocity_errors = []\n", "\n", "# Loop through valid GNSS stations\n", - "for i, site_name in enumerate(site_names):\n", + "for i, gnss_stn in enumerate(gnss_stns.values()):\n", " # Retrieve station information\n", - " gnss_stn = gnss.get_gnss_class(gnss_source)(site = site_name)\n", " gnss_stn.get_los_displacement(insar_metadata,\n", " start_date=start_date,\n", " end_date=end_date)\n", @@ -768,9 +826,6 @@ " # Outlier detection and removal\n", " remove_gnss_outliers(gnss_stn, 'LOS', model=ts_functions,\n", " threshold=3, max_iter=2, verbose=False)\n", - "\n", - " # Record filtered station to dictionary\n", - " gnss_stns[site_name] = gnss_stn\n", " \n", " # Model outlier-removed time-series\n", " dis_los_hat, m_hat, mhat_se = model_gnss_timeseries(gnss_stn, 'LOS', ts_functions)\n", @@ -845,7 +900,10 @@ "source": [ "# reference GNSS stations to GNSS reference site\n", "print(\"Using GNSS reference station: \", sitedata['sites'][site]['gps_ref_site_name'])\n", - "ref_site_ind = site_names.index(sitedata['sites'][site]['gps_ref_site_name'])\n", + "if sitedata['sites'][site]['gps_ref_site_name'] == 'auto':\n", + " ref_site_ind = 0\n", + "else:\n", + " ref_site_ind = site_names.index(sitedata['sites'][site]['gps_ref_site_name'])\n", "gnss_velocities = gnss_velocities - gnss_velocities[ref_site_ind]\n", "\n", "# reference InSAR to GNSS reference site\n",