Commit 743c17e8 authored by Bryson Howell's avatar Bryson Howell

Cleaned up terrain collection code

parent bb787a40
......@@ -82,7 +82,7 @@ def grab_features(anchor_point, extent, sample_dist = 10, case_name = 'blah', he
query_cnt = 0
while type(q)==list and query_cnt <= 30: # have to do this because arcgis is sketchy as hell and doesnt always come back
try:
print("querying {} layer...".format(name_list[i]))
print(" querying {} layer...".format(name_list[i]))
query_starttime = time.time()
#Want to query where the geom filter aligns
......@@ -99,13 +99,13 @@ def grab_features(anchor_point, extent, sample_dist = 10, case_name = 'blah', he
gis = GIS(url = 'https://virginiatech.maps.arcgis.com')
lyr = FeatureLayer(url=url, gis=gis)
print("query time {}".format(query_endtime - query_starttime))
#print("query time {}".format(query_endtime - query_starttime))
if query_cnt > 30 and not q:
print("{} layer failed too many times, leaving empty".format(name_list[i]))
continue
print("{} layer sucessfully queried".format(name_list[i]))
print("%s layer has %d points" %(name_list[i],len(q.features)))
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()
......@@ -200,7 +200,7 @@ def grab_features(anchor_point, extent, sample_dist = 10, case_name = 'blah', he
y_pts_inac = test_pts[mask,1]
inac_endtime = time.time()
print("{} inac took {}".format(j,inac_endtime - inac_starttime))
#print("{} inac took {}".format(j,inac_endtime - inac_starttime))
pts_inac = np.stack([x_pts_inac,y_pts_inac]).T
......@@ -275,7 +275,7 @@ 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 probably 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,
......
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
0,0,31.39708335,-111.20643,31.406333,-111.19683
1,1,31.42903334,-110.29328,31.429517,-110.30412
2,2,31.50571667,-110.67622,31.51,-110.65333
3,3,31.72261664,-110.18783,31.741767,-110.19178
4,4,31.83521665,-110.3567,31.862267,-110.39655
5,5,31.90240002,-109.27847,31.8395,-109.27667
6,6,31.92245,-109.96732,31.921483,-110.03358
7,7,32.01236667,-109.31572,31.99935,-109.3078
8,8,32.33333333,-110.85283,32.3715,-110.86467
9,8,32.33583333,-110.91017,32.3665,-110.88017
10,8,32.33699999,-110.91667,32.359833,-110.89767
11,9,32.40916665,-110.70983,32.400833,-110.69233
12,9,32.41976668,-110.74733,32.3655,-110.77583
13,9,32.43543332,-110.78933,32.428833,-110.79433
14,10,32.82141666,-111.20212,32.824667,-111.22717
15,11,33.37055,-111.11522,33.375283,-111.08918
16,12,33.38884999,-111.36565,33.4012,-111.37517
17,12,33.39680001,-111.34805,33.4215,-111.3605
18,12,33.39705,-111.34787,33.409333,-111.31717
19,12,33.39733334,-111.348,33.441667,-111.36767
20,12,33.39750001,-111.34783,33.415667,-111.365
21,13,33.47801666,-111.43767,33.415667,-111.36478
22,14,33.60398331,-112.51512,33.590333,-112.52283
23,15,33.70541668,-111.33805,33.67915,-111.31695
24,16,34.0927,-111.42465,34.068767,-111.43592
25,16,34.09278333,-111.42098,34.102333,-111.49033
26,17,34.55,-111.63333,34.55,-111.61667
27,17,34.55976664,-111.65395,34.540067,-111.70348
28,18,34.6,-112.55,34.633333,-112.63333
29,18,34.63041668,-112.5553,34.627233,-112.54577
30,19,34.82166665,-111.80667,34.818333,-111.79833
31,19,34.82386665,-111.77513,34.831667,-111.743
32,19,34.86666667,-111.88333,34.866667,-111.8
33,19,34.88683332,-111.784,34.902667,-111.78683
34,19,34.89333331,-111.86333,34.906667,-111.87603
35,20,34.90286668,-111.81313,34.891833,-111.80482
36,20,34.91666667,-111.8,34.966667,-111.86667
37,20,34.92530003,-111.7341,34.897333,-111.74117
38,21,34.97868334,-111.89643,34.946833,-111.88967
39,22,35.1938,-114.05703,35.208717,-114.13195
40,23,35.24375,-111.59967,35.232517,-111.6015
41,24,35.33068333,-111.71108,35.334233,-111.6975
42,25,36.23878333,-112.6892,36.246433,-112.70008
43,26,41.54819,-80.33056,41.59006,-80.32647
44,27,42.0097,-74.42595,42.03544,-74.35565
45,27,42.02893,-74.33659,42.02902,-74.35191
46,28,42.17965,-74.21362,42.18345,-74.19585
47,29,42.31735,-76.47791,42.30271,-76.48976
48,30,42.34432,-77.47638,42.36077,-77.48593
49,31,42.74271,-73.45475,42.75163,-73.46259
50,32,43.06902,-74.48481,43.05475,-74.4839
51,33,43.42473,-73.73209,43.42878,-73.73981
52,34,43.42498,-74.41496,43.41592,-74.4142
53,34,43.42736,-74.4481,43.44347,-74.45004
54,34,43.4332,-74.41433,43.45684,-74.41519
55,34,43.4449,-74.4086,43.43407,-74.40747
56,35,43.51063,-74.57393,43.53014,-74.57169
57,36,43.65649,-76.00019,43.65566,-76.0115
58,37,43.73413,-74.25577,43.73701,-74.28305
59,38,43.8756,-74.43076,43.91132,-74.37437
60,39,43.95385,-75.15748,43.99774,-75.16148
61,40,44.16065,-73.85545,44.1547,-73.85942
62,41,44.19013,-74.81336,44.18772,-74.79719
63,42,44.28656,-74.61429,44.28923,-74.60208
64,43,48.1103,-121.4917,48.142,-121.4739
import numpy as np
import matplotlib
import matplotlib.pyplot as plt
import os
import os, sys, inspect
from os import path, getcwd
import pandas as pd
#Runs replicates of lost person model on data sets collected by lpm_maps.py
currentdir = os.path.dirname(os.path.abspath(inspect.getfile(inspect.currentframe())))
parentdir = os.path.dirname(currentdir)
sys.path.insert(0, parentdir)
from arcgis_terrain import get_terrain_map, lat_lon2meters
from LostPersonModel.main_hiker import run_replicate #change this import to where you have the LPM repo
#Runs replicates of lost person model on data sets collected by lpm_maps.py, producing a heatmap
#Bryson Howell, 7/24/24
def main(exp_name='test', n_envs=1, n_iter=10):
test_extent = 20000
download_extent = 40000
#Pull incident locations
listname = './lp_data/indexed_incidents.csv'
incidents = pd.read_csv(listname)
ipp_latlons = incidents[["IPP_lat","IPP_lon"]].to_numpy(dtype=np.float32())
find_latlons = incidents[["find_lat","find_lon"]].to_numpy(dtype=np.float32())
areas = incidents['area'].to_numpy(dtype=np.float32())
start = 0
n_envs = np.size(np.unique(areas))
for i in range(start, start+n_envs):
#Pull incident locations for environment
i_area = np.argwhere(areas==i)
i_area = np.transpose(i_area)[0]
s = 'Area {0} - indices are {1}'
lst = [i,i_area]
print(s.format(*lst))
#test - print start / end locations on grid
grid = np.zeros((3000,3000),dtype=np.int8())
for incident in range(0,len(i_area)):
#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)
#Convert
#Locate find position from IPP
#Convert both to meters
#Use extent / cell count to determine which cell we're in:
return
#Set up which environments to run through.
start = 7
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) + '/'
#Load datasets
......@@ -49,11 +111,21 @@ def main(exp_name='test', n_envs=1, n_iter=10):
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())
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()
#For large downloads, re-size and center on IPP
#if(len(i_area) > 1)):
#for j in range
#we know initial point is at 0,0 but where is the final point...?
return
......
......@@ -169,31 +169,17 @@ def investigate():
print(np.shape(elev)) #6000 x 6000. Got this using res of 25, not 10...
print(elev[0][0])
#Make sure we can translate start / end grid cells...
return
#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(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 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
#Looks through a list of IPP and find locations for SAR incidents
#Groups IPPs that can fit within a map of extent_download size, such that a sub-map of extent_test size is centered on each IPP
def group_IPP(listname, extent_download=40000, extent_test=20000):
#Parameters for feature_set
res = 25
folder = folder_name
#Load in IPP and Find locations from csv file
incidents = pd.read_csv(listname)
latlons = incidents[["IPP_lat","IPP_lon"]].to_numpy(dtype=np.float32())
ipp_x, ipp_y = lat_lon2meters(latlons[:,0],latlons[:,1])
#Collection of map centers, in meters. Use to check if we need to download a new area
ipp_list = []
......@@ -202,48 +188,6 @@ def collect_terrain(folder_name='test',start_idx=0,max_collect=1):
#Used to group search incidents in the same area together. Added to DF at end
keys = []
listname = './lp_data/incident_locations.csv'
#samples = 1
#incidents = pd.read_csv(listname,nrows=samples)
incidents = pd.read_csv(listname)
samples = incidents.shape[0]
try_kmeans = 0
if(try_kmeans == 1):
latlons = incidents[["IPP_lat","IPP_lon"]].to_numpy(dtype=np.float32())
#latlons = np.concatenate(([incidents['IPP_lat']],[incidents['IPP_lon']]),axis=0)
print(latlons[0])
#latlons = whiten(latlons)
#clusters, distortions = kmeans(whiten(latlons),4)
clusters, distortions = kmeans(latlons,n_clusters)
plt.scatter(latlons[:,0],latlons[:,1])
plt.scatter(clusters[:,0],clusters[:,1],c='r')
plt.title('Lat/Lon Coordinates')
plt.show()
#Convert coordinates to x/y
clusters_x, clusters_y = lat_lon2meters(clusters[:,0],clusters[:,1])
ipp_x, ipp_y = lat_lon2meters(latlons[:,0],latlons[:,1])
print(clusters_x)
fig = plt.figure()
ax = fig.add_subplot(111)
plt.scatter(ipp_x,ipp_y)
plt.scatter(clusters_x,clusters_y,c='r')
#try drawing boxes (in meters) around the areas we want to download
for i in range(0, len(clusters_x)):
rect = matplotlib.patches.Rectangle((clusters_x[i]-0.5*extent_download,clusters_y[i]-0.5*extent_download),extent_download,extent_download,facecolor="#86cecb",alpha=0.5)
ax.add_patch(rect)
plt.title('Meters Coordinates')
plt.show()
return
latlons = incidents[["IPP_lat","IPP_lon"]].to_numpy(dtype=np.float32())
ipp_x, ipp_y = lat_lon2meters(latlons[:,0],latlons[:,1])
#This gives a loose clustering algorithm, such that we know maps can be made from groups
#Loop through search incidents
for index, row in incidents.iterrows():
......@@ -278,18 +222,17 @@ def collect_terrain(folder_name='test',start_idx=0,max_collect=1):
else:
keys.append(closest)
print("Found %d maps to collect." % len(ipp_list))
if(max_collect < 0):
max_down = len(ipp_list)
print(keys)
#Add group keys to dataframe
incidents.insert(1,'area',keys)
#Fix matlab indexing
incidents['incident_index'] = incidents['incident_index'] - 1
#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.
#For visualizing collected data on a map. Todo: show world map as background (https://stackoverflow.com/questions/53233228/plot-latitude-longitude-from-csv-in-python-3-6)
draw = False
fig = plt.figure()
ax = fig.add_subplot(111)
......@@ -307,20 +250,104 @@ def collect_terrain(folder_name='test',start_idx=0,max_collect=1):
plt.title('Experiment Areas')
plt.show()
#Attempt to do kmeans to group areas by location.
#Problem is, this does not guarantee there is an area of test_extent size around each IPP
try_kmeans = 0
if(try_kmeans == 1):
latlons = incidents[["IPP_lat","IPP_lon"]].to_numpy(dtype=np.float32())
#latlons = np.concatenate(([incidents['IPP_lat']],[incidents['IPP_lon']]),axis=0)
print(latlons[0])
#latlons = whiten(latlons)
#clusters, distortions = kmeans(whiten(latlons),4)
clusters, distortions = kmeans(latlons,n_clusters)
plt.scatter(latlons[:,0],latlons[:,1])
plt.scatter(clusters[:,0],clusters[:,1],c='r')
plt.title('Lat/Lon Coordinates')
plt.show()
#Convert coordinates to x/y
clusters_x, clusters_y = lat_lon2meters(clusters[:,0],clusters[:,1])
ipp_x, ipp_y = lat_lon2meters(latlons[:,0],latlons[:,1])
print(clusters_x)
fig = plt.figure()
ax = fig.add_subplot(111)
plt.scatter(ipp_x,ipp_y)
plt.scatter(clusters_x,clusters_y,c='r')
#try drawing boxes (in meters) around the areas we want to download
for i in range(0, len(clusters_x)):
rect = matplotlib.patches.Rectangle((clusters_x[i]-0.5*extent_download,clusters_y[i]-0.5*extent_download),extent_download,extent_download,facecolor="#86cecb",alpha=0.5)
ax.add_patch(rect)
plt.title('Meters Coordinates')
plt.show()
return
return
#From the set of LPM coordinates, download and save map/linear features for areas as necessary
def collect_terrain(folder_name='test',start_idx=0,max_collect=1,format_ipps=False):
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 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 = folder_name
#Load in IPP and Find locations from csv file
listname = './lp_data/indexed_incidents.csv'
incidents = pd.read_csv(listname)
ipp_latlons = incidents[["IPP_lat","IPP_lon"]].to_numpy(dtype=np.float32())
find_latlons = incidents[["find_lat","find_lon"]].to_numpy(dtype=np.float32())
areas = incidents['area'].to_numpy(dtype=np.float32())
#Limit what we download
if(max_down < 0):
max_down = np.size(np.unique(areas))
end_idx = start_idx + max_down
if(end_idx > np.size(np.unique(areas))):
end_idx = np.size(np.unique(areas))
#Decide if we want to save pics, and get the elev map (takes a while)
do_plot = True
do_elev = False
#Loop through unique areas to download
for i in range(start_idx, end_idx):
i_area = np.argwhere(areas==i)
i_area = np.transpose(i_area)[0]
#s = 'Area {0} - has points {1}'
#lst = [i,i_area]
#rint(s.format(*lst))
#Now download area centered on first index of area
print("Collecting GIS data %d/%d" % ((i-start_idx)+1,(end_idx-start_idx)))
#Download local map
if(len(i_area) > 0):
size = extent_test
else:
size = extent_download
ipp = ipp_latlons[i_area[0]]
s = ' Getting Area {0} with IPP(s) {1} - center is {2}'
lst = [i,i_area,ipp]
print(s.format(*lst))
folder_i = folder + '/' + folder + "_" + str(i)
grab_features(ipp, size, sample_dist = res,
case_name = folder_i, save_files = False, save_to_folder = True,
file_id = str(i), plot_data=do_plot, get_elev=do_elev)
print(" Saved terrain in %s" % folder_i)
return
#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_[]'
#Limit what we download
if(max_down < 0):
max_down = len(ipp_list)
#iterate through files in directory
#for file in os.listdir(directory):
#see if latitude matches
# for ilat in range(0,latlons.shape(0)):
# print(latlons[i][0])
#See if longitude matches
latlons = incidents[["IPP_lat","IPP_lon"]].to_numpy(dtype=np.float32())
#Now, if there's only one key in an area we can avoid downloading a lot.
......@@ -339,8 +366,8 @@ def collect_terrain(folder_name='test',start_idx=0,max_collect=1):
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)
folder_i = folder + "_" + str(i)
grab_features(ipp_lat_list[i], size, sample_dist = res, case_name = folder_i, save_files = False, save_to_folder = True, file_id = str(i), plot_data=do_plot, get_elev=do_elev)
return
......@@ -360,23 +387,19 @@ def collect_terrain(folder_name='test',start_idx=0,max_collect=1):
#[e,e_interp,x,y,data,lat_lon] = get_terrain_map(ipp_point, sample_dist = res, extent = size, show_plot = True)
#Load Linear Features. Probably need to convert some Matlab code...
#return
#Collect GIS maps and layers for SAR initial positions
def main():
listname = './lp_data/incident_locations.csv' #Collection of search incidents
extent = 3000 #size in meters of map area
#incident_list()
#gis_test()
collect_terrain(folder_name='big',start_idx=7,max_collect=1)
#group_IPP(listname) #Format incident data correctly
#Friday, 7/26 - make sure we get window size of 625 with larger environment.
#Also, try reducing size of downloads...
collect_terrain(folder_name='test',start_idx=0,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