Latest Contributions by XanderBakker

Hi Krzysztof, Clint already indicated some flaws in the code. Another one is the fact that you do a loop over your features, but nothing is actually done with the feature itself. What you should do in the loop is create a feature layer of the single feature and use that feature layer in the zonal statistics tool. See below an example of how this could be done: import arcpy import os # Set environment settings outWorkspace = r"C:\Project\_Forums\_json\fgdb\crawlAHN2.gdb" arcpy.env.workspace = outWorkspace arcpy.env.overwriteOutput = True # reference to the data fc = os.path.join(outWorkspace, "crowns") raster = os.path.join(outWorkspace, "AHN2r") zoneField = "ID_OK" # temporal output ZS table will be (over)written to in memory workspace outTable = os.path.join("IN_MEMORY", "stat_zs") # name of feature layer (lives in memory) flyr_name = 'a_crown' # get list of zones (this is used for the loop) lst_zones = [row[0] for row in arcpy.da.SearchCursor(fc, (zoneField))] # create an empty dictionary to store the statistics results dct_stats = {} # create a dictionary that holds the statistics type and the corresponding output field names dct_stat_fld = {'MEAN': 'cr_mean',                 'MAX':  'cr_max',                 'MIN': 'cr_min',                 'STD': 'cr_std'} # list of fields in ZS table that hold the statistics I'm interested in ['MEAN','MIN','MAX','STD'] lst_flds = [fld for fld in dct_stat_fld.iterkeys()] # get sa license arcpy.CheckOutExtension("Spatial") # now loop through zones and make a featurelayer with only that zone feature (crown) cnt = 0 for zone in lst_zones:     cnt += 1     # give some feedback on progress     if cnt % 10 == 0:         print "Processing crown: {0} ({1})".format(cnt, zone)     # create the where clause for that crown     where = "{0} = '{1}'".format(arcpy.AddFieldDelimiters(fc, zoneField), zone)     # make the featurelayer of current crown (zone)     arcpy.MakeFeatureLayer_management(fc, flyr_name, where)     # calculate the zonal statistics for that crown     arcpy.sa.ZonalStatisticsAsTable(flyr_name, zoneField, raster, outTable, "NODATA", "ALL")     # read the ZS table and write to dictionary     with arcpy.da.SearchCursor(outTable, lst_flds) as curs:         for row in curs:             i = 0             for fld in lst_flds:                 key = (zone, fld) # looks like ('crown ID', 'statistics type like MIN')                 dct_stats[key] = row                 i += 1 # return sa license arcpy.CheckInExtension("Spatial") # create update field collection (field  for the update cursor) flds = [zoneField] flds.extend([dct_stat_fld[fld] for fld in lst_flds]) # do an update cursor to write the results from dictionary to crown featureclass with arcpy.da.UpdateCursor(fc, flds) as curs:     for row in curs:         zone = row[0]         i = 1         for fld in lst_flds:             key = (zone, fld)             if key in dct_stats:                 row = dct_stats[key]             i += 1         curs.updateRow(row) I hope this is useful for you. Kind regards, Xander ... View more

Dear Kim and Xander, thank you very much for your responses, both worked a treat. I wish the ESRI online help was that clear and succinct. Xander, I'm afraid you've started something. I do have 10.1 (10.2 actually), and I tried using arcpy.da.SearchCursor as you suggested. I followed the esri online help and examples as best I could and ended up with this: datasetname = lyr.dataSource     #this is a shapefile fldname = "FIRECODE" lst_values = ["ARC", "CEA", "DAB", "MGD"] where = "{0} in ('{1}')".format(arcpy.AddFieldDelimiters(lyr1, fldname), "','".join(lst_values)) with arcpy.da.SearchCursor (datasetname,fldname,where) as cur:   for eachrecord in cur:     regcode = eachrecord.getValue(fldname) . . At this point I get the error message "AttributeError: 'tuple' object has no attribute 'getValue'" Any suggestions? Cheers, Damian Hi Damian, You are close. The values are retrieved from the row by using an index. If you only provide 1 field as fields in the da.SearchCursor, your value will be at index 0, hence the line: regcode = eachrecord[0]. Just look at the code below where I have introduced more fields to clarify: datasetname = lyr.dataSource fldname = "FIRECODE" lst_values = ["ARC", "CEA", "DAB", "MGD"] where = "{0} in ('{1}')".format(arcpy.AddFieldDelimiters(lyr1, fldname), "','".join(lst_values)) # define a tuple (or list) of field names to use in the cursor flds = (fldname, 'someOtherFieldname', 'andYetAnotherFieldName') with arcpy.da.SearchCursor(datasetname, flds, where) as cur:     for eachrecord in cur:         regcode = eachrecord[0] # firecode         someOtherValue = eachrecord[1]         andYetAnotherValue = eachrecord[2] Since you are only using 1 field in a search cursor, I assume you are trying to determine some statistic. If that is the case there maybe faster ways to do that (list comprehensions of using a dictionary). Just let me know what you want to do in the loop and I can give you some pointers. Another powerful feature of the da cursors are the build-in fields like "SHAPE@" for the geometry, "OID@" for ObjectID, etc. Kind regards, xander ... View more

Hi Damian, Kim already suggested a valid syntax for your where clause, but I would like to add something for you to consider. Since the syntax of the where clause depends on the dataset, it might be better to write code that is independent of the type of data you are querying. For this purpose arcpy provides the AddFieldDelimiters functionality. Have a look at the code below: fldname = "FIRECODE" lst_values = ["ARC", "CEA", "DAB", "MGD"] where = "{0} in ('{1}')".format(arcpy.AddFieldDelimiters(lyr, fldname), "','".join(lst_values)) cur = arcpy.SearchCursor(lyr, where) In this case I have the field name defined as a variable "fldname" and I also have a list of values to match it to ("lst_values"). I create a variable "where" containing the where clause and pass that variable in the search cursor (BTW if you have access to 10.1 you should try the arcpy.da.SearchCursor which is much faster, but uses a different syntax). Let's take the where clause construction apart to explain what happens: The AddFieldDelimiter adds the quotes, square brackets or nothing depending on the data source: print arcpy.AddFieldDelimiters(lyr, fldname) # >> "FIRECODE" If you use the join method on a string a trow in a list, it will return the list of items where the string is placed between each item: print "','".join(lst_values) # >> ARC','CEA','DAB','MGD If you use the format on a string you can replace the {0}, {1} by values you pass in, allowing better readability. In this case the list is enclosed by curved brackets and single quotes: print "('{0}')".format("','".join(lst_values)) # >> ('ARC','CEA','DAB','MGD') Putting it all together gives: where = "{0} in ('{1}')".format(arcpy.AddFieldDelimiters(lyr, fldname), "','".join(lst_values)) print where # >> "FIRECODE" in ('ARC','CEA','DAB','MGD') Kind regards, Xander ... View more

Hi Anne, That is a relatively simple thing to do with a ArcGIS toolbox tool. I advise you to read the following Help page: http://resources.arcgis.com/en/help/main/10.2/index.html#//00150000000n000000 If you scroll to the very end of this page you will see a header "obtained from". This is what it is all about including the "Filter" option and multi value setting. Create a new toolbox, add a script (name it appropriately), select the location of the script, and in the next screen you have to define your parameters. In this screen select feature layer as a parameter type (let's call it "point features"). In the Filter option (lowe part of dialog) set it to Feature Class and select only point features. The user will not be able to select a layer that isn't a point featureclass. The next parameter would be of the type field (let's call this one "IDW fields"). Now in the lower part of the dialog set the following things: MultiValue: Yes (this means more than 1 field van be chosen) Filter: Field (in the dialoog with field type select only the appropriate field types for IDW: like Float and Double), so only these field will be shown in your script dialog Obtained from: select your "point features". This means that the field will be obtained from the point feature layer the user selected With a few simple clicks you have set up (part of) the scrip dialog. Now when you read the fields inside the script (arcpy.GetParameter() or arcpy.GetParameterAsText()), I advise you to read is as text. It will be a single string. Split is using the ';' as split character to obtain a list of the field names. flds = arcpy.GetParameterAsText(1) lst_flds = flds.split(';') Now you have a list of the chosen field names that can be used in the IDW loop. See attached images to clarify. [ATTACH=CONFIG]32136[/ATTACH] Defining the field parameter [ATTACH=CONFIG]32137[/ATTACH] The result in the script dialog Kind Regards, Xander ... View more

Hi Jeffrey, It's a little difficult to help you without knowing what data was used for the process and what you specified in the dialog. If it's possible to attach a part of you data, what you specified in the dialog and to indicate the version you are using, I can see if I can reproduce it. Otherwise contact Esri support. Kind regards, Xander ... View more

Hi Abdel, I suppose the speed would not matter as long as the line from point Q extent outside the convex hull polygon. So with that in mind, take a look at the Python code below. It works with a list of points (hard coded in this snippet) and an angle (could be arithmetic or geographic). The following steps are carried out: the points are added to a multipoint object the convex hull polygon is calculated from the multipoint the boundary is extracted from the polygon the centroid of the polygon is in this case used for point Q to be sure that the line starting in point Q will exit the polygon, I use a "large" length (so no speed used) the point Q2 is determined using the length and angle a line is created using Q and Q2 the line is intersected with the boundary to obtain the desired intersection point the X and Y of the resulting point is printed def main():     import arcpy, os     # input points     pnts = [[200,800],[100, 500],[300,100],[1000,200],[1200,700]]     # input angle     angle_geo = 100 # geographic degrees     angle_ari = arithmatic_geographic(angle_geo) # = 350 in arithmetic degrees     # create a multi point from input points     arr_mp = arcpy.Array()     for p in pnts:         pnt = arcpy.Point(p[0], p[1])         arr_mp.add(pnt)     mp = arcpy.Multipoint(arr_mp)     # determine convex hull of input points     pol = mp.convexHull()     # get boundary of convex hull polygon     bnd = pol.boundary()     # let's take centroid for Q     q = pol.centroid     # determine line length (enough to intersect boundary)     length = pol.extent.width + pol.extent.height     # create second point (Q2)     dx, dy = LengthDir(length, angle_ari)     q2 = arcpy.Point(q.X + dx, q.Y + dy)     # make line of Q and Q2     arr_line = arcpy.Array()     arr_line.add(q)     arr_line.add(q2)     line = arcpy.Polyline(arr_line)     # intersect (other, dimension), dimension can be 1 or 2     pnt_int = line.intersect(bnd, 1) # returns mp     # print resulting intersection coordinates     for p in pnt_int:         print p.X, p.Y def arithmatic_geographic(angle):     new_angle = angle + 270     while new_angle > 360:         new_angle -= 360     return 360 - new_angle def LengthDir(length, angle):     import math     radian_angle = angle * math.pi / 180     return length * math.cos(radian_angle), length * math.sin(radian_angle) if __name__ == '__main__':     main() If you can explain what your input data looks like, the code can be adapted. Kind regards, Xander ... View more

Hi, Sorry it took so long to say thank you. After working through the tables and the code I finally got it to run properly without any errors. It was amazing. Thank you. I put together this code for the IDW so that someone can select several different fields from the same table and run IDW on the data for that set of points. It's really crude but it works I suppose. Again I want to thank you for all your help. If you have any suggestions on the stuff at the bottom let me know. import arcpy from arcpy import env from arcpy.sa import * arcpy.CheckOutExtension("Spatial") env.workspace = "C:/Users/Anne/Desktop/Python/finalproj" #Set Variables points = "Export_Output.shp" zfield1 = "pH" zfield2 = "Pphate" zfield3 = "K" zfield4 = "CA" zfield5 = "Mg" #First Run outidw1 = Idw(points, zfield1) outidw1.save("C:/Users/Anne/Desktop/Python/finalproj/data/outputs/ab10ph") #Second Run outidw2 = Idw(points, zfield2) outidw2.save("C:/Users/Anne/Desktop/Python/finalproj/data/outputs/ab10pphate") #Third Run outidw3 = Idw(points, zfield3) outidw3.save("C:/Users/Anne/Desktop/Python/finalproj/data/outputs/ab10ll") #Fourth Run outidw4 = Idw(points, zfield4) outidw4.save("C:/Users/Anne/Desktop/Python/finalproj/data/outputs/ab10ca") #Fifth Run outidw5 = Idw(points, zfield5) outidw5.save("C:/Users/Anne/Desktop/Python/finalproj/data/outputs/ab10mg") Hi Anne, Glad to hear you worked it out. I would like to tell you something on using dictionaries and looping through data. In your sample where you calculate the IDW for 5 fields, you could have used a dictionary containing the input field names and output raster names. import arcpy, os arcpy.CheckOutExtension("Spatial") ws_out = "C:/Users/Anne/Desktop/Python/finalproj/data/outputs" env.workspace = ws_out # set it to your output workspace arcpy.env.overwriteOutput = True #Set Variables points = "C:/Users/Anne/Desktop/Python/finalproj/Export_Output.shp" # dictionary with fieldnames and output raster names d = {'pH': 'ab10ph', 'Pphate': 'ab10pphate', 'K': 'ab10ll',      'CA': 'ab10ca', 'Mg': 'ab10mg'} # loop through the dictionary for fieldname, rasname in d.items():     # are you sure you want to accept default cellsize, power and search radius?     outidw = arcpy.sa.Idw(points, fieldname)     outidw.save(rasname) # no path included, so saved to (output) workspace # return the sa license arcpy.CheckInExtension("Spatial") In case you want to use the field name in the output raster name you could a list as well: l = ['pH', 'Pphate', 'K', 'CA', 'Mg'] ras_prefix = 'ab10' for name in l:     outidw = arcpy.sa.Idw(points, name)     outidw.save('{0}{1}'.format(ras_prefix, name)) Another thing to remember is that you are not specifying any cellsize, power and search radius settings. If you have your environment settings setup correctly and your points are also located outside your study area, this might give appropriate results. If your area is cutoff (since it will by default take the extent of the points) you may want to define things like arcpy.env.snapRaster and arcpy.env.extent. Kind regards, Xander ... View more

I second Richard on this. What does the comparison of linguistic signatures has to do with GIS? In the recent past I answered a post similar to yours, here: http://forums.arcgis.com/threads/103785-Hey-I-was-wondering-about-this-function-problem-for-a-while Since all you posts are not related to GIS, you can probably better take your questions to a forum like stackoverflow? ... View more

Hi Jennifer, First of all, I'm glad that things work for you, but somehow there is something that doesn't make sense to me. The IDW interpolation is available in Spatial Analyst, 3D Analyst and Geostatistical Analyst. The IWD in ga is different, since it also requires a layer, so let's take this one out of the equation. If you manually start the tool (in 10.2) and copy the Python snippet after a successful run, you can see the code generated. In case of 3D Analyst it call the tool from; arcpy.Idw_3d In case of Spatial Analyst it does this: arcpy.gp.Idw_sa Apparently it calls it from the old (?) gp module. That's odd... If you go the the Help page on Idw in Spatial Analyst and scroll down to the code sample you see they are using; arcpy.sa.Idw The 3D Analyst IDW code sample in the Help is using the same arcpy.Idw_3d. You are using arcpy.Idw_sa. In Python you can print the help from any function like this: help(arcpy.sa.Idw) It prints the Help on the function. Just have a closer look at what happens: arcpy.Idw_3d Help on function Idw in module arcpy.ddd: Idw(in_point_features=None, z_field=None, out_raster=None, cell_size=None, power=None, search_radius=None, in_barrier_polyline_features=None)     arcpy.ddd.Idw Help on function Idw in module arcpy.ddd: Idw(in_point_features=None, z_field=None, out_raster=None, cell_size=None, power=None, search_radius=None, in_barrier_polyline_features=None) arcpy.sa.Idw Help on function Idw in module arcpy.sa.Functions: Idw(in_point_features, z_field, cell_size='#', power='#', search_radius='#', in_barrier_polyline_features='#') arcpy.gp.Idw_sa Help on function <lambda> in module arcpy.geoprocessing._base: <lambda> lambda *args arcpy.Idw_sa AttributeError: 'module' object has no attribute 'Idw_sa'     So, the 3D IDW does have an output raster parameter. The Idw_sa (inside gp or in directly from arcpy) do not show any help. However from the python snippet obtained from manually executing the Spatial Analyst IDW we know this: arcpy.gp.Idw_sa("myPointFeatureClass","FieldName","myOutputRaster","cellSize","power","some other parameters","#") It does have an output raster parameter. Many tools have different ways to execute them: #this: arcpy.management.AddField # is the same as: arcpy.AddField_management What you should take into account for future development is that within Spatial Analyst there is a trend to use the structure: outputraster = tool(parameters) This creates a temporal raster object which can be saved by calling the save methode on the object; outputraster.save(outPathAndName) So use arcpy.sa.Exp and arcpy.sa.Idw (instead of arcpy.Exp_sa and arcpy.Idw_sa). It would be nice if Jason Scheirer (from the Python development team at Esri) could jump in the conversation and shed his light on this, but I suppose he is too busy with the DevSummit preps which is starting today... Kind regards, Xander ... View more

Hi Jennifer, There are a few suggestion I have when I look at your second snippet of code: use os.path.join(pathname, filename) to create a valid output (requires "import os") use "IDW_{0}".format(COCLayer) to create the name access the Idw tool through arcpy.sa.Idw note that the output is not a parameter of IDW the result of a arcpy.sa tool is normally a raster object (in memory) use the .save methode on the object to save the raster to a location import arcpy, os geoDB = r'C:\Path\to\Folder\With\GeoDB.gdb' COCLayer = 'aLayerName' COCField = 'aFieldName' IDWCellSize = 10 IDWPower = 2 IDWRadius = 100 IDWBarrier = r'C:\Path\to\Folder\With\GeoDB.gdb\myBarriers' # use os.path.join to combine path and file names and use string.format to format names idw_out = os.path.join(geoDB, "IDW_{0}".format(COCLayer)) # Idw (in_point_features, z_field, {cell_size}, {power}, {search_radius}, {in_barrier_polyline_features}) fc_points = os.path.join(geoDB, 'temp_data') myIDWraster = arcpy.sa.Idw(fc_points, COCField, IDWCellSize, IDWPower, IDWRadius, IDWBarrier) myEXPraster = arcpy.sa.Exp(myIDWraster) myFinalRaster = myEXPraster / 1000 # arcpy.Raster(myEXPraster) / 1000 myFinalRaster.save(idw_out) Kind regards, Xander BTW see samples in Help: Exp (Spatial Analyst) IDW (Spatial Analyst) ... View more

Hi Jen, You're close. What is going wrong in your script is that in the definition of the  function you provide your field name, which contains a dot. What you  should provide is the name (any name) for a variable. In my example I  use "building_height" and "store_height". As pre-logic script code you can use this (this includes logic for buildings that have more than 2 stories and being able to define the storey height): def getStories(building_height, store_height):     return "{0} storey".format(int(building_height / store_height) + 1) The function is called like this: getStories( !YourFieldName!, 4) Please not that in this code 4 height units will return "1 storey", but 4.1 height units, will return "2 storey". I don't know if that is entirely correct... Kind regards, Xander ... View more

Hi Mark, I had some problems with the "custom" coordinate system, but I have attached a toolbox with the functionality and also an updated file geodatabase. Unzip the toolbox ("toolshare2.zip") and navigate with the ArcCatalog window to the unzip location. Try the tools inside. Hope they work for you. The SL distance (length) that is returned follows the linear unit of the input featureclass. In your case it is US_Foot. Kind regards, Xander ... View more

Hi Mark, When I was asking for a few streams, I was refering to a featureclass. Maybe you can attach a few streams as featureclass in a zipped fgdb or as a zipped shapefile (although some fieldnames will be renamed when exporting to shapefile). This way I can look at the geometry itself and the spatial reference and test my script. Kind regards, Xander ... View more

Hi Mark, If your data is located somewhere at the border of Pennsylvania and West Virginia, the coordinates are in decimal degrees (not miles). The coordinates seem to correspond to a stream, just next to the Golden Oaks rd (PA 18). Using decimal degrees to calculate the angle and the straight distance in the scrips will not give correct answers. The coordinates should be projected first. I don't know what the reason is for the calculation to trow a 99999 error. If you are willing to share a small part of your data (streams) I can have a look. Kind regards, Xander ... View more