Merge branch 'main' of github.com:brendanjmeade/celeri

celeri/celeri.py

@@ -3089,7 +3089,6 @@ def get_strain_rate_displacements(
     lon_obs = np.deg2rad(lon_obs)
     lat_obs = latitude_to_colatitude(lat_obs)
     lat_obs = np.deg2rad(lat_obs)
-
     # Calculate displacements from homogeneous strain
     u_up = np.zeros(
         lon_obs.size
@@ -3118,7 +3117,6 @@ def get_block_strain_rate_to_velocities_partials(block, station, segment):
             block_centroid_lon, block_centroid_lat = get_block_centroid(
                 segment, strain_rate_block_idx[i]
             )
-
             # Find stations on current block
             station_idx = np.where(station.block_label == strain_rate_block_idx[i])[0]
             stations_block_lon = station.lon[station_idx].to_numpy()
@@ -3165,17 +3163,17 @@ def get_block_strain_rate_to_velocities_partials(block, station, segment):
                 strain_rate_lon_lat=1,
             )
             block_strain_rate_operator[3 * station_idx, 3 * i] = vel_east_lon_lon
-            block_strain_rate_operator[3 * station_idx, 3 * i + 1] = vel_east_lat_lat
+            block_strain_rate_operator[3 * station_idx, 3 * i + 1] = -vel_east_lat_lat
             block_strain_rate_operator[3 * station_idx, 3 * i + 2] = vel_east_lon_lat
             block_strain_rate_operator[3 * station_idx + 1, 3 * i] = vel_north_lon_lon
             block_strain_rate_operator[3 * station_idx + 1, 3 * i + 1] = (
-                vel_north_lat_lat
+                -vel_north_lat_lat
             )
             block_strain_rate_operator[3 * station_idx + 1, 3 * i + 2] = (
                 vel_north_lon_lat
             )
             block_strain_rate_operator[3 * station_idx + 2, 3 * i] = vel_up_lon_lon
-            block_strain_rate_operator[3 * station_idx + 2, 3 * i + 1] = vel_up_lat_lat
+            block_strain_rate_operator[3 * station_idx + 2, 3 * i + 1] = -vel_up_lat_lat
             block_strain_rate_operator[3 * station_idx + 2, 3 * i + 2] = vel_up_lon_lat
     return block_strain_rate_operator, strain_rate_block_idx
 
@@ -4840,9 +4838,7 @@ def write_output(
 
     # Write all major variables to .pkl file in output folder
     if bool(command.pickle_save):
-        with open(
-            os.path.join(command.output_path, "output" + ".pkl"), "wb"
-        ) as f:
+        with open(os.path.join(command.output_path, "output" + ".pkl"), "wb") as f:
             pickle.dump([command, estimation, station, segment, block, meshes], f)
 
 
@@ -6003,13 +5999,14 @@ def latitude_to_colatitude(lat):
     """
     if lat.size == 1:  # Deal with the scalar case
         if lat >= 0:
-            lat = 90.0 - lat
+            colat = 90.0 - lat
         elif lat < 0:
-            lat = -90.0 - lat
+            colat = -90.0 - lat
     else:  # Deal with the array case
-        lat[np.where(lat >= 0)[0]] = 90.0 - lat[np.where(lat >= 0)[0]]
-        lat[np.where(lat < 0)[0]] = -90.0 - lat[np.where(lat < 0)[0]]
-    return lat
+        colat = np.zeros_like(lat)
+        colat[np.where(lat >= 0)[0]] = 90.0 - lat[np.where(lat >= 0)[0]]
+        colat[np.where(lat < 0)[0]] = -90.0 - lat[np.where(lat < 0)[0]]
+    return colat
 
 
 def get_transverse_projection(lon0, lat0):

notebooks/celeri_western_north_america_dense_tricoup.ipynb

@@ -384,6 +384,52 @@
     "# print(estimation.tensile_slip_rates[csi])\n",
     "print(segment.azimuth[csi])\n"
    ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Debugging strain rate calculations"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "strain_rate_block_index[0]\n",
+    "\n",
+    "centroid_lon = 133.3318559095709\n",
+    "centroid_lat = 34.42960594634139\n",
+    "\n",
+    "lonrange = np.arange(centroid_lon - 1, centroid_lon + 1, 0.05)\n",
+    "latrange = np.arange(centroid_lat - 1, centroid_lat + 1, 0.05)\n",
+    "\n",
+    "fslon, fslat = np.meshgrid(lonrange, latrange)\n",
+    "\n",
+    "fakesta = pd.DataFrame(columns=[\"lon\", \"lat\", \"block_label\"])\n",
+    "fakesta.lon = np.reshape(fslon, -1)\n",
+    "fakesta.lat = np.reshape(fslat, -1)\n",
+    "fakesta.block_label = 10*np.ones_like(fakesta.lon)\n",
+    "\n",
+    "strain_partials, _ = celeri.get_block_strain_rate_to_velocities_partials(block, fakesta, segment)\n",
+    "\n",
+    "strain_comp = [0, 1, 0]\n",
+    "\n",
+    "strainvels = strain_partials@strain_comp\n",
+    "\n",
+    "fig, ax = plt.subplots()\n",
+    "ax.quiver(fakesta.lon, fakesta.lat, strainvels[0::3], strainvels[1::3])\n",
+    "plt.show()\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": []
   }
  ],
  "metadata": {