Latest Contributions by XanderBakker

An example to show the power of List Comprehensions. Some additional resources: https://docs.python.org/2/tutorial/datastructures.html#list-comprehensions http://www.pythonforbeginners.com/lists/list-comprehensions-in-python/ ... View more

Some explanation of using python (arcpy) and matplotlib to create profile graphs. Enjoy, Xander ... View more

General snippets Split path and filename filePathName = 'C:/Folder/SubFolder/afile.ext' filePath,fileName=os.path.split(filePathName) # Results in: filePath = 'C:/Folder/SubFolder' fileName = 'afile.ext'            Check if field exists def FieldExist(tbl, fieldname):     """Check if a field exists, return boolean"""     return bool(arcpy.ListFields(tbl, fieldname))          Formatting leading zero's # format a number to 3 digits (with leading zero's) text = "%03d" % (value,)         or following the suggestion of Neil Ayres‌: length = 7 value = 123 text ="{0}".format(value).zfill(length) # or indicating the number of decimals: length = 10 value = 123.4567890 decimals = 2 text = "{0}".format(round(value, decimals)).zfill(length) Also have a look at this website: Python String Format Cookbook | mkaz.com  as was suggested by Anthony Giles in this thread: What is label expression for formatting a number to have thousands separator and decimals? Working with dictionaries # Create a dictionary: myDict = dict() myDict = {}            # Check if key exists if myKey in myDict:     # do something            # Read the value corresponding to a key if myKey in myDict:     myVal = myDict[myKey]            # add a new key, value pair: myDict.update({aKey: aValue})            # update an existing key, value pair if myKey in myDict:     myVal = statsDict[myKey]     # do something with value     myVal += 1     # update key, value pair in dictionary     myDict.update({myKey: myVal}) # better and shorter is: if myKey in myDict:     myDict[myKey] += 1            Dictionaries are very useful to store all kinds of information as Blake T‌ indicated in a comment below: layers = {'myOutputPoints': 'sourcePoints',           'myOutputLines': 'sourceLines',           'myOutputPolys': 'sourcePolys'} for out_lyr, in_lyr in layers.items():     # do something with the layer names Personally, I store my database names and connection strings in dictionaries. To sort a dictionary by value (not by key) use a lambda function: for key, val in sorted(dct.items(), key=lambda x: x[1]):     # do something with key and value pair            To create a dictionary from OID and an arbitrary field: fc = r'C:\path\to\your\featureclass\or\table' fld_oid = arcpy.Describe(fc).OIDfieldname fld_other = 'afieldname' d = {r[0]: r[1] for r in arcpy.da.SearchCursor(fc, (fld_oid, fld_other))}            To conditionally add items to dictionary, use list comprehensions: d = {r[0]: r[1] for r in arcpy.da.SearchCursor(fc, (fld_oid, fld_other)) if r[1] > somevalue}            Working with lists Convert string to list and back myListText = 'A;B;C;1;blah' myList = myListText.split(';') myListText =  ';'.join(myList) # convert it back print myList[0] # prints 'A'         Convert all items in a list to uppercase This construct is called a list comprehension and is more efficient than creating a list with a loop. uppList = [x.upper() for x in lowList]            Compare lists and dictionaries # Compare two lists and get the items in list 1 but not in list 2 listDifference = list(set(myList1) - set(myList2)) # Do the same, but sort the list on he fly listDifference = sorted(list(set(myList1) - set(myList2))) # Compare two lists and get the items that are in both lists listSame = list(set(myList1) & set(myList2)) # Do the same, but sort the list on he fly listSame = sorted(list(set(myList1) & set(myList2)))            Combining lists # suppose you have two lists: lstOne = ['A','B'] lstTwo = ['C','D','E','F'] # and you want to make a list with the items from both lists, using append lstOne.append(lstTwo) # ... will result in lstOne = ['A','B',['C','D','E','F']] # It adds the list as item 3 (nested list)! # you can loop through list and add the items: lstOne = ['A','B'] lstTwo = ['C','D','E','F'] for itemTwo in lstTwo:     lstOne.append(itemTwo) # ... will result in lstOne = ['A','B','C','D','E','F'] # you can also sum the two lists: lstOne = ['A','B'] lstTwo = ['C','D','E','F'] lstThree = lstOne + lstTwo # ... will result in lstThree = ['A','B','C','D','E','F'] # or you can use the 'extend': lstOne = ['A','B'] lstTwo = ['C','D','E','F'] lstOne.extend(lstTwo) # ... will result in lstOne = ['A','B','C','D','E','F']            Note that append and extend do not return any value and change the list itself. Using the + operator will create a new list Enumerate a list The examples below show what you can achieve by using the enumerate function on a list: mylist = ["a","b","c","d"] print list(enumerate(mylist)) # [(0, 'a'), (1, 'b'), (2, 'c'), (3, 'd')] print [(i, j) for i, j in enumerate(mylist)] # [(0, 'a'), (1, 'b'), (2, 'c'), (3, 'd')] print [pair for pair in enumerate(mylist)] # [(0, 'a'), (1, 'b'), (2, 'c'), (3, 'd')] # which basically does the same as: print [(mylist.index(a), a) for a in mylist] # [(0, 'a'), (1, 'b'), (2, 'c'), (3, 'd')]     An example of using the enumerate function to write a data table to Excel (using the xlwt module) was contributed by Neil Ayres (Thanx Neil): tbl = "Table to export" fldOb = arcpy.ListFields(tbl) FieldNames = [f.name for f in fldOb] # set up excel output outXL = tbl + ".xls" wb = xlwt.Workbook() # add first sheet SheetNum = 1 ws = wb.add_sheet("Data_" + str(SheetNum)) # first row in sheet, header line r = 0 [ws.write(r, c, head) for c, head in list(enumerate(FieldNames))] with arcpy.da.SearchCursor(tbl, FieldNames) as cur:     for row in cur:         # if row limit is reached, increment and insert another sheet         if r > 65534:             SheetNum += 1             ws = wb.add_sheet("Data_" + str(SheetNum))             r = 0             [ws.write(r, c, head) for c, head in list(enumerate(FieldNames))]         r += 1           [ws.write(r, c, val) for c, val in list(enumerate(row))] wb.save(outXL) del wb     Using arcpy.da with lists Get a list of unique values from a table import arcpy set_one = set(r[0] for r in arcpy.da.SearchCursor(FC_or_BL, (fieldName)))         Extract a field value from a table using arcpy.da def GetFieldValue(FClassOrTable, field):     return arcpy.da.SearchCursor(FClassOrTable, field).next()[0] def GetMaxFieldValue(FClassOrTable, field):     return sorted(arcpy.da.SearchCursor(FClassOrTable, field), reverse=True)[0][0] def GetMinFieldValue(FClassOrTable, field):     return sorted(arcpy.da.SearchCursor(FClassOrTable, field))[0][0]          Determine the sum of a column - using arcpy.da and Statistics def GetSum(FClassOrTable, fldName):     import arcpy     fields = [fldName]     total = 0     with arcpy.da.SearchCursor(FClassOrTable, fields) as cursor:         for row in cursor:             val = row[0]             total += val     return total          If there are many rows, this may not be the fastest way. In that case, you might consider using Statistics_analysis: import arcpy def GetSum(FClassOrTable, fldName):      stats = r"C:\temp\Default.gdb\stats" # Point this to a temporary table     arcpy.Statistics_analysis(FClassOrTable, stats, [[fldName, "SUM"]])     row = arcpy.SearchCursor(stats).next()     arcpy.Delete_management(stats)           return row.getValue("SUM_" + fldName)          Arcpy.da cursor tuples and lists considerations A da search cursor returns a tuple (row), which is immutable. You can, though, convert the tuple to a list, using tmpList = list(inRow)            If you know the position (index) of the item in the tuple (now a mutable list), you can update it by its index. tmpList[indexPosition] = valueToAssign            The list can be converted back into a tuple outRow = tuple(tmpList)            ... and use the cursor to insert the tuple as a new row: inCur.insertRow(outRow)            Using eval in Python The eval command is very powerful. It allows you to construct a command as text and execute it. Let's have a look at the example below: import arcpy import os dct_conn = {"DEV 10.2": r"C:\Users\xbakker\AppData\Roaming\ESRI\Desktop10.2\ArcCatalog\Dev 10.2.sde",             "TST 10.2": r"C:\Users\xbakker\AppData\Roaming\ESRI\Desktop10.2\ArcCatalog\Tst 10.2.sde",             "ACC 10.2": r"C:\Users\xbakker\AppData\Roaming\ESRI\Desktop10.2\ArcCatalog\Acc 10.2.sde",             "PRO 10.2": r"C:\Users\xbakker\AppData\Roaming\ESRI\Desktop10.2\ArcCatalog\Pro 10.2.sde"} fds = "GEOGENESIS.PROY_INFRAESTRUCTURA" lst_props = ['canVersion', 'datasetType', 'DSID', 'extent', 'isVersioned',              'MExtent', 'spatialReference.name', 'spatialReference.factoryCode', 'ZExtent'] for amb, con in dct_conn.items():     arcpy.env.workspace = con     lst_fc = arcpy.ListFeatureClasses(feature_dataset=fds)     if lst_fc:         fc = lst_fc[0]         desc = arcpy.Describe(fc)         for prop in lst_props:             try:                 print "'{0}'={1}".format(prop, eval("desc.{0}".format(prop)))             except:                 pass     In this case, the code loops through 4 connection files and lists the featureclasses stored in a specific feature dataset. For each featureclass it creates a describe object. To list a number of properties (defined in a list called "lst_props") it constructs a string and executes it using the eval command. The magic happens on line 24: eval("desc.{0}".format(prop))     This will format a string. If the property is "spatialReference.name", the string will be "desc.spatialReference.name". When it is evaluated using eval, it will return the name of the spatial reference. This result is used in another string format function to show: >> 'spatialReference.name'=Some Spatial Reference Name     General advice Always use functions (when practical) If you want to know why a function with code is faster than the same code outside a function, read this topic: http://stackoverflow.com/questions/11241523/why-does-python-code-run-faster-in-a-function Clean up after  yourself using a finally block Python scripts that do not delete layers they create create "memory leaks" that can slow down and break ArcGIS! # temp dataset and layer variables tmp1, lyr1 = [None]*2 try:      tmp1 = arcpy.CreateScratchName("","", "table")      lyr1 = "tmplyr"      lyr1 = arcpy.MakeFeatureLayer(tmp1, lyr1) except:     # error handling finally:     for k in [lyr1, tmp1]:         if k: arcpy.Delete_management(k)          Arcpy and License levels When you create a script and execute it inside the Python Window of ArcMap, the process will use the same license level as defined in the host (ArcMap). In case of working with stand alone scripts, and you have a pool of licenses (arcview-basic / arceditor-standard / arcinfo-advanced), the script may not take the same license level as you specified in the ArcGIS Administrator, but claim the highest level. There is a function in arcpy called SetProduct that would seem to allow you to set the license level, but that only works with the previous argisscripting (9.x) and not with arcpy. Once you import arcpy the license level is set and cannot be changed. So if you are using ArcMap with a Basic license, your stand alone script might be claiming an additional Standard or Advanced license on the same machine! To force the use of a license level you should import the desired license level before importing arcpy. Valid license levels are: 'arcview', 'arceditor' and 'arcinfo' (and 'engine', 'enginegeodb', 'arcserver') To force a Basic (arcview) license, do this: import arcview import arcpy # print the license level print "Initial ProductInfo: {0}".format(arcpy.ProductInfo())    Please be aware of this, since you might be claiming two licenses on one machine! Beware, at 10.3 ArcGIS will consume the highest level available and the import statement of a lower level will not force a lower license level! See thread: import arceditor does not set proper license level Working with domains Below a example of obtaining the domain description when reading attributes of a feature. There is no error handling in this snippet, so be aware... import arcpy, os def main():     # input data     fgdb = r"D:\Xander\Genesis\Procedimientos\CargarDatos\actualizacion.gdb"     fc_name = "PROY_Proyecto_Gen"     fld_name = "FASE_PROY"     # get featureclass     fc = os.path.join(fgdb, fc_name)     # get field and the domain name     fld = arcpy.ListFields(fc, wild_card=fld_name)[0]     dom_name = fld.domain # returns domain name     # list the domains and get the domain that we need     doms = arcpy.da.ListDomains(fgdb)     dom = GetDomainOnName(dom_name, doms)     # let's loop through the data and return the code and the description     with arcpy.da.SearchCursor(fc, ("OID@", fld_name)) as curs:         for row in curs:             oid = row[0]             code = row[1]             description = GetDescription(code, dom)             print "oid {0}, code={1}, description={2}".format(oid, code, description) def GetDescription(code, dom):     if dom.domainType == 'CodedValue':         coded_values = dom.codedValues         if code in coded_values:             return coded_values         else:             return "code not found in domain" def GetDomainOnName(name, doms):     dom_found = None     for dom in doms:         if dom.name == name:             dom_found = dom             break     return dom_found if __name__ == '__main__':     main()   in my case this yields something like this: oid 867, code=3, description=Design oid 868, code=3, description=Design oid 869, code=9, description=Operation oid 870, code=9, description=Operation oid 871, code=9, description=Operation   You will find an even better example at the ArcPy Café written by ArcGIS Team Python: Get coded-value descriptions when accessing data with cursors | ArcPy Café Please note that when exporting a featureclass of table, normally the domain code (the actual content of the field) will be exported. To include the domain description in the output, there is a setting to achieve this: Menu "Geoprocessing", select "Environments...". Next Expand "Fields", and switch on the option "Transfer field domain descriptions". In code this is done as follows: import arcpy from arcpy import env # settings outLocation = r"path you output folder" inFeatures = r"path to input featureclass with code value domains" # configurar setting env.transferDomains = True # or "TRANSFER_DOMAINS" is also allowed # convert to shapefile outName = "my_output_name.shp" arcpy.conversion.FeatureClassToFeatureClass(inFeatures, outLocation, outName) ... View more

When you create a dictionary in Python the order will not necessary be the order in which you create it. There is a solution for this and that is the use of an OrderedDict (can be found in the collections module). Take the following example; you have a dictionary you want to use to classify an angle (0-360) into a cardinal direction. You define the dictionary, but if you print the dictionary directly afterwards you will notice that the order is not necessary the same as defined: dctCard = {22.5:"E", 67.5: "NE", 112.5: "N",            157.5: "NW", 202.5: "W", 247.5: "SW",            292.5: "S", 337.5: "SE", 360: "E"} print dctCard results in: >> {22.5: 'E', 67.5: 'NE', 247.5: 'SW', 202.5: 'W', 292.5: 'S', 112.5: 'N', 360: 'E', 337.5: 'SE', 157.5: 'NW'} So if you would use code like below to obtain the cardinal direction based on an angle, this will give erroneous results: dctCard = {22.5:"E", 67.5: "NE", 112.5: "N",            157.5: "NW", 202.5: "W", 247.5: "SW",            292.5: "S", 337.5: "SE", 360: "E"} myangle = 75 for angle, cardinal in dctCard.items():     if myangle < angle:         result = cardinal         break print "Angle {0} gives {1}".format(myangle, result) could result in: >> Angle 75 gives SW Ordered dictionary It is possible to solve this by using collections.OrderedDict(). Please be aware that the code below will still result in a randomly sorted dictionary: from collections import OrderedDict dctCard = OrderedDict({22.5:"E", 67.5: "NE", 112.5: "N",            157.5: "NW", 202.5: "W", 247.5: "SW",            292.5: "S", 337.5: "SE", 360: "E"}) The reason is, that first the dictionary is created in a random order and than that random order is stored as "ordered". There are two ways of solving this. Either you loop through a sorted list of keys (the angle values): dctCard = {22.5:"E", 67.5: "NE", 112.5: "N",            157.5: "NW", 202.5: "W", 247.5: "SW",            292.5: "S", 337.5: "SE", 360: "E"} myangle = 75 for angle in sorted(dctCard.keys()):     if myangle < angle:         result = dctCard[angle]         break Or you create a list of key value pairs and convert that to OrderedDict from collections import OrderedDict lstCardinal = [(22.5,"E"), (67.5, "NE"), (112.5, "N"),                (157.5, "NW"), (202.5, "W"), (247.5, "SW"),                (292.5, "S"), (337.5, "SE"), (360, "E")] dctCard2 = OrderedDict(lstCardinal) myangle = 75 for angle, cardinal in dctCard2.items():     if myangle < angle:         result = cardinal         break ... View more

Today I stumbled upon a tweet showing how children in Japan are using line intersections to calculate the multiplication of two numbers. I was amazed by the beauty of the method which can be extended to larger numbers too.  Just take a look at the example below: Source: http://www.pinpopular.com/japanese-multiplication-trick-pinterest-education/ It is basically counting intersections of lines. Wait ... intersections of lines ... sounds to me that arcpy could do that too. So, why not write some arcpy / python code to multiply 2 numbers? I know, what you're thinking ... why? Because it's fun and completely useless! So this is what I came up with: # Japanese visual multiplication with arcpy import arcpy import os ws = r"C:\Forum\JapaneseMultiply\fgdb\test.gdb" arcpy.env.workspace = ws arcpy.env.overwriteOutput = True sr = arcpy.SpatialReference(3857) arcpy.env.outputCoordinateSystem = sr # values to multiply: value1 = 47 # (1 -99) value2 = 63 # (1 -99) fc_points = "points" fc_lines_h = "lines_h" fc_lines_v = "lines_v" fc_polygones = "polygons" fc_result = "result" size = 10 fld_factor = "factor" val_h1 = value1 / 10 val_h2 = value1 % 10 val_v1 = value2 / 10 val_v2 = value2 % 10 # create horizontal lines arcpy.CreateFeatureclass_management(ws, fc_lines_h, "POLYLINE", "#", "DISABLED", "DISABLED", sr) with arcpy.da.InsertCursor(fc_lines_h, ("SHAPE@")) as rows:     for y in range(1, val_h1 + 1):         x1 = -1 * size         x2 = size         array = arcpy.Array([arcpy.Point(x1, y),                      arcpy.Point(x2, y)])         polyline = arcpy.Polyline(array)         rows.insertRow([polyline])     for y in range(1, val_h2 + 1):         x1 = -1 * size         x2 = size         array = arcpy.Array([arcpy.Point(x1, -1 * y),                      arcpy.Point(x2, -1 * y)])         polyline = arcpy.Polyline(array)         rows.insertRow([polyline]) # create vertical lines arcpy.CreateFeatureclass_management(ws, fc_lines_v, "POLYLINE", "#", "DISABLED", "DISABLED", sr) with arcpy.da.InsertCursor(fc_lines_v, ("SHAPE@")) as rows:     for x in range(1, val_v1 + 1):         y1 = -1 * size         y2 = size         array = arcpy.Array([arcpy.Point(-1 * x, y1),                      arcpy.Point(-1 * x, y2)])         polyline = arcpy.Polyline(array)         rows.insertRow([polyline])     for x in range(1, val_v2 + 1):         y1 = -1 * size         y2 = size         array = arcpy.Array([arcpy.Point(x, y1),                      arcpy.Point(x, y2)])         polyline = arcpy.Polyline(array)         rows.insertRow([polyline]) # Create polygon cuadrants arcpy.CreateFeatureclass_management(ws, fc_polygones, "POLYGON", "#", "DISABLED", "DISABLED", sr) arcpy.AddField_management(fc_polygones, fld_factor, "DOUBLE") with arcpy.da.InsertCursor(fc_polygones, ("SHAPE@", fld_factor)) as rows:     # left top factor 100     factor = 100     array = arcpy.Array([arcpy.Point(-1 * size, 0),                          arcpy.Point(-1 * size, size),                          arcpy.Point(0, size),                          arcpy.Point(0, 0),                          arcpy.Point(-1 * size, 0)])     polygon = arcpy.Polygon(array)     rows.insertRow([polygon, factor])     # right top factor 10     factor = 10     array = arcpy.Array([arcpy.Point(0, 0),                          arcpy.Point(0, size),                          arcpy.Point(size, size),                          arcpy.Point(size, 0),                          arcpy.Point(0, 0)])     polygon = arcpy.Polygon(array)     rows.insertRow([polygon, factor])     # left bottom factor 10     factor = 10     array = arcpy.Array([arcpy.Point(-1 * size, -1 * size),                          arcpy.Point(-1 * size, 0),                          arcpy.Point(0, 0),                          arcpy.Point(0, -1 * size),                          arcpy.Point(-1 * size, -1 * size)])     polygon = arcpy.Polygon(array)     rows.insertRow([polygon, factor])     # right bottom factor 1     factor = 1     array = arcpy.Array([arcpy.Point(0, 0),                          arcpy.Point(0, -1 * size),                          arcpy.Point(size, -1 * size),                          arcpy.Point(size, 0),                          arcpy.Point(0, 0)])     polygon = arcpy.Polygon(array)     rows.insertRow([polygon, factor]) # start analysis: intersect horizontal and vertical lines arcpy.Intersect_analysis([fc_lines_h, fc_lines_v], fc_points, "ALL", "", "POINT") # Next intersect the polygon with the points (intersection of lines) arcpy.Intersect_analysis([fc_polygones, fc_points], fc_result) # now use Search cursor to get result res = 0 with arcpy.da.SearchCursor(fc_result, (fld_factor)) as rows:     for row in rows:         factor = row[0]         res += factor print "{0} * {1} = {2}".format(value1, value2, res) print "verify: {0} * {1} = {2}".format(value1, value2, value1 * value2) In the example above the following result is printed: 47 * 63 = 2961.0 verify: 47 * 63 = 2961 Displaying the featureclasses gives the following result: I hope my next blog post will be of more use. Xander ... View more