# Final Project: To find suitable land to create vacation home # Class: GEO376 # File name: hp_finalproject.py # Created by: Helen Peng # Contact info: huierppp@yahoo.com ''' Brief Description: the Python script is to create vacation-home footprints in a suitable land where can see rivers and access by roads. It divides into two parts: one is to find suitable lands, the other is to create footprints polygons in the suitable land. In the first part, it defines three functions. One function is to select proper roads and buffer roads with a certain distance. The second function is to run viewshed analysis for a stream data and convert raster data to polygons. The third function is to intersect road buffer and viewshed polygon, then select polygons with larger areas by calculating areas In the second part, it reads centroids from the suitable land polygon shapefile, and create footprints polygons from centroids. ''' # Import system modules import sys, string, os, arcgisscripting, traceback # Create the Geoprocessor object gp = arcgisscripting.create() # Check out any necessary licenses gp.CheckOutExtension("spatial") # Set toolboxes gp.AddToolbox("C:/Program Files/ArcGIS/ArcToolbox/Toolboxes/Spatial Analyst Tools.tbx") gp.AddToolbox("C:/Program Files/ArcGIS/ArcToolbox/Toolboxes/Spatial Statistics Tools.tbx") gp.AddToolbox("C:/Program Files/ArcGIS/ArcToolbox/Toolboxes/Conversion Tools.tbx") gp.AddToolbox("C:/Program Files/ArcGIS/ArcToolbox/Toolboxes/Analysis Tools.tbx") # Path and workspace gp.workspace = "C:\\temp\\GEO376\\Final_project\\" # Set data paths input_path = "C:\\temp\\GEO376\\Final_project\\" output_path = "C:\\temp\\GEO376\\Final_project\\output\\" # Set input data variables boundary_clip = input_path + "boundary_clip.shp" tn_roads = input_path + "tn_roads.shp" tn_streams = input_path + "tn_streams.shp" dem1 = input_path + "ned_clip1" # Set output variables roads_shp = output_path + "roads_clip.shp" road_select = output_path + "road_select.shp" roads_buffer = output_path + "roads_buffer.shp" streams_clip = output_path + "streams_clip.shp" streams_select = output_path + "streams_select.shp" viewshed_strm = output_path + "viewshed_strm" setnull_view = output_path + "setnull_view" reclass_view = output_path + "reclass_view" viewshed_poly = output_path + "viewshed_poly.shp" roads_buffer = output_path + "roads_buffer.shp" view_intersect = output_path + "view_intersect.shp" view_calculateareas = output_path + "view_calculateareas.shp" # Set variable for new footprints shapefile inputShp = output_path + "view_select1.shp" outputShp = output_path + "footprints1.shp" # Open a new text file to write to txtFile = open("C:\\temp\\GEO376\\Final_project\\output\\footprints_coords1.csv", "w") #if the output shapefile exists, delete it if gp.Exists(outputShp): gp.Delete_Management(outputShp) # Define function to clip and select proper roads, then buffer roads with a certain distance def RoadBuffer(boundary_clip, tn_roads, roads_buffer, rds_select_querry, buffer_distance): # Process: Clip gp.Clip_analysis(tn_roads, boundary_clip, roads_shp, "") # Process: Select gp.Select_analysis(roads_shp, road_select, rds_select_querry) # Process: Buffer gp.Buffer_analysis(road_select, roads_buffer, buffer_distance, "FULL", "ROUND", "NONE", "") print "Finished road buffer" # Define function to run viewshed analysis for a stream data and convert raster data to polygons def StreamViwshed(boundary_clip, tn_streams, dem1, strm_select_querry, viewshed_poly): # Process: Clip... gp.Clip_analysis(tn_streams, boundary_clip, streams_clip, "") # Process: Select... gp.Select_analysis(streams_clip, streams_select, strm_select_querry) # Process: Viewshed... gp.Viewshed_sa(dem1, streams_select, viewshed_strm, "1", "FLAT_EARTH", "") print "Finished viewshed analysis" # Process: Set Null... gp.SetNull_sa(viewshed_strm, viewshed_strm, setnull_view, "VALUE =0") # Process: Reclassify... gp.Reclassify_sa(setnull_view, "VALUE", "1 180 1", reclass_view, "DATA") # Process: Raster to Polygon... gp.RasterToPolygon_conversion(reclass_view, viewshed_poly, "SIMPLIFY", "VALUE") print "Convert to polygon" # Define function to intersect road buffer and viewshed polygon, then select polygons with larger areas by calculating areas def IntersectCalSel (roads_buffer, viewshed_poly, area_querry, inputShp) # Process: Intersect... gp.Intersect_analysis(roads_buffer + ';' + viewshed_poly, view_intersect, "ALL", "", "INPUT") print "intersect finished" # Process: Calculate Areas gp.CalculateAreas_stats(view_intersect, view_calculateareas) print "calculate area" # Process: Select polygons with larger area gp.Select_analysis(view_calculateareas, inputShp, area_querry) print "select suitable land" try: # Call RoadBuffer function to buffer selected roads rds_select_querry = "FCC = 'A20' OR FCC = 'A30' OR NAME ='GREAT SMOKY MTNS %'" buffer_distance = "150 Meters" RoadBuffer(boundary_clip, tn_roads, roads_buffer, rds_select_querry, buffer_distance) # Call function of viewshed analysis for a stream data strm_select_querry = "\"STREETNAME\" LIKE '% Creek' OR \"STREETNAME\" = 'Little Pigeon River'" StreamViwshed(boundary_clip, tn_streams, dem1, strm_select_querry, viewshed_poly) # Call function of IntersectCalSel for roads_buffer, viewshed_poly area_querry = "\"F_AREA\" >= 7.43E-06" IntersectCalSel (roads_buffer, viewshed_poly, area_querry, inputShp) # Create a new feature class gp.CreateFeatureClass_management(os.path.dirname(outputShp), os.path.basename(outputShp), "Polygon", "", "", "", inputShp) # Create search cursor cur = gp.SearchCursor(inputShp) row = cur.Next() # Describe the new feature class desc = gp.Describe(outputShp) shpFld = desc.ShapeFieldName # Create Array, Point, and Cursor objects pnt = gp.CreateObject("Point") curr = gp.InsertCursor(outputShp) a = 1 # Enter while loop for each feature/row while row: # read a centroid from a polygon provided feat = row.shape centroid = feat.centroid print "The midpoint/centroid is: "+ centroid #edited--split string, get x,y values ctr = centroid.split(' ') x = float(ctr[0]) y = float(ctr[1]) print a, x, y #create array and draw the square from centroids #write x, y to excel/textfile array = gp.CreateObject("Array") pnt.x = x - 0.001 pnt.y = y - 0.001 array.Add(pnt) txtFile.write(str(a) + "," + str(x) + "," + str(y)+ "\n") pnt.x = x - 0.001 pnt.y = y + 0.001 array.Add(pnt) txtFile.write(str(a) + "," + str(x) + "," + str(y)+ "\n") pnt.x = x + 0.001 pnt.y = y + 0.001 array.Add(pnt) txtFile.write(str(a) + "," + str(x) + "," + str(y)+ "\n") pnt.x = x + 0.001 pnt.y = y - 0.001 array.Add(pnt) txtFile.write(str(a) + "," + str(x) + "," + str(y)+ "\n") pnt.x = x - 0.001 pnt.y = y - 0.001 array.Add(pnt) rowr = curr.NewRow() rowr.Shape = array curr.InsertRow(rowr) print "polygon done" del array a = a + 1 row = cur.Next() except: tb = sys.exc_info()[2] tbinfo = traceback.format_tb(tb)[0] pymsg = "PYTHON ERRORS:\nTraceback Info:\n" + tbinfo + "\nError Info:\n " + str(sys.exc_type) + ": " + str(sys.exc_value) + "\n" msgs = "GP ERRORS:\n" + gp.GetMessages(2) + "\n" gp.AddError(msgs) gp.AddError(pymsg) print msgs print pymsg gp.AddMessage(gp.GetMessages(1)) print gp.GetMessages(1)