Hikersim correctly iterates through all collected terrain data. Added a search...

Hikersim correctly iterates through all collected terrain data. Added a search incident visualization function too.

lpm_hikersim.py

@@ -4,6 +4,7 @@ import matplotlib.pyplot as plt
 import os, sys, inspect
 from os import path, getcwd
 import pandas as pd
+import time
 
 #Add parent directory to path, where LPM code is
 currentdir = os.path.dirname(os.path.abspath(inspect.getfile(inspect.currentframe())))
@@ -21,6 +22,7 @@ from LostPersonModel.main_hiker import run_replicate   #change this import to wh
 
 def main(exp_name='test', start_idx=0, n_envs=1, n_iter=10):
 
+    #Size in meters of GIS datasets
     test_extent = 20000
     download_extent = 40000
 
@@ -43,40 +45,54 @@ def main(exp_name='test', start_idx=0, n_envs=1, n_iter=10):
     scale_factor = 3/20 # factor to get 6.66667m mapping from 1m mapping (1/6.6667)
     test_extent = np.ceil(scale_factor*test_extent).astype(np.int32)
     download_extent = np.ceil(scale_factor*download_extent).astype(np.int32)
+    half_test = np.ceil(test_extent*0.5).astype(np.int32)
+    half_download = np.ceil(download_extent*0.5).astype(np.int32)
 
-
-    for i in range(start_idx, end_idx):
+    #Iterate through downloaded GIS datasets
+    for env in range(start_idx, end_idx):
         #Pull incident locations for environment
-        i_area = np.argwhere(areas==i)
+        i_area = np.argwhere(areas==env)
         i_area = np.transpose(i_area)[0]
 
-        #test - print start / end locations on grid
-        grid = np.zeros((extent,extent),dtype=np.int8())
+        #Load datasets
+        dirn = './map_layers/' + exp_name + '/' + exp_name + '_' + str(env) + '/'
+        elev = np.load(dirn+'elevation.npy')
+        lf = np.load(dirn+'linear_feats.npy')
+        inac = np.load(dirn+'inac_layers.npy')
 
+        #Iterate through search incidents
         for incident in range(0,len(i_area)):
+            cur_ipp = ipp_latlons[i_area[incident]]
+            ipp_xy = lat_lon2meters(cur_ipp[0],cur_ipp[1])
+            find_ll = find_latlons[i_area[incident]]
+            find_xy = lat_lon2meters(find_ll[0],find_ll[1])
             #Dealing with large areas. Find the IPP point and resize around it
-            if(len(i_area) > 0):
-                for subzone in range(0,len(i_area)):
-                    #First lat/lon in area is at the center
-                    if(subzone == 0):
-                        sub_center = ipp_latlons[i_area[incident]]
-                        center_xy = lat_lon2meters(sub_center[0],sub_center[1])
-                    #Determine position of other IPPs from the first
-                    else:
-                        cur_ipp = ipp_latlons[i_area[incident+subzone]]
-                        ipp_xy = lat_lon2meters(cur_ipp[0],cur_ipp[1])
-                        print(cur_ipp)
-                        #Distance calculation
-                        x_dist = cur_ipp[0] - sub_center[0]
-                        y_dist = cur_ipp[1] - sub_center[1]
-                        x_index = np.ceil(0.5*extent + scale_factor*x_dist).astype(np.int32)
-                        y_index = np.ceil(0.5*extent + scale_factor*y_dist).astype(np.int32)
-                        if(x_index >= extent or y_index >= extent or x_index)
-                        #Grab subset of big grid
-                        sub_grid =
-
-
-
+            if(len(i_area) > 1):
+                #First lat/lon in area is at the center
+                if(incident == 0):
+                    center_xy = ipp_xy
+                #Find array index of IPP
+                x_dist = ipp_xy[0] - center_xy[0]
+                y_dist = ipp_xy[1] - center_xy[1]
+                x_index = (0.5*download_extent + scale_factor*x_dist).astype(np.int32)
+                y_index = (0.5*download_extent + scale_factor*y_dist).astype(np.int32)
+                #print('Getting sub area centered at x=%d y=%d' % (x_index,y_index))
+                #Reduce size of large array
+                if(x_index >= download_extent or y_index >= download_extent or x_index < 0 or y_index < 0):
+                    print('Error - Out of bounds!!')
+                    print('x = %d\ny = %d' % (x_index,y_index))
+                    return
+                sub_lf = lf[x_index-half_test:x_index+half_test,y_index-half_test:y_index+half_test,:]
+                sub_inac = inac[x_index-half_test:x_index+half_test,y_index-half_test:y_index+half_test,:]
+                sub_elev = elev[x_index-half_test:x_index+half_test,y_index-half_test:y_index+half_test]
+                map_data = lf_format(sub_inac,sub_lf,sub_elev)
+                plot_env(map_data,ipp_xy,find_xy,exp_name,env,incident)
+                #run LPM
+                #run_replicate(ipp, find, map, T, p_behavior, alpha, LL)
+            else:
+                map_data = lf_format(inac,lf,elev)
+                plot_env(map_data,ipp_xy,find_xy,exp_name,env,incident)
+                #run_replicate(ipp, find, map, T, p_behavior, alpha, LL)
 
         #Locate find position from IPP
         #Convert both to meters
@@ -104,22 +120,105 @@ def main(exp_name='test', start_idx=0, n_envs=1, n_iter=10):
         #But what I can do is walk back through the steps, adding to a 'heatmap' by doing +1 on cells that are visited
         #Then, we normalize it all to 0-1
 
-        run_replicate(ipp, find, map, T, p_behavior, alpha, LL)
+        #run_replicate(ipp, find, map, T, p_behavior, alpha, LL)
 
 
 
     return
 
+#Combines Linear and Inacessible features for all map layers
+#Returns tuple of map data for LPM input
+def lf_format(inac, lf, elev):
+
+    #Size error checking
+    extent = np.shape(elev)[0]
+    if(np.shape(lf)[0] != extent or np.shape(inac)[0] != extent):
+        print("Error - Terrain files are not the same size")
+        return
+    if(np.shape(lf)[0] != np.shape(lf)[1] or np.shape(inac)[0] != np.shape(inac)[1] or np.shape(elev)[0] != np.shape(elev)[1]):
+        print("Error - Collected terrain is not square")
+        return
+
+    #Put Linear features and inacessible layers into correct format. Takes 10-30s
+    print("Formatting terrain data...")
+    start_t = time.time()
+    bw_lf = np.zeros((extent,extent),dtype=np.int8())
+    bw_inac = np.zeros((extent,extent),dtype=np.int8())
+    inac_ct = 0
+    lf_ct = 0
+    for i in range(0,extent):
+        for j in range(0,extent):
+            lf_check = np.count_nonzero(lf[i][j])
+            if(lf_check > 0):
+                bw_lf[i][j] = 1
+                lf_ct += 1
+            inac_check = np.count_nonzero(inac[i][j])
+            if(inac_check > 0):
+                bw_inac[i][j] = 1
+                inac_ct +=1
+    end_t = time.time()
+    end_t = end_t - start_t
+    print("We have %d linear features and %d inacessible areas, it took %.1f seconds" % (lf_ct,inac_ct,end_t))
+    map_data = (bw_inac, bw_lf, elev)
+
+    return map_data
+
+#Creates a plot of the environment from saved data
+def plot_env(map_data, ipp, find, exp_name='test', env=0, incident=0):
+
+    scale_factor = 3/20
+    plot_dir = ''
+
+    lf = map_data[1]
+    inac = map_data[0]
+    elev = map_data[2]
+
+    extent = np.shape(lf)[0]
+
+    #Distance calculation
+    x_dist = find[0] - ipp[0]
+    y_dist = find[1] - ipp[1]
+    x_index = np.ceil(0.5*extent + scale_factor*x_dist).astype(np.int32)
+    y_index = np.ceil(0.5*extent + scale_factor*y_dist).astype(np.int32)
+
+    #For visualizing all terrain features
+    t = 'Area ' + str(env) + ' Incident ' + str(incident) + ' - Linear Features'
+    plt.title(t)
+    plt.imshow(lf)
+    #Show initial position and find point
+    plt.plot(1500,1500,'g*',label='IPP', ms=10)
+    plt.plot(x_index,y_index,'yX',label='Find Position', ms=10)
+    #plot_name = plot_dir + file_id + '_' + name_list[i]
+    #plt.savefig(plot_name)
+    plt.legend()
+    plt.show()
+    #plt.close()
+    t = 'Area ' + str(env) + 'Incident ' + str(incident) + ' - Inacessible Areas'
+    plt.title(t)
+    plt.imshow(inac)
+    plt.plot(1500,1500,'g*',label='IPP', ms=10)
+    plt.plot(x_index,y_index,'yX',label='Find Position', ms=10)
+    #plot_name = plot_dir + file_id + '_' + name_list[i]
+    #plt.savefig(plot_name)
+    plt.legend()
+    plt.show()
+
+    return
+
+#TODO: Delete this, not needed any more
+#Combines Linear and Inacessible features for all map layers
+def test_lf_format(exp_name='test',force_save='False'):
 
     #Set up which environments to run through.
-    start = 7
+    start = 8
     n_envs = 1
     if(n_envs == -1):
         n_envs = np.size(np.unique(areas))
 
     #Iterate through environment datasets
-    for i in range(start, start+n_envs):
-        dir = './map_layers/' + exp_name + '_' + str(i) + '/'
+    for env in range(start, start+n_envs):
+        dir = './map_layers/' + exp_name + '/' + exp_name + '_' + str(env) + '/'
+        print(dir)
         #Load datasets
         elev = np.load(dir+'elevation.npy')
         lf = np.load(dir+'linear_feats.npy')
@@ -136,30 +235,54 @@ def main(exp_name='test', start_idx=0, n_envs=1, n_iter=10):
 
         #Put Linear features and inacessible layers into correct format
         print("Formatting terrain data...")
+        start_t = time.time()
         bw_lf = np.zeros((extent,extent),dtype=np.int8())
         bw_inac = np.zeros((extent,extent),dtype=np.int8())
         inac_ct = 0
         lf_ct = 0
         for i in range(0,extent):
             for j in range(0,extent):
-                lf_check = np.sum(lf[i][j])
+                lf_check = np.count_nonzero(lf[i][j])
                 if(lf_check > 0):
                     bw_lf[i][j] = 1
                     lf_ct += 1
-                inac_check = np.sum(inac[i][j])
+                inac_check = np.count_nonzero(inac[i][j])
                 if(inac_check > 0):
                     bw_inac[i][j] = 1
                     inac_ct +=1
-        print("We have %d linear features and %d inacessible areas" % (lf_ct,inac_ct))
+        end_t = time.time()
+        end_t = end_t - start_t
+        print("We have %d linear features and %d inacessible areas, it took %.1f seconds" % (lf_ct,inac_ct,end_t))
         map_data = (bw_inac, bw_lf, elev)
 
+        #For visualizing all terrain features
+        t = exp_name + ' Area ' + str(env) + ' - Linear Features'
+        plt.title(t)
+        plt.imshow(bw_lf)
+        #Show initial position and find point
+        plt.plot(1500,1500,'g*',label='IPP', ms=10)
+        plt.plot(400,600,'yX',label='Find Position', ms=10)
+        #plot_name = plot_dir + file_id + '_' + name_list[i]
+        #plt.savefig(plot_name)
+        plt.legend()
+        plt.show()
+        #plt.close()
+        t = exp_name + ' Area ' + str(env) + ' - Inacessible Areas'
+        plt.title(t)
+        plt.imshow(bw_inac)
+        #plot_name = plot_dir + file_id + '_' + name_list[i]
+        #plt.savefig(plot_name)
+        plt.legend()
+        plt.show()
+
+
         #Pull incident locations for environment
-        i_area = np.argwhere(areas==i)
-        i_area = np.transpose(i_area)[0]
-        for incident in range(0,len(i_area)):
+        #i_area = np.argwhere(areas==i)
+        #i_area = np.transpose(i_area)[0]
+        #for incident in range(0,len(i_area)):
             #s = 'Environment {0} - start {1} find {2}'
             #s_lst = [i,]
-            print()
+        #    print()
             #For large downloads, re-size and center on IPP
             #if(len(i_area) > 1)):
 
@@ -180,8 +303,10 @@ def main(exp_name='test', start_idx=0, n_envs=1, n_iter=10):
 
 if __name__ == "__main__":
 
-    start = 8       #Environment to start at (0 is first, 8 is first large map)
-    n = 1           #Number of environments to test
+    start = 7       #Environment to start at (0 is first, 8 is first large map)
+    n = 2           #Number of environments to test
     mc = 10         #Monte Carlo iterations
 
-    main(exp_name='hiker',start_idx=start,n_envs=n,n_iter=mc)
+
+    #test_lf_format()
+    main(exp_name='trust',start_idx=start,n_envs=n,n_iter=mc)