changed terrain interpolation method to avoid edge artifacts (hopefully)

arcgis_terrain.py

@@ -10,6 +10,7 @@ import math
 from pyproj import Transformer, Proj, transform
 from scipy.interpolate import griddata
 from scipy import interpolate
+from scipy import ndimage
 import time
 
 # This function computes the factor of the argument passed
@@ -208,6 +209,16 @@ def get_terrain_map(lat_lon = [0,0], sample_dist = 10, extent = 100, heading = 0
         y = np.append(y, y_row, axis=0)
         e = np.append(e, e_row, axis=0)
 
+    # flip elevation data up to down to match other layers
+    e = np.flipud(e)
+    # interpolate terrain to match size/resolution of other layers
+    c = np.int(extent/sample_dist)
+    scale_factor = 3/20 # factor to get 6.66667m mapping from 1m mapping (1/6.6667)
+    scaled_extent = np.ceil(scale_factor*extent).astype(np.int)
+
+    factor = scaled_extent/e.shape[0]
+    e_interp = ndimage.zoom(e, factor, order = 3)
+
     if show_plot:
         # Attaching 3D axis to the figure
         grid_x, grid_y = np.mgrid[min(x):max(x):100j, min(y):max(y):100j]
@@ -228,7 +239,7 @@ def get_terrain_map(lat_lon = [0,0], sample_dist = 10, extent = 100, heading = 0
         print("y max: {}".format(np.max(y)))
         print("y max - min: {}".format(np.max(y)-np.min(y)))
 
-    return [e,x,y,data,lat_lon]
+    return [e,e_interp,x,y,data,lat_lon]
     # savedimg = elv_layer.export_image(bbox=elv_layer.extent, size=[3840,2160], f='image', save_folder='.', save_file='testerino.jpg')
 
 if __name__ == "__main__":
@@ -237,36 +248,51 @@ if __name__ == "__main__":
     # anchor_point = [float(ics_pt[0]), float(ics_pt[1])]
     anchor_point = [42.17965, -74.21362]
     extent = 20e3
-    sample_dist = int(extent/100)
+    sample_dist = int(extent/10)
     heading = 0
     start_time = time.time()
     # save terrain as csv file (this method is pretty slow, but can compensate with interp)
-    [e,x,y,data,ll_pt] = get_terrain_map(lat_lon=anchor_point,
+    [e,e_interp,x,y,data,ll_pt] = get_terrain_map(lat_lon=anchor_point,
                                          sample_dist = sample_dist,
                                          extent = extent,
                                          heading = -heading) # because flipping
-    # flip elevation data up to down to match other layers
-    e = np.flipud(e)
 
-    scale_factor = 3/20 # factor to get 6.66667m mapping from 1m mapping (1/6.6667)
-    scaled_extent = np.ceil(scale_factor*extent).astype(np.int)
+    plt.imshow(e_interp)
+    plt.title('e interp')
+    plt.show()
 
-    # interpolate terrain to match size/resolution of other layers
-    c = np.int(extent/sample_dist)
-    f = interpolate.interp2d(np.linspace(0, extent, c), np.linspace(0, extent, c), e, kind='cubic')
+    gdt = np.gradient(e_interp)
+    plt.imshow(np.sqrt(gdt[0]**2 + gdt[1]**2))
+    plt.title('e interp grad')
+    plt.show()
+
+    c = np.int(extent / sample_dist)
+    scale_factor = 3 / 20  # factor to get 6.66667m mapping from 1m mapping (1/6.6667)
+    scaled_extent = np.ceil(scale_factor * extent).astype(np.int)
+    x_start = np.linspace(0, extent, c)
+    y_start = np.linspace(0, extent, c)
+    X_start, Y_start = np.meshgrid(x_start, y_start)
+    Z_start = np.zeros_like(X_start)
+
+    f = interpolate.Rbf(X_start, Y_start, Z_start, e, function='multiquadric')
     x_temp = np.linspace(0,extent,scaled_extent) # get correct size of terrain map
     y_temp = np.linspace(0,extent,scaled_extent)
-    e_interp = f(x_temp, y_temp)
+    X_temp, Y_temp = np.meshgrid(x_temp, y_temp)
+    Z_temp = np.zeros_like(X_temp)
+    e_interp = f(X_temp, Y_temp, Z_temp)
 
     plt.imshow(e_interp)
+    plt.title('e interp old')
     plt.show()
-    plt.imshow(e)
-    plt.show()
-    plt.imshow(x)
-    plt.show()
-    plt.imshow(y)
+
+    gdt = np.gradient(e_interp)
+    plt.imshow(np.sqrt(gdt[0]**2 + gdt[1]**2))
+    plt.title('e interp grad old')
     plt.show()
 
+
+    print('done')
+
     # elv_filename = "map_layers\\elv_data_"+file_id+".csv"
     # if save_files:
     #     np.savetxt(elv_filename,e_interp,delimiter=",", fmt='%f')

feature_set.py

@@ -71,8 +71,9 @@ def grab_features(anchor_point, extent, sample_dist = 10, heading = 0, save_file
                 print("querying {} layer...".format(name_list[i]))
                 q = lyr.query(return_count_only=False, return_ids_only=False, return_geometry=True,
                               out_sr='3857', geometry_filter=geom_filter)
-            except json.decoder.JSONDecodeError as e:
-                query_cnt = query_cnt + 1
+            except (json.decoder.JSONDecodeError, TypeError) as e:
+                if type(e) != TypeError:
+                    query_cnt = query_cnt + 1
                 print("error on query: {}".format(e))
                 print("{} layer failed on query, trying again ...".format(name_list[i]))
                 gis = GIS(username="larkinheintzman",password="Meepp97#26640") # linked my arcgis pro account
@@ -218,44 +219,6 @@ def grab_features(anchor_point, extent, sample_dist = 10, heading = 0, save_file
                         inac_bin_map[pts_inac[:,1], pts_inac[:,0]] = 1  # set indices to 1
                         # print("looped(tm) inac calculation time = {} sec".format(time.time() - s_time))
 
-
-                        #-----------------------------
-
-                        # # print("points not looped")
-                        # # do boundary calculation for binary matrix (slow for large bounaries but whatever)
-                        #
-                        # # test_pts is the rectangular matrix covering ring for boundary calculation
-                        # x_test, y_test = np.meshgrid(np.arange(np.min(x_pts), np.max(x_pts), 1) , np.arange(np.min(y_pts), np.max(y_pts), 1))
-                        # test_pts = np.array([x_test.flatten(), y_test.flatten()]).T
-                        # mask = np.zeros(test_pts.shape[0])
-                        # core_pts = np.stack([x_pts,y_pts]).T
-                        #
-                        # for pt in core_pts:
-                        #     test_dists = np.sqrt(np.square(test_pts[:,0] - pt[0]) +
-                        #                          np.square(test_pts[:,1] - pt[1]))
-                        #     mask[test_dists<=1] = 1
-                        #
-                        # # instead of filling gaps, we want to save filled in areas separately
-                        # # so we need to re-create the bin_map here but on inac. points
-                        # x_pts_inac = test_pts[np.where(mask),0].flatten()
-                        # y_pts_inac = test_pts[np.where(mask),1].flatten()
-                        # pts_inac = np.stack([x_pts_inac,y_pts_inac]).T
-                        #
-                        # # remove points being used as linear features
-                        # for pt in core_pts:
-                        #     pts_inac = np.delete(pts_inac, np.where(np.equal(pt,pts_inac).all(1)), axis = 0)
-                        #
-                        # # binarization step
-                        # pts_inac = np.round(pts_inac).astype(np.int)
-                        # # flip y axis
-                        # pts_inac[:,1] = inac_bin_map.shape[1] - pts_inac[:,1]
-                        # # remove any points outside limits of binary map (fixes round versus ceil issues)
-                        # rm_mask = np.logical_or(pts_inac[:,0] < 0, pts_inac[:,0] >= inac_bin_map.shape[1])
-                        # rm_mask = np.logical_or(rm_mask, np.logical_or(pts_inac[:,1] < 0, pts_inac[:,1] >= inac_bin_map.shape[0]))
-                        # pts_inac = pts_inac[np.invert(rm_mask),:]
-                        # inac_bin_map[pts_inac[:,1], pts_inac[:,0]] = 1  # set indices to 1
-                        # print("non looped inac calculation time = {} sec".format(time.time() - s_time))
-
                 # binarization step
                 x_pts_idx = np.round(x_pts).astype(np.int)
                 y_pts_idx = np.round(y_pts).astype(np.int)
@@ -284,19 +247,11 @@ def grab_features(anchor_point, extent, sample_dist = 10, heading = 0, save_file
             np.savetxt(fn,bin_map,delimiter=",", fmt='%f')
 
     # save terrain as csv file (this method is pretty slow, but can compensate with interp)
-    [e,x,y,data,ll_pt] = get_terrain_map(lat_lon=anchor_point,
+    [e,e_interp,x,y,data,ll_pt] = get_terrain_map(lat_lon=anchor_point,
                                          sample_dist = sample_dist,
                                          extent = extent,
                                          heading = -heading) # because flipping
-    # flip elevation data up to down to match other layers
-    e = np.flipud(e)
 
-    # interpolate terrain to match size/resolution of other layers
-    c = np.int(extent/sample_dist)
-    f = interpolate.interp2d(np.linspace(0, extent, c), np.linspace(0, extent, c), e, kind='cubic')
-    x_temp = np.linspace(0,extent,scaled_extent) # get correct size of terrain map
-    y_temp = np.linspace(0,extent,scaled_extent)
-    e_interp = f(x_temp, y_temp)
 
     elv_filename = "map_layers\\elv_data_"+file_id+".csv"
     if save_files:
@@ -352,7 +307,7 @@ if __name__ == "__main__":
         anchor_point = [float(ics_pt[0]), float(ics_pt[1])]
         extent = 20e3
         save_flag = True
-        plot_flag = True
+        plot_flag = False
         file_extension = 'temp'
 
         sample_dist = int(extent/100)