90 likes | 206 Vues
This project, submitted by Robert Abbott for GEOG 375, utilizes ArcPy to manage and analyze geospatial data. It begins by allowing user input for shapefiles and setting local variables for input and output features. The script checks for existing shapefiles and deletes them as needed. Using spatial join functionality, the project combines various shapefiles to create output features while adding specific fields including GIS_ID, latitude, and longitude. Error handling is integrated to ensure smooth execution, with proper field management and cleanup of temporary files.
E N D
GEOG 375 Final Project Robert Abbotts Spring 2013
Code Explanation Set Input shapefile from user input • inFeatures = arcpy.GetParameterAsText(0) Set local variable to required shapefiles • join_features = "//sac/GIS_Script_Shapefiles/24kquad_trim.shp" • plss="//sac/GIS_Script_Shapefiles/pls_fill.shp" • county="//sac/GIS_Script_Shapefiles/Counties.shp" Get input folder • inputfolder = os.path.dirname(inFeatures) Check if shapefiles exist, if they are, delete them • if os.path.isfile(inputfolder+"\out.shp"): • arcpy.Delete_management(inputfolder+"\out.shp") • if os.path.isfile(inputfolder+"\outcounty.shp"): • arcpy.Delete_management(inputfolder+"\outcounty.shp") • if os.path.isfile(inputfolder+"\PExport.shp"): • arcpy.Delete_management(inputfolder+"\PExport.shp") Set local variables for shapefiles • out_feature_class = inputfolder+"\out.shp" • out_county_class = inputfolder+"\outcounty.shp" • outFeatures = inputfolder+"\PExport.shp"
Code Explanation Use Spatial join geoprocessing tools to create output shapefiles • arcpy.SpatialJoin_analysis(inFeatures, join_features, out_feature_class) • arcpy.SpatialJoin_analysis(out_feature_class, plss,out_county_class) • arcpy.SpatialJoin_analysis(out_county_class, county, outFeatures)
Code Explanation Set local variables for table fields • fieldName1 = "GIS_ID" • fieldLength = 16 • fieldName2 = "DDLAT" • fieldPrecision2 = 9 • fieldName3 = "DDLON" • fieldPrecision3 = 9 • fieldName4 = "TEALE_X" • fieldPrecision4 = 9 • fieldName5 = "TEALE_Y" • fieldPrecision5 = 9
Code Explanation Add GIS_ID, DDLAT, DDLON, TEALE_X, TEALE_Y fields • arcpy.AddField_management(outFeatures, fieldName1, "TEXT", "", "", fieldLength) Calculate fields for update data to dispplay date, GPS date and GPS person in the table arcpy.CalculateField_management(outFeatures, fieldName1, • '"%s%s%s%s" % (!Datafile![0:8], !GPS_Date![-4:-1],!GPS_Date![-1], !GPS_Person!)', • "PYTHON") • arcpy.AddField_management(outFeatures, fieldName2, "DOUBLE", "", "","","","NULLABLE","") • arcpy.AddField_management(outFeatures, fieldName3, "DOUBLE", "", "","","","NULLABLE","") • arcpy.AddField_management(outFeatures, fieldName4, "DOUBLE", "", "","","","NULLABLE","") • arcpy.AddField_management(outFeatures, fieldName5, "DOUBLE", "", "","","","NULLABLE","") Add X Y fields and calculate fields for X, Y from Point_X and Point_Y • arcpy.AddXY_management(outFeatures) • arcpy.CalculateField_management(outFeatures, fieldName4, • '!POINT_X!', "PYTHON") • arcpy.CalculateField_management(outFeatures, fieldName5, • '!POINT_Y!', "PYTHON")
Code Explanation Convert X Y to Lat/Lon using update cursor using for loop • desc = arcpy.Describe(outFeatures) • shapefieldname = desc.ShapeFieldName • try: • row, rows = None, None • rows = arcpy.UpdateCursor(outFeatures, r'', \ • r'GEOGCS["GCS_WGS_1984",' + \ • 'DATUM["D_WGS_1984",' + \ • 'SPHEROID["WGS_1984",6378137,298.257223563]],' + \ • 'PRIMEM["Greenwich",0],' + \ • 'UNIT["Degree",0.017453292519943295]]') • for row in rows: • feat = row.getValue(shapefieldname) • pnt = feat.getPart() • row.DDLat = pnt.Y • row.DDLon = pnt.X • rows.updateRow(row)
Code Explanation Except clause if error occurred and delete table rows and temp shapefiles • except: • if not arcpy.GetMessages() == "": • arcpy.AddMessage(arcpy.GetMessages(2)) • finally: • # Regardless of whether the script succeeds or not, delete • # the row and cursor • # • if row: • del row • if rows: • del rows Delete fields and shapefile • dropFields = ["POINT_X", "POINT_Y", "Join_Count", "TARGET_FID","Join_Cou_1","TARGET_F_1","Join_Cou_2","TARGET_F_2","AREA","PERIMETER", "DRGINDEX_","DRGINDEX_I", "OCODE", "USGS100","USGS250","CD","AREA_1","PERIMETE_1","PLSFILL_","PLSFILL_ID","X_COORD","Y_COORD","AREA_12","PERIMETE_2","COUNTY_","COUNTY_ID","Range_1","Staff"] • arcpy.DeleteField_management(outFeatures, dropFields) • arcpy.Delete_management(inputfolder+"\out.shp") • arcpy.Delete_management(inputfolder+"\outcounty.shp") Get the map document mxd • mxd = arcpy.mapping.MapDocument("CURRENT") Get the data frame • df = arcpy.mapping.ListDataFrames(mxd,"*")[0] Create a new layer • newlayer = arcpy.mapping.Layer(outFeatures) Add the layer to the map at the bottom of the TOC in data frame 0 • arcpy.mapping.AddLayer(df, newlayer,"TOP") Refresh arcmap view, table of content • arcpy.RefreshActiveView() • arcpy.RefreshTOC() • del mxd, df, newlayer
Output • outFeatures is PExport.shp • If Arcmap is open when the script is run using ArcToolbar, a temp new layer (outFeatures) will be displayed on top of table of content in ArcMap.
Summary • First, the python script is to get user input shapefile from arctoolbox prompt. • Second, the input shapefiles join other required shapefiles to create result output shapefile. • Third, check if empty or none data fields exist in table using searchcursor function. • Fourth, using add fields (add GIS_ID, DDLAT, DDLON, TEALE_X, TEALE_Y fields and calculate fields to populate the table fields with updated data using updatecursor function. • Then, delete temporary shapefiles and table rows. • Finally, add a new layer (outFeatures shapefile) to ArcMap on top of the table of content. And refresh the ArcMap document.