Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

GNSS station checks #65

Merged
merged 1 commit into from
Sep 11, 2024
Merged

GNSS station checks #65

Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
130 changes: 94 additions & 36 deletions methods/secular/Secular_Requirement_Validation.ipynb
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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": [
"<div class=\"alert alert-warning\">\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",
"</div>"
]
},
{
"cell_type": "code",
"execution_count": null,
Expand All @@ -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"
]
},
{
Expand Down Expand Up @@ -664,6 +684,13 @@
"### 4.4. Get GNSS Position Time Series <a id='gps_ts'></a>"
]
},
{
"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,
Expand All @@ -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"
]
},
{
Expand Down Expand Up @@ -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",
Expand All @@ -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",
Expand Down Expand Up @@ -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",
Expand Down
Loading