|
DOC
|
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
07-25-2014
01:27 PM
|
3
|
0
|
3984
|
|
POST
|
Hi Todd, Maybe a little Python can help: Look at the code below:
def GetPixelValueForPercentile(dctper, percentile):
"""will return last pixel value
where percentile LE searched percentile."""
for k in sorted(dctper.keys()):
perc = dctper
if perc <= percentile:
pix_val = k
else:
break
return pix_val
import arcpy
# path to raster (has to be integer!)
ras = r"C:\Path\To\Folder\Or\Geodatabase\With\Raster"
# create dictionary with value vs count
flds = ("VALUE", "COUNT")
dct = {row[0]:row[1] for row in arcpy.da.SearchCursor(ras, flds)}
# calculate number of pixels in raster
cnt_sum = sum(dct.values())
dct_per = {}
cnt_i = 0
# loop through dictionary and create new dictionary with val vs percentile
for val in sorted(dct.keys()):
cnt_i += dct[val]
dct_per[val] = cnt_i / cnt_sum
# use dct_per to determine percentiles
perc = 0.05 # (5%)
pixval = GetPixelValueForPercentile(dct_per, perc)
print "Found for percentile {0}, pixel value {1} with real percentile {2}".format(perc, pixval, dct_per[pixval])
perc = 0.95 # (95%)
pixval = GetPixelValueForPercentile(dct_per, perc)
print "Found for percentile {0}, pixel value {1} with real percentile {2}".format(perc, pixval, dct_per[pixval])
It will return the following information:
>> Found for percentile 0.05, pixel value 401 with real percentile 0.049861616298
>> Found for percentile 0.95, pixel value 2350 with real percentile 0.949769065993
The code requires that your raster is of type integer (to make sure it has an attribute table with pixelvalue and count. This is what happens in the code: on line 16 a path to your integer raster is defined (e.g. ras = r'D:\Xander\tmp\ras.gdb\DEMsmallInt') on line 20 the attribute table is used by a da search cursor to create a dictionary (this helps to speed up the proces) on line 23 the sum of the count is calculated (the number of pixel with a value in the raster) on line 29 a dictionary is created that holds the pixel value and corresponding percentile on line 35 and 39 the function (defined at the top) is used to search the pixel value for a percentile value. Hope this is of any use. Kind regards, Xander
... View more
07-25-2014
12:45 PM
|
2
|
47
|
7736
|
|
DOC
|
Some explanation of using python (arcpy) and matplotlib to create profile graphs. Enjoy, Xander
... View more
07-24-2014
03:29 PM
|
5
|
7
|
10601
|
|
POST
|
So, if I understand it correctly your folderstructure does not contain all raster in the same folder (workspace), but in subfolders of the workspace (so the right part if the image below): In that case the arcpy.da.Walk option (available from 10.1 SP 1) as suggested by Ian Murray is the right way to go. The code would look like this (once again, the code is not tested):
import arcpy,os
from arcpy import env
from arcpy.sa import *
#To overwrite output
arcpy.env.overwriteOutput = True
#Set environment settings
ws = "C:/Subhasis/Test/Neshanic_Python"
env.workspace = ws
outws="C:/Subhasis/Test/Neshanic_Python/extract"
threshold = 0.75
#checkout ArcGIS spatial analyst extension license
arcpy.CheckOutExtension("Spatial")
# create a list of all rasters in current workspace
dct_ras = {}
# use arcpy.da.Walk (available from 10.1 SP1) to generate dictionary of rasters
for dirpath, dirnames, filenames in arcpy.da.Walk(ws,
topdown=True,
datatype="Raster"):
for filename in filenames:
ras = os.path.join(dirpath, filename)
dct_ras[ras] = filename
# loop through dictionary and extract path to input ras and rasname
for ras, rasname in dct_ras.items():
flds = ("VALUE", "COUNT")
dct = {row[0]:row[1] for row in arcpy.da.SearchCursor(ras, flds)}
sumcnt = sum(dct.values())
threscnt = sumcnt * threshold
cnt, result = 0, 0
for k in sorted(dct):
cnt += dct
# print " - Value = {0}, Sum = {1}".format(k, cnt)
if cnt < threscnt:
result = int(k)
else:
break
# print "Threshold = {0}".format(result)
attExtract = ExtractByAttributes(ras, "VALUE>={0}".format(result))
outname = os.path.join(outws,"{0}{1}".format(rasname, result))
attExtract.save(outname)
# return SA extension
arcpy.CheckInExtension("Spatial")
In this case I do not write the rasters to a list with simply the raster names, since they reside in different workspaces. The rasters path is used as dictionary key and the raster name is used as value. Do notice that in your structure the raster names can repeat and create output rasters with same name (in case the threshold value is equal too). Be ware of that. Kind regards, Xander
... View more
07-24-2014
02:47 PM
|
0
|
1
|
3782
|
|
DOC
|
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
07-23-2014
08:03 PM
|
40
|
47
|
26765
|
|
BLOG
|
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
07-23-2014
07:22 PM
|
3
|
8
|
1922
|
|
DOC
|
Nice document that will get your started with rasters in arcpy
... View more
07-23-2014
04:00 PM
|
0
|
0
|
1053
|
|
POST
|
Hi Subhasis, I suppose the code you posted did not work, right? The list you created to test the code is actually a list with only one element. The list should have been created as (notice the missing quotes): inraster = ["01367620-r-r","01367700-r-r","01367770-r-r"] Don't worry about the folders, each folder contains an Esri grid and the info folder is a common storage for those grids. The other part that wasn’t optimal was this block:
attExtract = ExtractByAttributes(str(i), "VALUE>=" + str("{0}".format(result)))
outname=os.path.join(outws,str(i), str("{0}".format(result)))
attExtract.save(outname
Apart from the missing bracket at the end, you don’t need to cast the name of the raster to string (it is already string). The strength of using “format” is that you can combine strings and values together in a much easier way. I haven't tested the code below:
import arcpy,os
from arcpy import env
from arcpy.sa import *
#To overwrite output
arcpy.env.overwriteOutput = True
#Set environment settings
env.workspace = "C:/Subhasis/Test/Neshanic_Python"
outws="C:/Subhasis/Test/Neshanic_Python/extract"
#checkout ArcGIS spatial analyst extension license
arcpy.CheckOutExtension("Spatial")
# create a list of all rasters in current workspace
lst_ras = arcpy.ListRasters()
# set local variable (only presented 3 out of 70 for testing)
# inraster = ["01367620-r-r,01367700-r-r,01367770-r-r"]
threshold = 0.75
for ras in lst_ras:
flds = ("VALUE", "COUNT")
dct = {row[0]:row[1] for row in arcpy.da.SearchCursor(i, flds)}
sumcnt = sum(dct.values())
threscnt = sumcnt * threshold
# print "Total number of pixels with a value = {0}".format(sumcnt)
# print "Total number of pixels within threshold = {0}".format(threscnt)
cnt, result = 0, 0
for k in sorted(dct):
cnt += dct
# print " - Value = {0}, Sum = {1}".format(k, cnt)
if cnt < threscnt:
result = int(k)
else:
break
# print "Threshold = {0}".format(result)
attExtract = ExtractByAttributes(ras, "VALUE>={0}".format(result))
outname = os.path.join(outws,"{0}{1}".format(ras, result))
attExtract.save(outname)
arcpy.CheckInExtension("Spatial")
Kind regards, Xander
... View more
07-23-2014
11:00 AM
|
1
|
5
|
3782
|
|
POST
|
Actually, there is a lot more going on in the script that wont work: the nested loop and reusing the same variable names (row, rows, inraster) the where clause will not be filled with the value you're reading the print statement at the end will never show a result since the "row" is printed as text you're mixing the "normal" SearchCursor and da.SearchCursor syntax Furthermore, there is actually no need to create rasters, and therefore you don't need a SA license. Take a look at the code below:
import arcpy
inraster ="C:/Subhasis/Test/Neshanic_Python/01367620-r-r"
threshold = 0.25
flds = ("VALUE", "COUNT")
dct = {row[0]:row[1] for row in arcpy.da.SearchCursor(inraster, flds)}
sumcnt = sum(dct.values())
threscnt = sumcnt * threshold
print "Number of pixels with a value = {0}".format(sumcnt)
print "Maximum number of pixels within threshold = {0}".format(threscnt)
cnt, result = 0, 0
for k in sorted(dct):
cnt += dct
print " - Value = {0}, Sum = {1}".format(k, cnt)
if cnt < threscnt:
result = int(k)
else:
break
print "Threshold = {0}".format(result)
On line 7 a dictionary is created using the da.SeachCursor. A dictionary is an object that hold key (=VALUE), value (=COUNT) pairs. Insert a print statement in the code to print the "dct" object to see how it is filled. The sum of all the counts is calculated on line 8 On line 9 a variable is introduced that holds the maximum number of pixels based on the specified threshold. Line 14 till 21 loop through the sorted dictionary (it is sorted on the pixel values) and calculates the sum of the COUNT of all the pixel values less than or equal to the current pixel value. In case this count is less than the maximum number of pixels calculated on line 9, the result is overwritten with the new best result. Have a look at the output generated with a test raster: Number of pixels with a value = 219996 Maximum number of pixels within threshold = 54999.0 - Value = 0, Sum = 10649 - Value = 1, Sum = 10721 - Value = 2, Sum = 10807 - Value = 3, Sum = 10896 << other values omitted>> - Value = 38, Sum = 53141 - Value = 39, Sum = 54281 - Value = 40, Sum = 55675 Threshold = 39 Kind regards, Xander
... View more
07-18-2014
11:15 PM
|
0
|
7
|
3782
|
|
BLOG
|
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
07-03-2014
07:51 PM
|
3
|
11
|
3418
|
|
POST
|
Hi Will, As far as I can see it is not possible to do this using the interface of ArcGIS. This normally means that you won't be able to do so using Python code. The easiest way would be to add the M and Z values as attribute. This way you can still use your code. Another possibility is to use matplotlib (pyplot), but this requires an additional installation (http://matplotlib.org/downloads.html, which is free). To give you an idea of how the code would look, see the (simplified) sample below: import arcpy
import os
import matplotlib.pyplot as plt
from arcpy import env
infc = "\HaulRd.gdb\Treatment\XS_2009BEBRklnTerr_interpZ_Rte_isect"
flds = ('SHAPE@Z','SHAPE@M')
x = [] # Z
y = [] # M
outFolder = r'C:\path\to\output\folder'
outFile = os.path.join(outFolder, "NameForThisProfile.png")
fig = plt.figure(figsize=(10, 5)) # size in inches
with arcpy.da.SearchCursor(infc, flds) as rows:
for row in rows:
x.append(row[0]) # Z
y.append(row[1]) # M
plt.plot(x,y,'r',linewidth=1.0)
plt.plot(x,y,'ro',alpha=0.5)
plt.xlabel('x label text')
plt.ylabel('y label text')
plt.title('Give the graph a title')
fig.savefig(outFile, dpi=300) Kind regards, Xander
... View more
06-27-2014
11:29 AM
|
0
|
0
|
1015
|
|
POST
|
Hi Mohamad, If you create a Toolbox and want to add a script, you can use the "Weighted Sum" Data Type for the parameter to collect the necessary information. This will return a string with the following structure: rastername1 fieldname1 weight1;rastername2 fieldname2 weight2 In the Help on Weighted Sum you can see in the example code that (like Mark explains) a WSTable can be created: wstbl = WSTable([["snow", "VALUE", 0.25], ["land", "VALUE", 0.25],["soil", "VALUE", 0.5]]) If you execute this python code: import arcpy
from arcpy.sa import *
arcpy.CheckOutExtension("Spatial")
wstbl = WSTable([["snow", "VALUE", 0.25], ["land", "VALUE", 0.25], ["soil", "VALUE", 0.5]])
print wstbl It will return: 'snow' VALUE 0.25; 'land' VALUE 0.25; 'soil' VALUE 0.5 ... which is a similar structure as the parameter provides. Normally to understand what a tool needs to execute, you can perform it manually and from the Results screen Copy the code as Python Snippet. Although the Help uses a syntax like: output_raster = arcpy.sa.WeightedSum(wstbl) The copied Python snippet will use: arcpy.gp.WeightedSum_sa("snow VALUE 0.25;land VALUE 0.25;soil VALUE 0.5",output_raster) I would recommend to use the parameter type "Weighted Sum" and provide that as text to the "arcpy.gp.WeightedSum_sa" tool. import arcpy
arcpy.CheckOutExtension("Spatial")
wstbl = arcpy.GetParameterAsText(0) # list of rasters and weights
out_ras = arcpy.GetParameterAsText(1) # output raster
arcpy.gp.WeightedSum_sa(wstbl,out_ras) Kind regards, Xander
... View more
06-27-2014
08:36 AM
|
0
|
0
|
919
|
|
POST
|
Hi Gaston, In theory it is not very difficult to transform the code from an "arcpy.da" search cursor to the "arcpy" cursor. The problem is that the arcpy.SearchCursor does not seem to support raster fields. The arcpy.da method that works: import arcpy, os
output_path = r'C:\Forum\Raster attribute\output'
fc = r'C:\Forum\Raster attribute\test.gdb\rasatt01'
fld_raster = 'Raster'
fld_name = 'Name'
flds = (fld_raster, fld_name)
with arcpy.da.SearchCursor(fc, flds) as cursor:
for row in cursor:
filename = "{0}.tif".format(row[1])
row[0].save(os.path.join(output_path, filename)) The arcpy method that doesn't work: import arcpy, os
output_path = r'C:\Forum\Raster attribute\output'
fc = r'C:\Forum\Raster attribute\test.gdb\rasatt01'
fld_raster = 'Raster'
fld_name = 'Name'
cursor = arcpy.SearchCursor(fc)
for row in cursor:
filename = "{0}.tif".format(row.getValue(fld_name))
row.getValue(fld_raster).save(os.path.join(output_path, filename))
del cursor, row This will throw an error "AttributeError: 'NoneType' object has no attribute 'save'". You might be interested in a similar thread regarding Blob fields and Raster fields: http://forums.arcgis.com/threads/108402-How-to-load-image-to-Raster-attribute-field-with-python?p=387169&viewfull=1#post387169 Kind regards, Xander
... View more
06-10-2014
05:10 PM
|
1
|
0
|
644
|
|
POST
|
Hi David, Indeed the assumption "firstpnt = part[len(part) + 0]" is wrong. The "part" represents a list of points. Between the square brackets you put the index to the point you're interested in. The list index starts with 0 (for the first point) and ends with the number of points minus 1. Therefor, in case of the last point we used "len(part) - 1", which is the number of points - 1. For the first point you just use index 0. import arcpy
fc = r'c:\yourFGDB.gdb\yourLines'
# update cursor (be careful with this and try it on a copy of your data)
fields = ["SHAPE@","XCentroid","YCentroid"]
with arcpy.da.UpdateCursor(fc,fields) as curs:
for row in curs:
geom = row[0]
x_new = row[1]
y_new = row[2]
array = arcpy.Array()
for part in geom:
firstpnt = part[0]
firstpnt.X = x_new
firstpnt.Y = y_new
part[0] = firstpnt
array.add(part)
curs.updateRow([arcpy.Polyline(array), x_new, y_new])
Hope this explanation helps you to achieve your goal. Kind regards, Xander
... View more
06-06-2014
07:15 AM
|
0
|
0
|
380
|
| Title | Kudos | Posted |
|---|---|---|
| 6 | 12-20-2019 08:41 AM | |
| 1 | 01-21-2020 07:21 AM | |
| 2 | 01-30-2020 12:46 PM | |
| 1 | 05-30-2019 08:24 AM | |
| 1 | 05-29-2019 02:45 PM |
| Online Status |
Offline
|
| Date Last Visited |
11-26-2025
02:43 PM
|