Can now download linear features, inac, and elevation data in desired format

arcgis_terrain.py

@@ -105,7 +105,8 @@ 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(api_key="AAPKe7cb1ab4f2ba44748cf53ac4f30d0caavu5t_-uMdWP8SRjtoea3s66-cRyFZmMTNx4rqy2w5sjwegU_tbJyBesd0LCMmdtV")
+    #gis = GIS(url="http://virginiatech.maps.arcgis.com")
+    gis = GIS(api_key="AAPKe7cb1ab4f2ba44748cf53ac4f30d0caavu5t_-uMdWP8SRjtoea3s66-cRyFZmMTNx4rqy2w5sjwegU_tbJyBesd0LCMmdtV") # linked Bryson's pro account
     if verbosity:
         #print("Successfully logged in as: " + gis.properties.user.username)
         print("Logged into GIS")
@@ -152,6 +153,7 @@ def get_terrain_map(lat_lon = [0,0], sample_dist = 10, extent = 100, heading = 0
     e = np.empty([0,sc*cc])
     data = []
     for j in range(cc): # outter loop for y values
+        print('Iteration %d/%d' % (j,cc))
         # lists of values for a single row, reset at every y iteration
         x_row = np.empty([sc,0])
         y_row = np.empty([sc,0])
@@ -182,7 +184,7 @@ def get_terrain_map(lat_lon = [0,0], sample_dist = 10, extent = 100, heading = 0
                 qs = np.stack([xs.reshape(len(elv_samples)), ys.reshape(len(elv_samples))]).T
                 es = es.reshape(len(elv_samples))
                 # data came back in weird ordering, need to re-order
-                print("re-ordering data ...")
+                #print("re-ordering data ...")
                 es_square = np.zeros([sc,sc])
                 for scx in range(sc):
                     for scy in range(sc): # maybe the least efficient way to do this
@@ -193,7 +195,7 @@ def get_terrain_map(lat_lon = [0,0], sample_dist = 10, extent = 100, heading = 0
                         except ValueError as e:
                             print("hellfire")
                 es = es_square
-                print("done re-ordering")
+                #print("done re-ordering")
 
             # then just the tuple of all
             data_temp = []

feature_set.py

@@ -29,7 +29,7 @@ def trim_extent(x_pts, y_pts, scaled_extent):
     y_pts_trimmed = y_pts[np.invert(rm_mask)] # trim points
     return [x_pts_trimmed, y_pts_trimmed]
 
-def grab_features(anchor_point, extent, sample_dist = 10, case_name = 'blah', heading = 0, save_files = False, save_to_folder = False, file_id = 'temp', plot_data = False):
+def grab_features(anchor_point, extent, sample_dist = 10, case_name = 'blah', heading = 0, save_files = False, save_to_folder = False, file_id = 'temp', plot_data = False, get_elev = False):
     roads_url = "https://carto.nationalmap.gov/arcgis/rest/services/transportation/MapServer/30"
     river_url = "https://hydro.nationalmap.gov/arcgis/rest/services/nhd/MapServer/6"
     riverw_url = "https://hydro.nationalmap.gov/arcgis/rest/services/nhd/MapServer/8"
@@ -46,15 +46,18 @@ def grab_features(anchor_point, extent, sample_dist = 10, case_name = 'blah', he
 
     ky = "AAPKe7cb1ab4f2ba44748cf53ac4f30d0caavu5t_-uMdWP8SRjtoea3s66-cRyFZmMTNx4rqy2w5sjwegU_tbJyBesd0LCMmdtV"
     tkn = "3NKHt6i2urmWtqOuugvr9ZEfiuh0-3Amg_JW92hWw6MrRjQKzJfjBbFhwGtOqYD3IEnYsljGDxUWTXa4b2HGR9WawiS9jXoVzqX7qF3-y5y3jPSTTL5TYKO8tHiitpqK"
-    gis = GIS(token=tkn,api_key=ky)
+    #gis = GIS(token=tkn,api_key=ky)
     #gis = GIS(api_key="AAPKe7cb1ab4f2ba44748cf53ac4f30d0caavu5t_-uMdWP8SRjtoea3s66-cRyFZmMTNx4rqy2w5sjwegU_tbJyBesd0LCMmdtV") # linked Bryson's pro account
+    gis = GIS(url = 'https://virginiatech.maps.arcgis.com') #New API does not work for api keys?
     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.int32)
+    #print("Making grid of size %d" % scaled_extent)
     viz_cnt = 0
-    viz_map = np.zeros([scaled_extent,scaled_extent,len(name_list)+len(inac_layers)])
+    viz_cnt_inac = 0
+    viz_map = np.zeros([scaled_extent,scaled_extent,len(name_list)])
+    viz_map_inac = np.zeros([scaled_extent,scaled_extent,len(inac_layers)])
 
     for i,url in enumerate(url_list):
 
@@ -62,6 +65,8 @@ def grab_features(anchor_point, extent, sample_dist = 10, case_name = 'blah', he
         bin_map = np.zeros([scaled_extent,scaled_extent])
         inac_bin_map = np.zeros([scaled_extent,scaled_extent])
 
+        #Array of rings and a spatial reference
+        #Docs - https://developers.arcgis.com/python/api-reference/arcgis.geometry.html
         geom = arcgis.geometry.Polygon({'spatialReference': {"wkid" : 3857},
                              'rings': [[
                                 [ap_meters[0] - (extent/2), ap_meters[1] + (extent/2)],
@@ -78,12 +83,19 @@ def grab_features(anchor_point, extent, sample_dist = 10, case_name = 'blah', he
         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(lyr)
                 query_starttime = time.time()
 
+                #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)
-
+                             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:
@@ -92,7 +104,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(api_key="AAPKe7cb1ab4f2ba44748cf53ac4f30d0caavu5t_-uMdWP8SRjtoea3s66-cRyFZmMTNx4rqy2w5sjwegU_tbJyBesd0LCMmdtV") # linked my arcgis pro account
+                gis = GIS(url = 'https://virginiatech.maps.arcgis.com')
                 lyr = FeatureLayer(url=url, gis=gis)
 
         print("query time {}".format(query_endtime - query_starttime))
@@ -114,10 +126,11 @@ def grab_features(anchor_point, extent, sample_dist = 10, case_name = 'blah', he
                         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()
+
         for j,feat in enumerate(query_dict['features']):
 
             # print("starting feature {} ...".format(j))
@@ -177,7 +190,7 @@ def grab_features(anchor_point, extent, sample_dist = 10, case_name = 'blah', he
 
             interp_endtime = time.time()
 
-            print("{} interpolation took {}".format(j,interp_endtime - interp_starttime))
+            #print("{} interpolation took {}".format(j,interp_endtime - interp_starttime))
 
             x_pts = pts_interp[0]
             y_pts = pts_interp[1]
@@ -242,47 +255,68 @@ def grab_features(anchor_point, extent, sample_dist = 10, case_name = 'blah', he
 
             # print("done with feature {}".format(j))
 
-        # add to viz map
+        # add inac feats to viz map
         if name_list[i] in inac_layers:
-            viz_map[:, :, viz_cnt] = inac_bin_map
-            viz_cnt = viz_cnt + 1
+            viz_map_inac[:, :, viz_cnt_inac] = inac_bin_map
+            viz_cnt_inac = viz_cnt_inac + 1
             if save_files:
                 if save_to_folder:
                     fn = "./map_layers/" + case_name + "/" + name_list[i] + "_inac_data_" + file_id + ".csv"
+                    os.makedirs('./map_layers/'+case_name+'/', exist_ok=True)
                     np.savetxt(fn, inac_bin_map, delimiter=",", fmt='%f')
                 else:
                     fn = "./map_layers/" + name_list[i] + "_inac_data_" + file_id + ".csv"
                     np.savetxt(fn, inac_bin_map, delimiter=",", fmt='%f')
 
-
+        #Add feats to viz map
         viz_map[:, :, viz_cnt] = bin_map
         viz_cnt = viz_cnt + 1
         if save_files:
             if save_to_folder:
                 fn = "./map_layers/" + case_name + "/" + name_list[i]+"_data_"+file_id+".csv"
+                os.makedirs('./map_layers/'+case_name+'/', exist_ok=True)
                 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')
 
+    # save linear features as numpy file
+    fn = './map_layers/' + case_name + '/linear_feats'
+    np.save(fn,viz_map)
+    print("Created new LF file at %s" % fn)
+    # save inac layers as numpy file
+    fn = './map_layers/' + case_name + '/inac_layers'
+    np.save(fn,viz_map_inac)
+    print("Created new inac file at %s" % fn)
+
 
     # save terrain as csv file (this method is pretty slow, but can compensate with interp)
-    [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
-
-    if save_files:
-        if save_to_folder:
-            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"
-            np.savetxt(elv_filename,e_interp,delimiter=",", fmt='%f')
+    if(get_elev):
+        print("Getting terrain map, this will 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
+        if save_files:
+            if save_to_folder:
+                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"
+                np.savetxt(elv_filename,e_interp,delimiter=",", fmt='%f')
 
+    #Save plots of layers
+    plot_dir = './map_layers/' + case_name + '/'
+    plot_name = ''
     if plot_data:
-        plt.imshow(e_interp)
-        plt.show()
+        if(get_elev):
+            t = 'Elevation ' + case_name
+            plt.title(t)
+            plt.imshow(e_interp)
+            plot_name = plot_dir + 'elevation_' + case_name
+            plt.savefig(plot_name)
+            plt.close()
+            #plt.show()
 
         # plt_list = []
         # fix stupid names
@@ -294,9 +328,20 @@ def grab_features(anchor_point, extent, sample_dist = 10, case_name = 'blah', he
             # flip y values
             # col_idx = viz_map.shape[0] - col_idx
             # plt_list.append(go.Scatter(x=row_idx, y=col_idx, mode='markers', name=name_list[i]))
+            t = name_list[i] + " " + case_name
+            plt.title(t)
             plt.imshow(viz_map[:,:,i])
-            # plt.title(name_list[i])
-            plt.show()
+            plot_name = plot_dir + name_list[i] + '_features'
+            plt.savefig(plot_name)
+            plt.close()
+
+        for i in range(viz_map_inac.shape[-1]):
+            t = "Inacessible " + inac_layers[i] + " " + case_name
+            plt.title(t)
+            plt.imshow(viz_map_inac[:,:,i])
+            plot_name = plot_dir + 'inac_' + inac_layers[i]
+            plt.savefig(plot_name)
+            plt.close()
 
         # fig = go.Figure(data=plt_list)
         # fig.show()

lpm_maps.py

 import numpy as np
 import matplotlib
 import matplotlib.pyplot as plt
+import os
+from os import path, getcwd
 from scipy.cluster.vq import kmeans, whiten
 from scipy.io import loadmat
 import pandas as pd
+import arcgis.geometry
+from arcgis.gis import GIS
+from arcgis.features import FeatureLayer
 
 from arcgis_terrain import get_terrain_map, lat_lon2meters
 from feature_set import grab_features
@@ -12,6 +17,66 @@ from feature_set import grab_features
 #Script made to create a set of SAR datasets.
 
 
+#Used to test map layers, make sure we can access content.
+def gis_test():
+
+    #Lost Person Model content
+    roads_url = "https://carto.nationalmap.gov/arcgis/rest/services/transportation/MapServer/30"
+    river_url = "https://hydro.nationalmap.gov/arcgis/rest/services/nhd/MapServer/6"
+    riverw_url = "https://hydro.nationalmap.gov/arcgis/rest/services/nhd/MapServer/8"
+    water_url = "https://hydro.nationalmap.gov/arcgis/rest/services/nhd/MapServer/9"
+    powerlines_url = "https://services1.arcgis.com/Hp6G80Pky0om7QvQ/ArcGIS/rest/services/Electric_Power_Transmission_Lines/FeatureServer/0"
+    railroads_url = "https://carto.nationalmap.gov/arcgis/rest/services/transportation/MapServer/35"
+    trails_url = "https://partnerships.nationalmap.gov/arcgis/rest/services/USGSTrails/MapServer/0"
+    url_list = [riverw_url, river_url, roads_url, water_url, powerlines_url, railroads_url, trails_url]
+
+
+    #ESRI account information...
+    ky = "AAPKe7cb1ab4f2ba44748cf53ac4f30d0caavu5t_-uMdWP8SRjtoea3s66-cRyFZmMTNx4rqy2w5sjwegU_tbJyBesd0LCMmdtV"
+    user = 'blhowell@vt.edu'
+    pw = '5ZWt5%#FEm5&)@v'
+    url = 'https://virginiatech.maps.arcgis.com'
+    c_id = 'WIXoel276NxsLWI5'
+
+
+    #gis = GIS(url, client_id=c_id)
+    #print("Successfully logged in as: " + gis.properties.user.username)
+    gis = GIS(url)
+    #Display map tutorial
+    map = gis.map()
+    map.basemap = "topo-vector"
+    map.center = [34.027, -118.805]
+    map.zoom = 13
+
+    #Load the feature layer
+    river_layer = FeatureLayer(river_url,gis=gis)
+    print(river_layer)
+    print(river_layer.properties.capabilities)
+    map.add_layer(river_layer)
+
+    #See fields
+    for f in river_layer.properties.fields:
+        print(f['name'])
+
+
+    #Now, query it...
+    print('Starting query....')
+    q = river_layer.query(where='lengthkm>50',out_fields='flowdir,gnis_name')
+    print(len(q.features))
+
+
+    #q = river_layer.query(return_count_only=False, return_ids_only=False, return_geometry=True,
+    #          out_sr='3857', result_record_count = 10)
+    print('query finished!!')
+
+    #Get map to export to html file. Unfortunately this is all useless...
+    export_dir = './arcgis_maps/display_test.html'
+    map.export_to_html(export_dir, title="Testing ArcGIS")
+    print("Map created at %s" % export_dir)
+
+    return
+
+
 #Formats data into a csv that can be loaded as a Pandas dataframe.
 #Finds incident keys for LPM fit data
 def incident_list():
@@ -48,14 +113,70 @@ def incident_list():
     #print(new_df.columns)
     return
 
+#Save collected map layers into one file
+def format_terrain():
+
+
+
+    #For LPM, we need a map with inacessible features, linear features, and elevation all of same size
+    #Save as one matrix, (inac, lf, elev)
+
+
+    return
+
 #My way of checking matlab datasets without using matlab
 def investigate():
 
-    dir = './matlab_data/'
+    name_list = ['rivers_bdd', 'rivers', 'roads', 'lakes', 'powerlines', 'railroads', 'trails']
+    inac_layers = ['rivers_bdd', 'lakes']
+
+    ##############  For comparing to previously collected Matlab data   #################
+    dir = '../ags_grabber/matlab_data/'
     savedir = './plots_test/'
     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
+    print(np.shape(elev))
+#############################################################################################
+
+
+    #Load file made by feature_set
+    lf_dir = './map_layers/test2/'
+    lf_file = 'linear_feats.npy'
+    inac_file = 'inac_layers.npy'
+    layers = np.load(lf_dir+lf_file)
+    inac = np.load(lf_dir+inac_file)
+    elev = np.loadtxt('./map_layers/trust/elv_data_1.csv',delimiter=',')
+
+
+    #Testing content of linear feature file
+    print("Layer matrix shape = ")
+    print(np.shape(layers))
+    print(layers[0][0][0])
+
+    print("Inacessible matrix shape = ")
+    print(np.shape(inac))
+    print(inac[0][0][0])
+
+    print('Elev matrix shape = ')
+    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
 
@@ -63,17 +184,18 @@ def investigate():
 #Also make images to show where areas are.
 def collect_terrain():
 
-    extent_download = 100000 #Size of maps to dowload, in meters (big...)
+    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 = 2          #Minimum number of points needed for a large download
+    count_thres = 3          #Minimum number of points needed for a large download
 
     #Parameters for feature_set
     res = 25
-    folder = 'trust'
+    folder = 'init'
 
     #Collection of map centers, in meters. Use to check if we need to download a new area
     ipp_list = []
+    ipp_lat_list = []
 
     #Used to group search incidents in the same area together. Added to DF at end
     keys = []
@@ -131,6 +253,7 @@ def collect_terrain():
         #See if we need to download a new map
         if(len(ipp_list) == 0):
             ipp_list.append(ipp_xy)
+            ipp_lat_list.append(ipp_point)
             keys.append(len(ipp_list)-1)
         else:
             #Find nearest downloaded center point
@@ -147,34 +270,53 @@ def collect_terrain():
             if(not(topedge and botedge and rightedge and leftedge)):
                 #print("Downloading new map")
                 ipp_list.append(ipp_xy)
+                ipp_lat_list.append(ipp_point)
                 keys.append(len(ipp_list)-1)
             #Add key to incident locations so we know the general area
             else:
                 keys.append(closest)
     print("Found %d maps to collect." % len(ipp_list))
     print(keys)
-    print(ipp_list)
-    print(ipp_list[0][0])
+    #print(ipp_list)
+    #print(ipp_list[0][0])
 
     #Add group keys to dataframe
     incidents.insert(0,'area',keys)
 
+
+    #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_[]'
+
+    #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
+
+
+
     #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 = 1    #Maximum number of things to download
+    max_down = 3    #Maximum number of things to download
     count = 0
     do_plot = True
-    for i in range(0,len(ipp_list)):
-        print("Collecting GIS data %d/%d" % (i,len(ipp_list)))
+    do_elev = True
+    for i in range(0,len(ipp_lat_list)):
+        print("Collecting GIS data %d/%d" % (i,len(ipp_lat_list)))
         #Download local map
         if(keys.count(i) <= count_thres):
             size = extent_test
-        else:
-            size = extent_download
             if(count < max_down):
                 count = count + 1
-                grab_features(ipp_list[i], extent_download, sample_dist = res, case_name = folder, save_files = True, save_to_folder = True, file_id = str(count), plot_data=do_plot)
-
+                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
 
     #For visualizing
     draw = True
@@ -226,7 +368,10 @@ def main():
     extent = 3000   #size in meters of map area
 
     #incident_list()
+    #gis_test()
     collect_terrain()
+    #investigate()
+
 
     return