POST
|
I have a script that looks in a folder for all the rasters it contains, performs simple operations on those rasters, and saves them as new rasters in a new folder. The script works fine in ArcGIS 10.2 with python version 2.7.3, but gives the following error in ArcGIS 10.3 with python 2.7.8. The script is embedded in a tool, i.e. not standalone in python. I've attached example input files, the script (entitled .NDVI_Calc3_Parameters_W.py), and its associated toolbox. The script fails at: arcpy.CalculateStatistics_management(inputRaster) Troubleshooting has lead me to believe that someone it no longer recognizes the input as a raster. I've tried various workarounds (including removing the build attribute table line before that one, creating a raster object, etc.) but I still can't get it to work. Any ideas? Thank you in advance for any help you provide. Sincerely, Amelie
... View more
04-29-2016
12:23 PM
|
0
|
2
|
2175
|
POST
|
yes, sorry. The text file is just the kernel they talked about here ArcGIS Desktop saved as a .txt. I'm attaching my kernel (it's a bit more complicated than the examples provided in the help). Amelie
... View more
10-15-2015
02:18 PM
|
0
|
1
|
606
|
POST
|
Hi, I'd like to get the first two numbers in the .txt file that is loaded as by NbrWeight (neighborhood object weights) for the focal statistics tool. I want to use them to calculate the area of the rectangle that is the neighborhood. How do I easily access the values? Here is what I have so far: inWeightFile = arcpy.GetParameterAsText(5) if inWeightFile == '#' or not inWeightFile: inWeightFile = "C:\\temp\\PPR\\Test\\bee_kernel_2_5km.txt" # provide a default value if unspecified # Create the Neighborhood Object for the buffer myNbrWeight = arcpy.sa.NbrWeight(inWeightFile) # Get number of cells in the buffer NbrCellsBuffer = 65 * 65 # XXX need to fix this. Should come from first two numbers in inWeightFile # Process: Focal Statistics - within buffer calculate sum of f_rc1_H1 with exponential decay function and divide by number of cells in buffer arcpy.AddMessage("Calculate proportional cover of high floral quality classes within buffer for first floral season") f_rc_f_H1 = (arcpy.sa.FocalStatistics(f_rc1_H1, myNbrWeight, "SUM", "NODATA") * outConstRaster) / NbrCellsBuffer Any help is greatly appreciated. Best, Amelie
... View more
10-15-2015
10:06 AM
|
0
|
4
|
2730
|
POST
|
OMG!! Check this out: http://www.exploringbinary.com/why-0-point-1-does-not-exist-in-floating-point/ http://docs.oracle.com/cd/E19957-01/806-3568/ncg_goldberg.html Armed with this knowledge, I changed the code to this:
import arcpy
a = arcpy.GetParameter(0)
b = arcpy.GetParameter(1)
c = arcpy.GetParameter(2)
total = ((a * 100) + (b * 100) + (c * 100))
arcpy.AddMessage(total)
if total == 100:
arcpy.AddMessage(" all good")
else:
arcpy.AddMessage(" all NOT good")
And now it works!!! yay! My hair stylist thanks you. No, seriously though, THANK YOU for helping me work this out. I thought I was going crazy. Amelie
... View more
06-10-2014
01:18 PM
|
0
|
0
|
231
|
POST
|
Must be some setting on my computer. Still get the same problem: Executing: Script1 0.6 0.3 0.1 Start Time: Tue Jun 10 17:00:52 2014 Running script Script1... 1 all NOT good Completed script Script1... Succeeded at Tue Jun 10 17:00:52 2014 (Elapsed Time: 0.01 seconds) These are the only lines of code: import arcpy
a = arcpy.GetParameter(0)
b = arcpy.GetParameter(1)
c = arcpy.GetParameter(2)
total = a + b + c
arcpy.AddMessage(total)
if total == 1.0:
arcpy.AddMessage(" all good")
else:
arcpy.AddMessage(" all NOT good") I only have three parameters and they are all set up to Double. yet the 0.6, 0.4, 0.1 combo doesn't work it works in yours?
... View more
06-10-2014
01:05 PM
|
0
|
0
|
231
|
POST
|
this is utterly bizarre. I tried if Tot == "1.0": and nothing worked (went into the else loop and sys.exit() ). So then, I kept the code with if Tot == 1.0: and if the values for Spring, Summer, and Fall in the GUI are set to 0.2,0.6,0.2, respectively program runs completely So does 0.4, 0.6, 0 This works too 0.6, 0, 0.4 but not 0.6, 0.3, 0.1 (it goes into the else loop and does sys.exit() ). I'm tearing my hair out.
... View more
06-10-2014
11:59 AM
|
0
|
0
|
812
|
POST
|
Sorry, it's not an error it's just not what I would expect. I've changed "if Tot == 1" to "if Tot == 1.0" If in the tool GUI I set Spring = 1, Summer = 0, and Fall = 0. Then I go into the if statement and the rest of the script runs (i get the "Sum of floral weights is equal to: 1.0 message). This is what I would expect. BUT if in the tool GUI I set Spring = 0.6, Summer = 0.3, and Fall = 0.1, Then I go into the else statement and the tool stops running even though I've added a "catch message" (arcpy.AddError("The sum of floral weights must equal to one, but it is equal to " + str(Tot))) and it tells me that Tot is equal to 1.0. What do you think? Amelie
... View more
06-10-2014
11:37 AM
|
0
|
0
|
812
|
POST
|
Hi! That does work - sorry. Now I have a float as data type but this doesn't work:
Tot = Spring + Summer + Fall
arcpy.AddMessage("Floral Weights are equal to: " + str(Tot))
if Tot == 1:
arcpy.AddMessage("Sum of floral weights is equal to: " + str(Tot))
else:
arcpy.AddError("The sum of floral weights must equal to one.")
sys.exit() # exit out of model if weights don't add up to 1
Any ideas? Amelie
... View more
06-10-2014
10:53 AM
|
0
|
0
|
814
|
POST
|
Hi! I have a python script that runs fine stand alone but when I try to make it GUI based by adding it to a toolbox as a script I get errors which I believe are related to how I am setting the parameters. I have three values that I need to be set as a real number with a decimal but they keep coming in as string. In ArcGIS toolbox, under properties of the model, and the parameters tab, I set those three parameter types to Double (I have also tried String but that doesn't work either). then in python, that part of the code looks like this: # -*- coding: utf-8 -*- ############################################################################ import sys # Import arcpy module import arcpy #from arcpy import env # Set this to False once testing is complete arcpy.env.overwriteOutput = True ############################################################################ # Script arguments # user defined weights for importance of up to three floral seasons Spring = float(arcpy.GetParameter(3)) arcpy.AddMessage("Floral Weight #1 is equal to: " + str(Spring)) if Spring == '#' or not Spring: Spring = 0.3 # provide a default value if unspecified check1 = arcpy.Parameter(Spring) arcpy.AddMessage("Data type is: " + str(check1.datatype)) Summer = float(arcpy.GetParameter(4)) if Summer == '#' or not Summer: Summer = 0.6 # provide a default value if unspecified Fall = float(arcpy.GetParameter(5)) if Fall == '#' or not Fall: Fall = 0.1 # provide a default value if unspecified Tot = Spring + Summer + Fall arcpy.AddMessage("Floral Weights are equal to: " + str(Tot)) if Tot == 1: arcpy.AddMessage("Sum of floral weights is equal to: " + str(Tot)) else: arcpy.AddError("The sum of floral weights must equal to one.") sys.exit() # exit out of model if weights don't add up to 1 ############################################################################ I've tried GetParameter and GetParameterAsText but the value always seems to be coming in as a string. I know this is probably something very basic but I just can't figure it out. I'm using ArcGIS 10.2. Thank you in advance for any assistance. Sincerely, Amelie
... View more
06-10-2014
08:37 AM
|
0
|
13
|
1147
|
POST
|
Thank you!! Any suggestions on how to make it run faster by any chance? I can definitely comment out the out.save throughout - I was just keeping them while testing the model but can't think of anything else... In any case it now works in ArcGIS 10.1 so big thank you.
... View more
12-14-2012
08:36 AM
|
0
|
0
|
163
|
POST
|
Hi all, I wrote a script which works in 10.0 (which I run on my computer) but doesn't work on my collaborators machine which runs ArcGIS 10.1. I don't have access to a version of 10.1 so it is hard to trouble shoot but I was wondering if anyone has run into this problem before (and tell me how they fixed it) or could just spot the problem in the code below. The script is part of a toolbox I sent him. He had to manually point to where the script was stored under script > properties > source tab and he gets the following error: " 000714 : Error in script <value>. " Does anyone see off the top of their head what the problem is? It's just very weird that the script works in 10.0 but not in 10.1. Also can anyone give me pointers on how to get this script to run faster? This second question is unrelated to the first but I figured I ask since I was already here and its something that has been bugging me. thank you in advance for your time, Amelie ############ 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, "")
... View more
12-13-2012
04:41 PM
|
0
|
2
|
352
|
POST
|
OK I think I figured out how to do it so I figure I would post the answer (even if I did start this thread...) I'm not an expert coder so there probably is a much more elegant/faster way to do this but at least it works. Hopefully this helps someone.
# ---------------------------------------------------------------------------
# Dist2LandCover_v1.py
# Created on: 2012-05-09
# By Amelie Davis
# Description:
# ---------------------------------------------------------------------------
############ IMPORTANT ########################
# Assumes that your land cover layer is projected AND cell size is 30 meters
# Uses NLCD 2006 classification codes
# output file must be a shapefile and the output name must have .shp at the end of it.
##############################################
# 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_calc_dist_from = arcpy.GetParameterAsText(3)
#Zone_value_you_want_to_calc_dist_from = "11" # provide a default value if unspecified
LC_codes = arcpy.GetParameterAsText(4)
Final_Output_shp = Output_Workspace + "\\" + arcpy.GetParameterAsText(5)
####################
# Process: Reclassify #1
#set to 1 LC code from which you want to calculate distance. Example: water = "11"
text_in = Zone_value_you_want_to_calc_dist_from+" "+Zone_value_you_want_to_calc_dist_from+" 1"
fromClassremapList = text_in.split()
fromClassremapList = [int(k) for k in fromClassremapList]
remap = arcpy.sa.RemapRange([fromClassremapList])
nlcd_wat = arcpy.sa.Reclassify(Input_Land_Cover__NLCD_2006_, "VALUE", remap, "NODATA")
#remap = arcpy.sa.RemapRange([[11,11,1],[12,95,"NODATA"]])
#nlcd_wat = arcpy.sa.Reclassify(Input_Land_Cover__NLCD_2006_, "VALUE", remap, "DATA")
out_rc1 = Output_Workspace + "\\out1"
nlcd_wat.save(out_rc1)
# Process: Euclidean Distance
out_euc = arcpy.sa.EucDistance(nlcd_wat, "", "30", "")
out_euc_s = Output_Workspace + "\\euc_dist_1"
out_euc.save(out_euc_s)
arcpy.AddMessage("finished calculating euclidian distance")
# Process: Int
out_Int = arcpy.sa.Int(out_euc)
out_Int_s = Output_Workspace + "\\water_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: Reclassify (2)
#set to 1 LC code for which you want to know the average distance from LC in Process: Reclassify #1
#remap2 = arcpy.sa.RemapRange([[11,94,0],[95,95,1]]) #old
#remap2 = arcpy.sa.RemapRange([[95,95,1]]) #works
LC_CODE = LC_codes.split()
for i in LC_CODE:
text_in2 = i+" "+i+" 1"
forClassremapList = text_in2.split()
forClassremapList = [int(k) for k in forClassremapList]
remap2 = arcpy.sa.RemapRange([forClassremapList])
nlcd_ag = arcpy.sa.Reclassify(Input_Land_Cover__NLCD_2006_, "VALUE", remap2, "NODATA")
out_rc2 = Output_Workspace + "\\out2"
nlcd_ag.save(out_rc2)
arcpy.AddMessage("finished reclassifying")
# Process: Times
# save output permanently?
out_times = arcpy.sa.Times(out_Int, nlcd_ag)
out_temp2 = Output_Workspace + "\\temp_ag"
out_times.save(out_temp2)
field_name = "XX_" + str(i)
# Process: Add field to initial fishnet shapefile
arcpy.AddField_management(Final_Output_shp, field_name, "SHORT", 6, "", "","refcode", "NULLABLE")
# Process: Calculate Field
arcpy.CalculateField_management(Final_Output_shp, field_name, "!FID!","PYTHON_9.3","")
arcpy.AddMessage("Starting zonal statistics")
# Process: Zonal Statistics as Table
outTable = Output_Workspace + "\\dist_per_cell.dbf"
outZSaT = arcpy.sa.ZonalStatisticsAsTable(Final_Output_shp, field_name, out_times, outTable, "DATA", "ALL")
arcpy.AddMessage("finished calculating mean +/- std distance per cell")
# Process: Join Field
fieldList = ["MEAN", "STD", "MIN"]
arcpy.JoinField_management(Final_Output_shp, field_name, outTable, field_name,fieldList)
arcpy.AddMessage("finished with land cover" + str(i))
# Delete intermediary files
##arcpy.Delete_management(out_rc1, "")
##arcpy.Delete_management(out_Int_s, "")
##arcpy.Delete_management(out_rc2, "")
##arcpy.Delete_management(out_euc_s, "")
... View more
05-11-2012
08:29 AM
|
0
|
0
|
174
|
POST
|
Hi! I have a raster (NLCD2006) and I would like to calculate the distance of the nearest pixel of one class from another class. I am struggling to think through what would be the best approach. I thought about reclassifying all the ag and urban pixels to one, setting the rest to NoData and running the euclidean distance tool but that isn't quite what I want since it might give me the nearest a pixel rather than the urban one. Then I thought I could turn the ag pixels to points, same for the urban pixels but in a separate shapefile and then run the near tool but I have too many cells/points for this approach (I'm running this for the entire US for 30 by 30m cells. I can't think of another approach? Does anyone have any suggestions? Any help is greatly appreciated, Amelie
... View more
05-05-2012
08:26 PM
|
0
|
1
|
552
|
POST
|
Fixed it - problem was not with my while loops NOR with the delete row function but rather with my logic and incorrect use of FIDs... Here is the correct code:
# ---------------------------------------------------------------------------
# IterativeNearAnalysis.py
# Created on: thu oct 19 2011 by Amelie Davis
# Takes all shapefiles in a folder
# if points are identical in terms of lat long, removes one point at random
# checks the distance between all points in that same shapefile is less than a specified distance
# deletes at random one of the two 'offending points'
# creates a new shapefile with only the points that are wanted
# checks the output
# final output is a .dbf of the attribute table for those points
# ---------------------------------------------------------------------------
# Import system modules
import arcpy
from arcpy import env
import os
from os import path
# input shapfile
inputfile = arcpy.GetParameterAsText(0)
# default inputfile = "d:/data/leap540/52spp_withid_1000m/BCN_Birds_P_top52spamecro.shp"
#distance to look for points
dist2 = arcpy.GetParameterAsText(1)
dist = dist2 + " Meters"
# default dist2 = 500
#workspace
wksp = arcpy.GetParameterAsText(2)
arcpy.AddMessage("Workspace is "+ wksp)
#set the workspace
arcpy.env.workspace = wksp
(dirName, fileName) = os.path.split(inputfile)
# outputfiles
# Name of output shapefile from which points that are located at the exact same spot are removed
out1 = fileName[:-4]+"_copy.shp"
# make a copy of inputfile and name it out1
arcpy.Select_analysis(inputfile,out1,"#")
# Names of output tables and shapefiles
out1_a = out1[:-4]+"_sort.shp"
out2 = out1_a[:-4] + dist2 + "_near.dbf"
out3 = out2[:-4] + "_sort2.dbf"
# final output
out4 = arcpy.GetParameterAsText(3)
# check output
out5 = out3[:-4] + "_check.dbf"
arcpy.AddField_management(out1,"RANDOM","SHORT","5","#","#","#","NULLABLE","NON_REQUIRED","#")
codeblock1= "dim max, min\nmax=99999\nmin=0\nx=(Int((max-min+1)*Rnd+min))"
arcpy.CalculateField_management(out1,"RANDOM","x","VB",codeblock1)
arcpy.AddField_management(out1,"FID_COPY","SHORT","6","#","#","#","NULLABLE","NON_REQUIRED","#")
arcpy.CalculateField_management(out1,"FID_COPY","[FID]","VB","#")
# delete records that are identical in latitute and longitude
# make sure the point to delete is selected at random by sorting by random
arcpy.Sort_management(out1,out1_a,"RANDOM ASCENDING","#")
fields = ["LATITUDE","LONGITUDE"]
arcpy.DeleteIdentical_management(out1_a,fields,"#","0")
# generate list of points that are within a certain distance (specified by 'dist') or each other
arcpy.GenerateNearTable_analysis(out1_a,out1_a,out2,dist,"NO_LOCATION","NO_ANGLE","ALL","0")
arcpy.AddField_management(out2,"RANDOM2","SHORT","5","#","#","#","NULLABLE","NON_REQUIRED","#")
arcpy.CalculateField_management(out2,"RANDOM2","x","VB",codeblock1)
arcpy.Sort_management(out2,out3,"RANDOM2 ASCENDING","#")
arcpy.DeleteIdentical_management(out3,"NEAR_DIST","#","0")
rows = arcpy.SearchCursor(out3)
i_fid = "NEAR_FID"
a = []
arcpy.MakeFeatureLayer_management(out1_a,"grumble")
row = rows.next()
while row:
i_fid1 = row.getValue(i_fid)
arcpy.AddMessage("FID of row to delete is "+ str(i_fid1))
rows2 = arcpy.UpdateCursor(out1_a)
row2 = rows2.next()
while row2:
i_fid2 = row2.getValue("FID")
cute = row2.getValue("L_ID")
#arcpy.AddMessage("FID of row we are on is "+str(i_fid2))
if i_fid1 == i_fid2:
b = "L_ID = " + str(cute)
a.append(b)
arcpy.AddMessage("deleting row"+ str(i_fid2)+" and original loc is " + str(cute))
row2 = rows2.next()
row = rows.next()
#arcpy.AddMessage("out of loop")
del row
del row2
cntr = 0
while cntr<len(a):
f_query = a[cntr]
arcpy.SelectLayerByAttribute_management("grumble","ADD_TO_SELECTION",f_query)
arcpy.AddMessage(f_query)
cntr +=1
arcpy.SelectLayerByAttribute_management("grumble","SWITCH_SELECTION","#")
arcpy.Select_analysis("grumble",out4,"#")
arcpy.GenerateNearTable_analysis(out4,out4,out5,dist,"NO_LOCATION","NO_ANGLE","ALL","0")
# delete out2 and out3
... View more
10-21-2011
07:44 PM
|
0
|
0
|
412
|
Online Status |
Offline
|
Date Last Visited |
11-11-2020
02:24 AM
|