1 / 9

GEOG 375 Final Project

This code explanation demonstrates how to process geospatial data by creating and updating shapefiles using arcpy in Python scripting. Input shapefiles are joined, fields are added, and X, Y coordinates converted to Lat/Lon. Error handling and field deletion are included.

esutherland
Download Presentation

GEOG 375 Final Project

An Image/Link below is provided (as is) to download presentation Download Policy: Content on the Website is provided to you AS IS for your information and personal use and may not be sold / licensed / shared on other websites without getting consent from its author. Content is provided to you AS IS for your information and personal use only. Download presentation by click this link. While downloading, if for some reason you are not able to download a presentation, the publisher may have deleted the file from their server. During download, if you can't get a presentation, the file might be deleted by the publisher.

E N D

Presentation Transcript

  1. GEOG 375 Final Project Robert Abbotts Spring 2013

  2. 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"

  3. 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)

  4. 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

  5. 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")

  6. 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)

  7. 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

  8. 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.

  9. 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.

More Related