Initial commit, copying from branch unity_test of Larkin's ags_grabber repo

arcgis_terrain.py

@@ -104,9 +104,11 @@ def get_terrain_map(lat_lon = [0,0], sample_dist = 10, extent = 100, heading = 0
 
     # gis = GIS("pro")
     # gis = GIS(url="http://virginiatech.maps.arcgis.com", client_id="rluxzSWjZS6TfeXs", username="hlarkin3_virginiatech", password="arcgisheintzman97#26640", verify_cert=False)
-    gis = GIS(username="larkinheintzman",password="Meepp97#26640")
+    #gis = GIS(username="larkinheintzman",password="Meepp97#26640")
+    gis = GIS(api_key="AAPKe7cb1ab4f2ba44748cf53ac4f30d0caavu5t_-uMdWP8SRjtoea3s66-cRyFZmMTNx4rqy2w5sjwegU_tbJyBesd0LCMmdtV")
     if verbosity:
-        print("Successfully logged in as: " + gis.properties.user.username)
+        #print("Successfully logged in as: " + gis.properties.user.username)
+        print("Logged into GIS")
     elv_map = gis.content.get('58a541efc59545e6b7137f961d7de883')
     elv_layer = elv_map.layers[0]
 
@@ -212,18 +214,21 @@ def get_terrain_map(lat_lon = [0,0], sample_dist = 10, extent = 100, heading = 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)
+    c = np.int32(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)
 
+    print(show_plot)
     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]
+        grid_x, grid_y = np.mgrid[np.min(x):np.max(x):100j, np.min(y):np.max(y):100j]
         points = np.array([x, y])
-        grid_z = griddata(points.transpose(), z, (grid_x, grid_y), method='cubic', fill_value=np.mean(z))
+        print(np.shape(points))
+        #grid_z = griddata(points.transpose(), e, (grid_x, grid_y), method='cubic', fill_value=np.mean(e))
+        grid_z = griddata(points.transpose(), e, (grid_x, grid_y), method='cubic', fill_value=np.mean(e))
         print(grid_z)
         fig = plt.figure()
         ax = fig.gca(projection='3d')
@@ -255,7 +260,9 @@ if __name__ == "__main__":
     [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
+                                         heading = -heading,
+                                         show_plot = True,
+                                         verbosity = True) # because flipping
 
     plt.imshow(e_interp)
     plt.title('e interp')
@@ -296,4 +303,3 @@ if __name__ == "__main__":
     # elv_filename = "map_layers\\elv_data_"+file_id+".csv"
     # if save_files:
     #     np.savetxt(elv_filename,e_interp,delimiter=",", fmt='%f')
-

feature_set.py

@@ -46,12 +46,12 @@ def grab_features(anchor_point, extent, sample_dist = 10, case_name = 'blah', he
     name_list = ['rivers_bdd', 'rivers', 'roads', 'lakes', 'powerlines', 'railroads', 'trails']
     inac_layers = ['rivers_bdd', 'lakes']
 
-    gis = GIS(username="larkinheintzman",password="Meepp97#26640") # linked my arcgis pro account
+    gis = GIS(api_key="AAPKe7cb1ab4f2ba44748cf53ac4f30d0caavu5t_-uMdWP8SRjtoea3s66-cRyFZmMTNx4rqy2w5sjwegU_tbJyBesd0LCMmdtV") # linked Bryson's pro account
     ap_meters = lat_lon2meters(anchor_point[0], anchor_point[1])
 
 
     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)
+    scaled_extent = np.ceil(scale_factor*extent).astype(np.int32)
     viz_cnt = 0
     viz_map = np.zeros([scaled_extent,scaled_extent,len(name_list)+len(inac_layers)])
 
@@ -90,7 +90,7 @@ def grab_features(anchor_point, extent, sample_dist = 10, case_name = 'blah', he
                 print("error on query: {}".format(e))
                 print("{} layer failed on query, trying again ...".format(name_list[i]))
                 # gis.
-                gis = GIS(username="larkinheintzman",password="Meepp97#26640") # linked my arcgis pro account
+                gis = GIS(api_key="AAPKe7cb1ab4f2ba44748cf53ac4f30d0caavu5t_-uMdWP8SRjtoea3s66-cRyFZmMTNx4rqy2w5sjwegU_tbJyBesd0LCMmdtV") # linked my arcgis pro account
                 lyr = FeatureLayer(url=url, gis=gis)
 
         print("query time {}".format(query_endtime - query_starttime))
@@ -99,16 +99,16 @@ def grab_features(anchor_point, extent, sample_dist = 10, case_name = 'blah', he
             print("{} layer failed too many times, leaving empty".format(name_list[i]))
             if save_files:
                 if save_to_folder:
-                    fn = "map_layers\\" + case_name + "\\"+name_list[i]+"_data_"+file_id+".csv"
+                    fn = "./map_layers/" + case_name + "/"+name_list[i]+"_data_"+file_id+".csv"
                     np.savetxt(fn,bin_map,delimiter=",", fmt='%f')
                     if name_list[i] in inac_layers:
-                        fn = "map_layers\\" + case_name + "\\" + name_list[i] + "_inac_data_" + file_id + ".csv"
+                        fn = "./map_layers/" + case_name + "/" + name_list[i] + "_inac_data_" + file_id + ".csv"
                         np.savetxt(fn, bin_map, delimiter=",", fmt='%f')
                 else:
-                    fn = "map_layers\\" + name_list[i]+"_data_"+file_id+".csv"
+                    fn = "./map_layers/" + name_list[i]+"_data_"+file_id+".csv"
                     np.savetxt(fn,bin_map,delimiter=",", fmt='%f')
                     if name_list[i] in inac_layers:
-                        fn = "map_layers\\" + name_list[i] + "_inac_data_" + file_id + ".csv"
+                        fn = "./map_layers/" + name_list[i] + "_inac_data_" + file_id + ".csv"
                         np.savetxt(fn, bin_map, delimiter=",", fmt='%f')
             continue
         print("{} layer sucessfully queried".format(name_list[i]))
@@ -213,7 +213,7 @@ def grab_features(anchor_point, extent, sample_dist = 10, case_name = 'blah', he
                     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)
+                pts_inac = np.round(pts_inac).astype(np.int32)
                 # 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)
@@ -226,8 +226,8 @@ def grab_features(anchor_point, extent, sample_dist = 10, case_name = 'blah', he
             # trim features to scaled extent
             [x_pts, y_pts] = trim_extent(x_pts, y_pts, scaled_extent)
             # binarization step
-            x_pts_idx = np.round(x_pts).astype(np.int)
-            y_pts_idx = np.round(y_pts).astype(np.int)
+            x_pts_idx = np.round(x_pts).astype(np.int32)
+            y_pts_idx = np.round(y_pts).astype(np.int32)
 
             # flip y axis because indices are fliped
             y_pts_idx = bin_map.shape[1] - y_pts_idx
@@ -246,10 +246,10 @@ def grab_features(anchor_point, extent, sample_dist = 10, case_name = 'blah', he
             viz_cnt = viz_cnt + 1
             if save_files:
                 if save_to_folder:
-                    fn = "map_layers\\" + case_name + "\\" + name_list[i] + "_inac_data_" + file_id + ".csv"
+                    fn = "./map_layers/" + case_name + "/" + name_list[i] + "_inac_data_" + file_id + ".csv"
                     np.savetxt(fn, inac_bin_map, delimiter=",", fmt='%f')
                 else:
-                    fn = "map_layers\\" + name_list[i] + "_inac_data_" + file_id + ".csv"
+                    fn = "./map_layers/" + name_list[i] + "_inac_data_" + file_id + ".csv"
                     np.savetxt(fn, inac_bin_map, delimiter=",", fmt='%f')
 
 
@@ -257,10 +257,10 @@ def grab_features(anchor_point, extent, sample_dist = 10, case_name = 'blah', he
         viz_cnt = viz_cnt + 1
         if save_files:
             if save_to_folder:
-                fn = "map_layers\\" + case_name + "\\" + name_list[i]+"_data_"+file_id+".csv"
+                fn = "./map_layers/" + case_name + "/" + name_list[i]+"_data_"+file_id+".csv"
                 np.savetxt(fn,bin_map,delimiter=",", fmt='%f')
             else:
-                fn = "map_layers\\" + name_list[i]+"_data_"+file_id+".csv"
+                fn = "./map_layers/" + name_list[i]+"_data_"+file_id+".csv"
                 np.savetxt(fn,bin_map,delimiter=",", fmt='%f')
 
 
@@ -272,10 +272,10 @@ def grab_features(anchor_point, extent, sample_dist = 10, case_name = 'blah', he
 
     if save_files:
         if save_to_folder:
-            elv_filename = "map_layers\\" + case_name + "\\elv_data_" + file_id + ".csv"
+            elv_filename = "./map_layers/" + case_name + "/elv_data_" + file_id + ".csv"
             np.savetxt(elv_filename,e_interp,delimiter=",", fmt='%f')
         else:
-            elv_filename = "map_layers\\elv_data_" + file_id + ".csv"
+            elv_filename = "./map_layers/elv_data_" + file_id + ".csv"
             np.savetxt(elv_filename,e_interp,delimiter=",", fmt='%f')
 
     if plot_data:
@@ -406,7 +406,8 @@ if __name__ == "__main__":
         # [48.64606, -122.4247,'[48.64606, -122.4247]']
         ]
 
-    base_dir = 'C:/Users/Larkin/ags_grabber'
+    #You should be running this script from the base dir
+    base_dir = './'
 
     # for item in ics:
     #     try:
@@ -427,14 +428,14 @@ if __name__ == "__main__":
         anchor_point = [float(ics_pt[0]), float(ics_pt[1])]
         extent = 10e3
         save_flag = True
-        plot_flag = False
+        plot_flag = True
         file_extension = 'temp'
 
         sample_dist = int(extent/100)
         heading = 0
         start_time = time.time()
         grab_features(anchor_point = anchor_point, extent = extent, sample_dist = sample_dist, case_name = str(ics_pt[2]),
-                      heading = heading, save_files = save_flag, save_to_folder = False, file_id = file_extension, plot_data = plot_flag)
+                      heading = heading, save_files = save_flag, save_to_folder = True, file_id = file_extension, plot_data = plot_flag)
         time.sleep(1) # wait for... files to settle?
 
         # run matlab

unity_heightmap.py

+import numpy as np
+import matplotlib.pyplot as plt
+from scipy.io import loadmat
+from PIL import Image
+from PIL import features
+#import libtiff
+import cv2
+
+#Loads matlab matrices created by feature_set.py as numpy arrays
+def main():
+
+    #Load the matlab matrix
+    dir = './matlab_data/'
+    savedir = './unity_data/'
+    filename = 'BW_LFandInac_Zelev_kentland.mat'
+    matdict = loadmat(dir+filename, appendmat=False)
+
+    #Testing stuff
+    k = matdict.keys()
+    print(k)
+
+    #These are square numpy arrays of equal size
+    elev = matdict['sZelev']        #Values are height of grid cell (in meters)
+    linfeat = matdict['BWLF']       #Values are 1 for linear feature present, 0 for none
+    obstacles = matdict['BWInac']   #Values are 1 for obstacle, 0 for clear
+
+    #Resize arrays to work with Unity
+    u_elev = clip_size(elev)
+    u_linfeat = clip_size(linfeat)
+
+    check_elev(u_elev)
+
+    make_featmap(u_linfeat,savedir,"clip_featmap.png")
+    make_heightmap(u_elev,savedir,"clip_heightmap.png")
+
+    return
+
+#Unity needs heightmap to be a square of size (power of two)+1
+#Easiest solution is to reduce size of existing data
+def clip_size(mat):
+
+    clip = 1025
+    full = np.shape(mat)[1]
+    diff = (full - clip) / 2
+    diff = np.int32(diff)
+
+    #print(np.shape(mat))
+    #print(diff)
+
+    clipped = mat[diff:diff+clip,diff:diff+clip]
+
+    return clipped
+
+
+#Returns the difference between the min and max elevation.
+#To be used as terrain height in Unity.
+def check_elev(elev):
+
+    diff = np.max(elev)-np.min(elev)
+    print(diff)
+    print("Height of environment is %f\n" % diff)
+
+    return
+
+#Saves image of linear features so it can be used as a mesh in Unity
+def make_featmap(linfeat,savedir,name):
+    path = savedir + name
+
+    #Convert to png format image
+    print(np.shape(linfeat))
+    feat_int8 = (linfeat*255).astype(np.uint8)
+    img = Image.fromarray(feat_int8)
+    img.save(path)
+
+    return
+
+#Converts the terrain height data in the .mat files into a grayscale heightmap
+#that can be loaded into Unity. Heightmaps are 16-bit grayscale images saved in
+#the RAW format.
+def make_heightmap(elev,savedir,name):
+
+    path = savedir + name
+
+
+    #Makes sure libtiff is working right
+    print(features.check('libtiff'))
+    #fig = plt.figure()
+    #plt.imshow(elev)
+    #plt.show()
+
+
+    #Normalize values to be between 0-1 and preview as grayscale image
+    elev_gray = np.float16((elev-np.min(elev))/(np.max(elev)-np.min(elev)))
+    #elev_gray = np.ushort(elev)
+    #elev_gray = np.uint16(elev)
+
+
+    #From array expects uint8
+
+
+    fig = plt.figure()
+    #plt.imshow(elev_gray)
+    #plt.savefig('./unity_data/test.png',pad_inches=0)
+    #plt.show()
+
+
+    #f = open(name,'w')
+
+
+    #Save as png
+    elev_int8 = (elev_gray*255).astype(np.uint8)
+    img = Image.fromarray(elev_int8)
+    img.save(path)
+
+
+    #Attempts to convert to raw image format through code
+    img = img.convert('L')
+    img.save('./unity_data/converted_kentland_gray.raw', 'TIFF', compression='raw')
+    #img = img.convert('1')
+    #img = img.convert('I;16')
+
+    img = img.convert('1')
+    img.save(name, 'TIFF', compression='tiff_raw_16')  #want tiff_raw_16
+    img.show()
+    #cv2.imshow('Kentland Heightmap',elev_gray)
+    #cv2.waitKey(0)
+
+    #cv2.imwrite('kentland_heightmap.raw',elev_gray)
+
+    return
+
+
+
+
+
+
+
+if __name__ == '__main__':
+    main()