Skip to content

A
ags_grabber
Commit bb787a40 authored by Bryson Howell's avatar Bryson Howell
Browse files
Options

Starting on running LPM replicates

parent 32b71c53
Hide whitespace changes
Inline Side-by-side
Showing with 195 additions and 62 deletions
arcgis_terrain.py
View file @ bb787a40
......@@ -139,6 +139,7 @@ def get_terrain_map(lat_lon = [0,0], sample_dist = 10, extent = 100, heading = 0
if verbosity:
print("total calls: {}".format(np.square(cc)))
print("For extent of %d, we're making %d calls with window %d" % (extent,cc,sample_extents))
# set max values
elv_layer.extent['xmin'] = lhc_pt[0]
......@@ -223,7 +224,7 @@ def get_terrain_map(lat_lon = [0,0], sample_dist = 10, extent = 100, heading = 0
factor = scaled_extent/e.shape[0]
e_interp = ndimage.zoom(e, factor, order = 3)
print(show_plot)
#print(show_plot)
if show_plot:
# Attaching 3D axis to the figure
grid_x, grid_y = np.mgrid[np.min(x):np.max(x):100j, np.min(y):np.max(y):100j]
......@@ -235,6 +236,9 @@ def get_terrain_map(lat_lon = [0,0], sample_dist = 10, extent = 100, heading = 0
fig = plt.figure()
ax = fig.gca(projection='3d')
ax.plot_surface(grid_x, grid_y, grid_z, cmap='viridis')
plt.savefig('./test_3d_plot')
plt.show()
plt.close()
if verbosity:
print("x min: {}".format(np.min(x)))
......
feature_set.py
View file @ bb787a40
......@@ -88,14 +88,6 @@ def grab_features(anchor_point, extent, sample_dist = 10, case_name = 'blah', he
#Want to query where the geom filter aligns
q = lyr.query(return_count_only=False, return_ids_only=False, return_geometry=True,
out_sr='3857', geometry_filter=geom_filter)
#Catch ArcGIS authentication errors (problem was actually with input coords)
#except:
# print("Could not load layer")
# return
#lyr2 = FeatureLayer(url = 'https://carto.nationalmap.gov/arcgis/rest/services/transportation/MapServer/30', gis = gis)
#q2 = lyr2.query(return_count_only=False, return_ids_only=False, return_geometry=True,
#out_sr='3857', geometry_filter=geom_filter)
query_endtime = time.time()
except (json.decoder.JSONDecodeError, TypeError) as e:
......@@ -111,25 +103,14 @@ def grab_features(anchor_point, extent, sample_dist = 10, case_name = 'blah', he
if query_cnt > 30 and not q:
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"
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"
np.savetxt(fn, bin_map, delimiter=",", fmt='%f')
else:
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"
np.savetxt(fn, bin_map, delimiter=",", fmt='%f')
continue
print("{} layer sucessfully queried".format(name_list[i]))
print("%s layer has %d points" %(name_list[i],len(q.features)))
# re-build into list of x-y values
# feat_points = []
query_dict = q.to_dict()
#print(query_dict.keys()) #features, displayFieldName, spatialReference, geometryType, fields
#We look throug the features.
for j,feat in enumerate(query_dict['features']):
......@@ -255,6 +236,7 @@ def grab_features(anchor_point, extent, sample_dist = 10, case_name = 'blah', he
# print("done with feature {}".format(j))
#maybe try making these sparse arrays, as they're mostly empty... (scipy.sparse.bsr_array is promising)
# add inac feats to viz map
if name_list[i] in inac_layers:
viz_map_inac[:, :, viz_cnt_inac] = inac_bin_map
......@@ -280,6 +262,7 @@ def grab_features(anchor_point, extent, sample_dist = 10, case_name = 'blah', he
fn = "./map_layers/" + name_list[i]+"_data_"+file_id+".csv"
np.savetxt(fn,bin_map,delimiter=",", fmt='%f')
os.makedirs('./map_layers/'+case_name+'/', exist_ok=True)
# save linear features as numpy file
fn = './map_layers/' + case_name + '/linear_feats'
np.save(fn,viz_map)
......@@ -292,11 +275,13 @@ def grab_features(anchor_point, extent, sample_dist = 10, case_name = 'blah', he
# save terrain as csv file (this method is pretty slow, but can compensate with interp)
if(get_elev):
print("Getting terrain map, this will take a while....")
print("Getting terrain map, this will probably take a while....")
[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=False) # because flipping
fn = './map_layers/' + case_name + '/elevation'
np.save(fn,e_interp)
if save_files:
if save_to_folder:
elv_filename = "./map_layers/" + case_name + "/elv_data_" + file_id + ".csv"
......
lp_data/indexed_incidents.csv 0 → 100644
View file @ bb787a40
incident_index,area,IPP_lat,IPP_lon,find_lat,find_lon
1,0,31.39708335,-111.20643,31.406333,-111.19683
2,1,31.42903334,-110.29328,31.429517,-110.30412
3,2,31.50571667,-110.67622,31.51,-110.65333
4,3,31.72261664,-110.18783,31.741767,-110.19178
5,4,31.83521665,-110.3567,31.862267,-110.39655
6,5,31.90240002,-109.27847,31.8395,-109.27667
7,6,31.92245,-109.96732,31.921483,-110.03358
8,7,32.01236667,-109.31572,31.99935,-109.3078
9,8,32.33333333,-110.85283,32.3715,-110.86467
10,8,32.33583333,-110.91017,32.3665,-110.88017
11,8,32.33699999,-110.91667,32.359833,-110.89767
12,9,32.40916665,-110.70983,32.400833,-110.69233
13,9,32.41976668,-110.74733,32.3655,-110.77583
14,9,32.43543332,-110.78933,32.428833,-110.79433
15,10,32.82141666,-111.20212,32.824667,-111.22717
16,11,33.37055,-111.11522,33.375283,-111.08918
17,12,33.38884999,-111.36565,33.4012,-111.37517
18,12,33.39680001,-111.34805,33.4215,-111.3605
19,12,33.39705,-111.34787,33.409333,-111.31717
20,12,33.39733334,-111.348,33.441667,-111.36767
21,12,33.39750001,-111.34783,33.415667,-111.365
22,13,33.47801666,-111.43767,33.415667,-111.36478
23,14,33.60398331,-112.51512,33.590333,-112.52283
24,15,33.70541668,-111.33805,33.67915,-111.31695
25,16,34.0927,-111.42465,34.068767,-111.43592
26,16,34.09278333,-111.42098,34.102333,-111.49033
27,17,34.55,-111.63333,34.55,-111.61667
28,17,34.55976664,-111.65395,34.540067,-111.70348
29,18,34.6,-112.55,34.633333,-112.63333
30,18,34.63041668,-112.5553,34.627233,-112.54577
31,19,34.82166665,-111.80667,34.818333,-111.79833
32,19,34.82386665,-111.77513,34.831667,-111.743
33,19,34.86666667,-111.88333,34.866667,-111.8
34,19,34.88683332,-111.784,34.902667,-111.78683
35,19,34.89333331,-111.86333,34.906667,-111.87603
36,20,34.90286668,-111.81313,34.891833,-111.80482
37,20,34.91666667,-111.8,34.966667,-111.86667
38,20,34.92530003,-111.7341,34.897333,-111.74117
39,21,34.97868334,-111.89643,34.946833,-111.88967
40,22,35.1938,-114.05703,35.208717,-114.13195
41,23,35.24375,-111.59967,35.232517,-111.6015
42,24,35.33068333,-111.71108,35.334233,-111.6975
43,25,36.23878333,-112.6892,36.246433,-112.70008
44,26,41.54819,-80.33056,41.59006,-80.32647
45,27,42.0097,-74.42595,42.03544,-74.35565
46,27,42.02893,-74.33659,42.02902,-74.35191
47,28,42.17965,-74.21362,42.18345,-74.19585
48,29,42.31735,-76.47791,42.30271,-76.48976
49,30,42.34432,-77.47638,42.36077,-77.48593
50,31,42.74271,-73.45475,42.75163,-73.46259
51,32,43.06902,-74.48481,43.05475,-74.4839
52,33,43.42473,-73.73209,43.42878,-73.73981
53,34,43.42498,-74.41496,43.41592,-74.4142
54,34,43.42736,-74.4481,43.44347,-74.45004
55,34,43.4332,-74.41433,43.45684,-74.41519
56,34,43.4449,-74.4086,43.43407,-74.40747
57,35,43.51063,-74.57393,43.53014,-74.57169
58,36,43.65649,-76.00019,43.65566,-76.0115
59,37,43.73413,-74.25577,43.73701,-74.28305
60,38,43.8756,-74.43076,43.91132,-74.37437
61,39,43.95385,-75.15748,43.99774,-75.16148
62,40,44.16065,-73.85545,44.1547,-73.85942
63,41,44.19013,-74.81336,44.18772,-74.79719
64,42,44.28656,-74.61429,44.28923,-74.60208
65,43,48.1103,-121.4917,48.142,-121.4739
lpm_hikersim.py 0 → 100644
View file @ bb787a40
import numpy as np
import matplotlib
import matplotlib.pyplot as plt
import os
from os import path, getcwd
import pandas as pd
#Runs replicates of lost person model on data sets collected by lpm_maps.py
#Bryson Howell, 7/24/24
def main(exp_name='test', n_envs=1, n_iter=10):
start = 7
for i in range(start, start+n_envs):
dir = './map_layers/' + exp_name + '_' + str(i) + '/'
#Load datasets
elev = np.load(dir+'elevation.npy')
lf = np.load(dir+'linear_feats.npy')
inac = np.load(dir+'inac_layers.npy')
#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
print("Formatting terrain data...")
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])
if(lf_check > 0):
bw_lf[i][j] = 1
lf_ct += 1
inac_check = np.sum(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))
map_data = (bw_inac, bw_lf, elev)
#Pull incident locations for environment
listname = './lp_data/indexed_incidents.csv'
incidents = pd.read_csv(listname)
latlons = incidents[["IPP_lat","IPP_lon"]].to_numpy(dtype=np.float32())
areas = incidents['area'].to_numpy(dtype=np.float32())
return
if __name__ == "__main__":
n_env = 1 #Number of environments to test
n_iter = 10
main()
lpm_maps.py
View file @ bb787a40
......@@ -182,16 +182,18 @@ def investigate():
#From the set of LPM coordinates, download and save map/linear features for areas as necessary
#Also make images to show where areas are.
def collect_terrain():
def collect_terrain(folder_name='test',start_idx=0,max_collect=1):
extent_download = 40000 #Size of maps to dowload, in meters (big...)
extent_test = 20000 #Size of experiment area. 20km = 3000 grid cells
n_clusters = 10 #Number of areas to download
count_thres = 3 #Minimum number of points needed for a large download
n_clusters = 10 #Number of areas to download for kmeans
count_thres = 3 #Minimum number of points needed for a large download (not used any more, hard to implement.)
max_down = max_collect #Maximum number of things to download. Set to -1 to download all environments
start_down = start_idx #Index to start downloading from
#Parameters for feature_set
res = 25
folder = 'init'
folder = folder_name
#Collection of map centers, in meters. Use to check if we need to download a new area
ipp_list = []
......@@ -276,18 +278,41 @@ def collect_terrain():
else:
keys.append(closest)
print("Found %d maps to collect." % len(ipp_list))
if(max_collect < 0):
max_down = len(ipp_list)
print(keys)
#print(ipp_list)
#print(ipp_list[0][0])
#Add group keys to dataframe
incidents.insert(0,'area',keys)
incidents.insert(1,'area',keys)
#Save new list as csv
new_incidents = True
if(new_incidents):
incidents.to_csv('./lp_data/indexed_incidents.csv',index=False)
return
#For visualizing collected data on a map.
draw = False
fig = plt.figure()
ax = fig.add_subplot(111)
if(draw):
#Plot points
plt.scatter(ipp_x,ipp_y)
#Draw region areas
for i in range(0,len(ipp_list)):
#Draw download area
if(keys.count(i) > 1):
rect = matplotlib.patches.Rectangle((ipp_list[i][0]-0.5*extent_test,ipp_list[i][1]-0.5*extent_test),extent_test,extent_test,facecolor="#ed3419",alpha=0.5)
else:
rect = matplotlib.patches.Rectangle((ipp_list[i][0]-0.5*extent_download,ipp_list[i][1]-0.5*extent_download),extent_download,extent_download,facecolor="#86cecb",alpha=0.5)
ax.add_patch(rect)
plt.title('Experiment Areas')
plt.show()
#GIS is not working. Instead, look through matlab_data for matching IPP
matdir = '../ags_grabber/matlab_data/'
directory = os.fsencode(matdir)
fname = 'BW_LFandInac-Zelev_[]'
#matdir = '../ags_grabber/matlab_data/'
#directory = os.fsencode(matdir)
#fname = 'BW_LFandInac-Zelev_[]'
#iterate through files in directory
#for file in os.listdir(directory):
......@@ -300,41 +325,23 @@ def collect_terrain():
#Now, if there's only one key in an area we can avoid downloading a lot.
#Need to ouput data in a way that indicates if an IPP is a big area or not...
max_down = 3 #Maximum number of things to download
count = 0
do_plot = True
do_elev = True
for i in range(0,len(ipp_lat_list)):
for i in range(start_down,len(ipp_lat_list)):
print("Collecting GIS data %d/%d" % (i,len(ipp_lat_list)))
#Download local map
if(keys.count(i) <= count_thres):
if(keys.count(i) < count_thres):
size = extent_test
if(count < max_down):
count = count + 1
print('Getting initial point - ')
print(ipp_lat_list[i])
folder = "init_" + str(i)
grab_features(ipp_lat_list[i], size, sample_dist = res, case_name = folder, save_files = True, save_to_folder = True, file_id = str(i), plot_data=do_plot, get_elev=do_elev)
else:
size = extent_download
if(count < max_down):
count = count + 1
print('Getting initial point - ')
print(ipp_lat_list[i])
folder = folder + "_" + str(i)
grab_features(ipp_lat_list[i], size, sample_dist = res, case_name = folder, save_files = False, save_to_folder = True, file_id = str(i), plot_data=do_plot, get_elev=do_elev)
#For visualizing
draw = True
fig = plt.figure()
ax = fig.add_subplot(111)
if(draw):
#Plot points
plt.scatter(ipp_x,ipp_y)
#Draw region areas
for i in range(0,len(ipp_list)):
#Draw download area
if(keys.count(i) <= count_thres):
rect = matplotlib.patches.Rectangle((ipp_list[i][0]-0.5*extent_test,ipp_list[i][1]-0.5*extent_test),extent_test,extent_test,facecolor="#ed3419",alpha=0.5)
else:
rect = matplotlib.patches.Rectangle((ipp_list[i][0]-0.5*extent_download,ipp_list[i][1]-0.5*extent_download),extent_download,extent_download,facecolor="#86cecb",alpha=0.5)
ax.add_patch(rect)
plt.title('Experiment Areas')
plt.show()
return
......@@ -369,7 +376,7 @@ def main():
#incident_list()
#gis_test()
collect_terrain()
collect_terrain(folder_name='big',start_idx=7,max_collect=1)
#investigate()
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment