############ IMPORTANT ######################## # Assumes that your land cover layer is projected in meters # Uses NLCD 2006 classification codes # output file must be a shapefile and the output name must have .shp at the end of it. # The Zone_value (parameter 3) must match what you reclass as one in parameter 4 ############################################## # Set the necessary product code # import arcinfo # Import arcpy module import arcpy from arcpy import env # Check out any necessary licenses arcpy.CheckOutExtension("spatial") ###################### Script arguments Input_Fishnet_with_unique_ID2_field = arcpy.GetParameterAsText(0) #Input_Fishnet_with_unique_ID2_field = "D:\\DATA\\USA\\MigBirdFishnet\\NLCDfishnet20mi_p.shp" # provide a default value if unspecified Input_Land_Cover__NLCD_2006_ = arcpy.GetParameterAsText(1) #Input_Land_Cover__NLCD_2006_ = "D:\\DATA\\NLCD\\nlcd2006_cook" # provide a default value if unspecified Output_Workspace = arcpy.GetParameterAsText(2) #Output_Workspace = "D:\\DATA\\DELETE" # provide a default value if unspecified # Set environment settings env.workspace = Output_Workspace env.mask = Input_Land_Cover__NLCD_2006_ Zone_value_you_want_to_shrink__LC_code_ = arcpy.GetParameterAsText(3) #Zone_value_you_want_to_shrink__LC_code_ = "11" # provide a default value if unspecified Land_Cover_Code_for_which_you_want_to_extract_distance_measure = arcpy.GetParameterAsText(4) #Land_Cover_Code_for_which_you_want_to_extract_distance_measure = "95" # provide a default value if unspecified Final_Output_shp = Output_Workspace + "\\" + arcpy.GetParameterAsText(5) #################### # could try this: tmp3 = Reclassify(Int(EucDistance(inR, "", "1000", "")), "VALUE", "0 NODATA", "DATA") #Get the minimum value in the land cover raster MinResult = arcpy.GetRasterProperties_management(Input_Land_Cover__NLCD_2006_, "MINIMUM") #Get the elevation standard deviation value from geoprocessing result object min_r = MinResult.getOutput(0) #arcpy.AddMessage("minimum value of land cover raster is " + str(min_r)) #Get the maximum value in the land cover raster MaxResult = arcpy.GetRasterProperties_management(Input_Land_Cover__NLCD_2006_, "MAXIMUM") #Get the elevation standard deviation value from geoprocessing result object max_r = MaxResult.getOutput(0) #arcpy.AddMessage("maximum value of land cover raster is " + str(max_r)) # Process: Reclassify, i.e. extract wetland=95 from NLCD value = int(Land_Cover_Code_for_which_you_want_to_extract_distance_measure) if value == int(min_r): remapList1 = [value,value,1] remapList2 = [value+1,max_r,"NODATA"] complete = [remapList1,remapList2] elif value == int(max_r): remapList1 = [min_r,value-1,"NODATA"] remapList2 = [value,value,1] complete = [remapList1,remapList2] elif value >= int(min_r) and value <= int(max_r): remapList1 = [min_r,value-1,"NODATA"] remapList2 = [value,value,1] remapList3 = [value+1,max_r,"NODATA"] complete = [remapList1,remapList2,remapList3] else: arcpy.AddMessage("need to pick a value within range of " + str(min_r) +" to "+ str(max_r)) remap = arcpy.sa.RemapRange(complete) nlcd_ag = arcpy.sa.Reclassify(Input_Land_Cover__NLCD_2006_, "VALUE", remap, "NODATA") out_rc1 = Output_Workspace + "\\out1" nlcd_ag.save(out_rc1) del value del complete # Process: Reclassify again, i.e. extract water=11 from NLCD (to build shoreline) value = int(Zone_value_you_want_to_shrink__LC_code_) if value == int(min_r): remapList1 = [value,value,1] remapList2 = [value+1,max_r,0] complete = [remapList1,remapList2] elif value == int(max_r): remapList1 = [min_r,value-1,0] remapList2 = [value,value,1] complete = [remapList1,remapList2] elif value >= int(min_r) and value <= int(max_r): remapList1 = [min_r,value-1,0] remapList2 = [value,value,1] remapList3 = [value+1,max_r,0] complete = [remapList1,remapList2,remapList3] else: arcpy.AddMessage("need to pick a value within range of " + str(min_r) +" to "+ str(max_r)) remap5 = arcpy.sa.RemapRange(complete) nlcd_line = arcpy.sa.Reclassify(Input_Land_Cover__NLCD_2006_, "VALUE", remap5, "DATA") out_rc5 = Output_Workspace + "\\out5" nlcd_line.save(out_rc5) # Process: Shrink zoneSet = [int(Zone_value_you_want_to_shrink__LC_code_)] out_shrink = arcpy.sa.Shrink(Input_Land_Cover__NLCD_2006_, 1, zoneSet) out_shrink_s = Output_Workspace + "\\shrink1" out_shrink.save(out_shrink_s) # Process: Reclassify (2) nlcd_water = arcpy.sa.Reclassify(out_shrink, "VALUE", remap5, "DATA") out_rc2 = Output_Workspace + "\\out2" nlcd_water.save(out_rc2) # Process: Plus out_plus_for = arcpy.sa.Plus(nlcd_line, nlcd_water) out_plus = Output_Workspace + "\\temp1" out_plus_for.save(out_plus) # Process: Reclassify (3) remap3 = arcpy.sa.RemapRange([[0,0,"NODATA"],[1,1,1],[2,2,"NODATA"]]) nlcd_for = arcpy.sa.Reclassify(out_plus_for, "VALUE", remap3, "DATA") out_rc3 = Output_Workspace + "\\out3" nlcd_for.save(out_rc3) # Process: Reclassify (4) remap4 = arcpy.sa.RemapRange([[0,0,0],[1,1,1000],[2,2,0]]) nlcd_new_pre = arcpy.sa.Reclassify(out_plus_for, "VALUE", remap4, "DATA") out_rc4 = Output_Workspace + "\\out4" nlcd_new_pre.save(out_rc4) arcpy.AddMessage("finished reclassifying") # Process: Plus # save output permanently out_plus2 = arcpy.sa.Plus(nlcd_new_pre,Input_Land_Cover__NLCD_2006_) out_rc6 = Output_Workspace + "\\new_nlcd" out_plus2.save(out_rc6) arcpy.AddMessage("created new nlcd layer") # Process: Euclidean Distance # Probably should make cell_sz_in_raster a model parameter cell_sz_in_raster = 30 out_euc = arcpy.sa.EucDistance(nlcd_for, "", cell_sz_in_raster, "") out_euc_s = Output_Workspace + "\\euc_dist1" out_euc.save(out_euc_s) arcpy.AddMessage("finished calculating euclidian distance") # Process: Times out_times = arcpy.sa.Times(out_euc, nlcd_ag) out_temp2 = Output_Workspace + "\\temp2" out_times.save(out_temp2) # Process: Int # save output permanently out_Int = arcpy.sa.Int(out_times) out_Int_s = Output_Workspace + "\\ag_dist" out_Int.save(out_Int_s) # Process: Copy initial fishnet shapefile so that original is not modified arcpy.Copy_management(Input_Fishnet_with_unique_ID2_field, Final_Output_shp, "") # Process: Add field to initial fishnet shapefile arcpy.AddField_management(Final_Output_shp, "XX_ID_XX", "SHORT", 6, "", "","refcode", "NULLABLE") #arcpy.AddField_management(Final_Output_shp,"XX_ID_XX","SHORT","6","#","#","#","NULLABLE","NON_REQUIRED","#") # Process: Calculate Field #arcpy.CalculateField_management(Final_Output_shp, "XX_ID_XX", "[FID]","VB","X") arcpy.CalculateField_management(Final_Output_shp, "XX_ID_XX", "!FID!","PYTHON_9.3","") # Process: Zonal Statistics as Table outTable = Output_Workspace + "\\ag_dist_per_cell.dbf" outZSaT = arcpy.sa.ZonalStatisticsAsTable(Final_Output_shp, "XX_ID_XX", out_Int_s, outTable, "DATA", "ALL") arcpy.AddMessage("finished calculating mean +/- std distance per cell") # Process: Tabulate Area (2) outTable_T = Output_Workspace + "\\new_nlcd_per_cell.dbf" arcpy.sa.TabulateArea(Final_Output_shp, "XX_ID_XX", out_plus2, "VALUE", outTable_T,"30") arcpy.AddMessage("finished calculating area of each LC class in each cell") # Process: Join Field fieldList = ["MEAN", "STD"] arcpy.JoinField_management(Final_Output_shp, "XX_ID_XX", outTable, "XX_ID_XX",fieldList) # Process: Join Field (2) arcpy.JoinField_management(Final_Output_shp, "XX_ID_XX", outTable_T, "XX_ID_XX","") # Delete intermediary files ##arcpy.Delete_management(out_rc1, "") ##arcpy.Delete_management(out_rc5, "") ##arcpy.Delete_management(out_shrink_s, "") ##arcpy.Delete_management(out_rc2, "") ##arcpy.Delete_management(out_plus, "") ##arcpy.Delete_management(out_rc3, "") ##arcpy.Delete_management(out_rc4, "") ##arcpy.Delete_management(out_euc_s, "") ##arcpy.Delete_management(out_temp2, "")
Solved! Go to Solution.