|
627 | 627 | "command = 'timeseries2velocity.py ' + timeseries_filename + function_str\n", |
628 | 628 | "process = subprocess.run(command, shell=True)\n", |
629 | 629 | "\n", |
630 | | - "# load velocity file\n", |
| 630 | + "# Load velocity file\n", |
631 | 631 | "insar_velocities,_ = readfile.read(vel_file, datasetName = 'velocity') # read velocity file\n", |
632 | 632 | "insar_velocities = insar_velocities * 1000. # convert InSAR velocities from m/yr to mm/yr\n", |
633 | 633 | "\n", |
634 | | - "# set masked pixels to NaN\n", |
| 634 | + "# Set masked pixels to NaN\n", |
635 | 635 | "msk,_ = readfile.read(msk_file)\n", |
636 | 636 | "insar_velocities[msk == 0] = np.nan\n", |
637 | 637 | "insar_velocities[insar_velocities == 0] = np.nan" |
|
943 | 943 | }, |
944 | 944 | "outputs": [], |
945 | 945 | "source": [ |
946 | | - "# reference GNSS stations to GNSS reference site\n", |
| 946 | + "# Reference GNSS stations to GNSS reference site\n", |
947 | 947 | "ref_site = sitedata['sites'][site]['gps_ref_site_name']\n", |
948 | 948 | "if sitedata['sites'][site]['gps_ref_site_name'] == 'auto':\n", |
949 | 949 | " ref_site = [*gnss_stns][0]\n", |
950 | 950 | "\n", |
951 | 951 | "print(\"Using GNSS reference station: \", ref_site)\n", |
952 | 952 | "\n", |
953 | | - "# re-reference GNSS sites\n", |
| 953 | + "# Re-reference GNSS sites\n", |
954 | 954 | "if 'vel_mm_yr' in gnss_stns[ref_site].vel_dict:\n", |
955 | 955 | " ref_vel = gnss_stns[ref_site].vel_dict['vel_mm_yr']\n", |
956 | 956 | "\n", |
|
963 | 963 | " # Remove non-referenced station data to avoid future confusion\n", |
964 | 964 | " del gnss_stn.vel_dict['vel_mm_yr']\n", |
965 | 965 | "\n", |
966 | | - "# reference InSAR to GNSS reference site\n", |
| 966 | + "# Reference InSAR to GNSS reference site\n", |
967 | 967 | "ref_site_lat, ref_site_lon = gnss_stns[ref_site].get_site_lat_lon()\n", |
968 | 968 | "ref_y, ref_x = ut.coordinate(insar_metadata).geo2radar(ref_site_lat, ref_site_lon)[:2]\n", |
969 | 969 | "insar_velocities = insar_velocities - insar_velocities[ref_y, ref_x]\n", |
970 | 970 | "\n", |
971 | | - "# plot GNSS stations on InSAR velocity field\n", |
| 971 | + "# Plot GNSS stations on InSAR velocity field\n", |
972 | 972 | "vmin, vmax = -25, 25\n", |
973 | 973 | "cmap = plt.get_cmap('RdBu_r')\n", |
974 | 974 | "\n", |
|
1028 | 1028 | }, |
1029 | 1029 | "outputs": [], |
1030 | 1030 | "source": [ |
1031 | | - "#Set Parameters\n", |
| 1031 | + "# Set Parameters\n", |
1032 | 1032 | "pixel_radius = 5 #number of InSAR pixels to average for comparison with GNSS\n", |
1033 | 1033 | "\n", |
1034 | | - "#Loop over GNSS station locations\n", |
| 1034 | + "# Loop over GNSS station locations\n", |
1035 | 1035 | "for site_name in site_names:\n", |
1036 | 1036 | " gnss_stn = gnss_stns[site_name]\n", |
1037 | 1037 | " \n", |
1038 | | - " # convert GNSS station lat/lon information to InSAR x/y grid\n", |
| 1038 | + " # Convert GNSS station lat/lon information to InSAR x/y grid\n", |
1039 | 1039 | " stn_lat, stn_lon = gnss_stn.get_site_lat_lon()\n", |
1040 | 1040 | " stn_y, stn_x = ut.coordinate(insar_metadata).geo2radar(stn_lat, stn_lon)[:2]\n", |
1041 | 1041 | " \n", |
1042 | | - " # get velocities and residuals\n", |
| 1042 | + " # Get velocities and residuals\n", |
1043 | 1043 | " gnss_site_vel = gnss_stn.vel_dict['gnss_site_vel']\n", |
1044 | | - " #Caution: If you expand the radius parameter farther than the bounding grid it will break. \n", |
1045 | | - " #To fix, remove the station in section 4 when the site_names list is filtered\n", |
| 1044 | + " # Caution: If you expand the radius parameter farther than the bounding grid it will break. \n", |
| 1045 | + " # To fix, remove the station in section 4 when the site_names list is filtered\n", |
1046 | 1046 | " vel_px_rad = insar_velocities[stn_y-pixel_radius:stn_y+1+pixel_radius, \n", |
1047 | 1047 | " stn_x-pixel_radius:stn_x+1+pixel_radius]\n", |
1048 | 1048 | " insar_site_vel = np.median(vel_px_rad)\n", |
1049 | 1049 | " residual = gnss_site_vel - insar_site_vel\n", |
1050 | 1050 | "\n", |
1051 | | - " # populate data structure\n", |
| 1051 | + " # Populate data structure\n", |
1052 | 1052 | " gnss_stn.vel_dict['insar_site_vel'] = insar_site_vel\n", |
1053 | 1053 | " gnss_stn.vel_dict['residual'] = residual\n", |
1054 | 1054 | "\n", |
|
1076 | 1076 | "stn_dist_list = []\n", |
1077 | 1077 | "n_gnss_sites = len(site_names)\n", |
1078 | 1078 | "\n", |
1079 | | - "# loop over stations\n", |
| 1079 | + "# Loop over stations\n", |
1080 | 1080 | "for i in range(n_gnss_sites - 1):\n", |
1081 | 1081 | " stn1 = gnss_stns[site_names[i]]\n", |
1082 | 1082 | " for j in range(i + 1, n_gnss_sites):\n", |
|
1195 | 1195 | "metadata": {}, |
1196 | 1196 | "outputs": [], |
1197 | 1197 | "source": [ |
1198 | | - "# use the assumed non-earthquake displacement as the insar_displacment for statistics and convert to mm\n", |
| 1198 | + "# Use the assumed non-earthquake displacement as the insar_displacment for statistics and convert to mm\n", |
1199 | 1199 | "insar_velocities,_ = readfile.read(vel_file, datasetName = 'velocity') #read velocity\n", |
1200 | 1200 | "velStart = sitedata['sites'][site]['download_start_date']\n", |
1201 | 1201 | "insar_velocities = insar_velocities * 1000. # convert velocities from m to mm\n", |
1202 | 1202 | "\n", |
1203 | | - "# set masked pixels to NaN\n", |
| 1203 | + "# Set masked pixels to NaN\n", |
1204 | 1204 | "msk,_ = readfile.read(msk_file)\n", |
1205 | 1205 | "insar_velocities[msk == 0] = np.nan\n", |
1206 | 1206 | "insar_velocities[insar_velocities == 0] = np.nan\n", |
1207 | 1207 | "\n", |
1208 | | - "# display map of velocity data after masking\n", |
| 1208 | + "# Display map of velocity data after masking\n", |
1209 | 1209 | "cmap = plt.get_cmap('RdBu_r')\n", |
1210 | 1210 | "fig, ax = plt.subplots(figsize=[18, 5.5])\n", |
1211 | 1211 | "img1 = ax.imshow(insar_velocities, cmap=cmap, vmin=-20, vmax=20, interpolation='nearest', extent=(W, E, S, N))\n", |
|
1438 | 1438 | }, |
1439 | 1439 | "outputs": [], |
1440 | 1440 | "source": [ |
1441 | | - "# grab the time-series file used for time function estimation given the template setup\n", |
| 1441 | + "# Grab the time-series file used for time function estimation given the template setup\n", |
1442 | 1442 | "template = readfile.read_template(os.path.join(mintpy_dir, 'smallbaselineApp.cfg'))\n", |
1443 | 1443 | "template = ut.check_template_auto_value(template)\n", |
1444 | 1444 | "insar_ts_file = TimeSeriesAnalysis.get_timeseries_filename(template, mintpy_dir)['velocity']['input']\n", |
1445 | 1445 | "\n", |
1446 | | - "# read the time-series file\n", |
| 1446 | + "# Read the time-series file\n", |
1447 | 1447 | "insar_ts, atr = readfile.read(insar_ts_file, datasetName='timeseries')\n", |
1448 | 1448 | "mask = readfile.read(os.path.join(mintpy_dir, 'maskTempCoh.h5'))[0]\n", |
1449 | 1449 | "print(f'reading timeseries from file: {insar_ts_file}')\n", |
|
1454 | 1454 | "date0, date1 = date_list[0], date_list[-1]\n", |
1455 | 1455 | "insar_dates = ptime.date_list2vector(date_list)[0]\n", |
1456 | 1456 | "\n", |
1457 | | - "# spatial reference\n", |
| 1457 | + "# Spatial reference\n", |
1458 | 1458 | "coord = ut.coordinate(atr)\n", |
1459 | 1459 | "ref_site = sitedata['sites'][site]['gps_ref_site_name']\n", |
1460 | 1460 | "ref_gnss_obj = gnss_stns[ref_site]\n", |
|
1470 | 1470 | "for i, site_name in enumerate(site_names):\n", |
1471 | 1471 | " prog_bar.update(i+1, suffix=f'{site_name} {i+1}/{num_site}')\n", |
1472 | 1472 | "\n", |
1473 | | - " ## read data\n", |
1474 | | - " # recall gnss station displacements with outliers removed\n", |
| 1473 | + " ## Read data\n", |
| 1474 | + " # Recall gnss station displacements with outliers removed\n", |
1475 | 1475 | " gnss_obj = gnss_stns[site_name]\n", |
1476 | 1476 | " gnss_lalo = (gnss_obj.site_lat, gnss_obj.site_lon)\n", |
1477 | 1477 | "\n", |
1478 | | - " # get relative LOS displacement on common dates\n", |
| 1478 | + " # Get relative LOS displacement on common dates\n", |
1479 | 1479 | " gnss_dates = np.array(sorted(list(set(gnss_obj.dates) & set(ref_gnss_obj.dates))))\n", |
1480 | 1480 | " gnss_dis = np.zeros(gnss_dates.shape, dtype=np.float32)\n", |
1481 | 1481 | " for i, date_i in enumerate(gnss_dates):\n", |
1482 | 1482 | " idx1 = np.where(gnss_obj.dates == date_i)[0][0]\n", |
1483 | 1483 | " idx2 = np.where(ref_gnss_obj.dates == date_i)[0][0]\n", |
1484 | 1484 | " gnss_dis[i] = gnss_obj.dis_los[idx1] - ref_gnss_obj.dis_los[idx2]\n", |
1485 | 1485 | " \n", |
1486 | | - " # shift GNSS to zero-mean in time [for plotting purpose]\n", |
| 1486 | + " # Shift GNSS to zero-mean in time [for plotting purpose]\n", |
1487 | 1487 | " gnss_dis -= np.nanmedian(gnss_dis)\n", |
1488 | 1488 | "\n", |
1489 | | - " # read InSAR\n", |
| 1489 | + " # Read InSAR\n", |
1490 | 1490 | " y, x = coord.geo2radar(gnss_lalo[0], gnss_lalo[1])[:2]\n", |
1491 | 1491 | " insar_dis = insar_ts[:, y, x] - ref_insar_dis\n", |
1492 | | - " # apply a constant shift in time to fit InSAR to GNSS\n", |
| 1492 | + " # Apply a constant shift in time to fit InSAR to GNSS\n", |
1493 | 1493 | " comm_dates = sorted(list(set(gnss_dates) & set(insar_dates)))\n", |
1494 | 1494 | " if comm_dates:\n", |
1495 | 1495 | " insar_flag = [x in comm_dates for x in insar_dates]\n", |
1496 | 1496 | " gnss_flag = [x in comm_dates for x in gnss_dates]\n", |
1497 | 1497 | " insar_dis -= np.nanmedian(insar_dis[insar_flag] - gnss_dis[gnss_flag])\n", |
1498 | 1498 | "\n", |
1499 | | - " ## plot figure\n", |
| 1499 | + " ## Plot figure\n", |
1500 | 1500 | " if gnss_dis.size > 0 and np.any(~np.isnan(insar_dis)):\n", |
1501 | 1501 | " fig, ax = plt.subplots(figsize=(10, 3))\n", |
1502 | 1502 | " ax.axhline(color='grey',linestyle='dashed', linewidth=2)\n", |
|
0 commit comments