<html xmlns:v="urn:schemas-microsoft-com:vml" xmlns:o="urn:schemas-microsoft-com:office:office" xmlns:w="urn:schemas-microsoft-com:office:word" xmlns:m="http://schemas.microsoft.com/office/2004/12/omml" xmlns="http://www.w3.org/TR/REC-html40">
<head>
<meta http-equiv="Content-Type" content="text/html; charset=us-ascii">
<meta name="Generator" content="Microsoft Word 14 (filtered medium)">
<style><!--
/* Font Definitions */
@font-face
{font-family:Calibri;
panose-1:2 15 5 2 2 2 4 3 2 4;}
/* Style Definitions */
p.MsoNormal, li.MsoNormal, div.MsoNormal
{margin:0in;
margin-bottom:.0001pt;
font-size:11.0pt;
font-family:"Calibri","sans-serif";}
a:link, span.MsoHyperlink
{mso-style-priority:99;
color:blue;
text-decoration:underline;}
a:visited, span.MsoHyperlinkFollowed
{mso-style-priority:99;
color:purple;
text-decoration:underline;}
p.MsoListParagraph, li.MsoListParagraph, div.MsoListParagraph
{mso-style-priority:34;
margin-top:0in;
margin-right:0in;
margin-bottom:0in;
margin-left:.5in;
margin-bottom:.0001pt;
font-size:11.0pt;
font-family:"Calibri","sans-serif";}
span.EmailStyle18
{mso-style-type:personal;
font-family:"Calibri","sans-serif";
color:windowtext;}
span.EmailStyle19
{mso-style-type:personal-reply;
font-family:"Calibri","sans-serif";
color:#1F497D;}
.MsoChpDefault
{mso-style-type:export-only;
font-size:10.0pt;}
@page WordSection1
{size:8.5in 11.0in;
margin:1.0in 1.0in 1.0in 1.0in;}
div.WordSection1
{page:WordSection1;}
/* List Definitions */
@list l0
{mso-list-id:1272585816;
mso-list-type:hybrid;
mso-list-template-ids:-432119446 67698705 67698713 67698715 67698703 67698713 67698715 67698703 67698713 67698715;}
@list l0:level1
{mso-level-text:"%1\)";
mso-level-tab-stop:none;
mso-level-number-position:left;
text-indent:-.25in;}
@list l0:level2
{mso-level-number-format:alpha-lower;
mso-level-tab-stop:none;
mso-level-number-position:left;
text-indent:-.25in;}
@list l0:level3
{mso-level-number-format:roman-lower;
mso-level-tab-stop:none;
mso-level-number-position:right;
text-indent:-9.0pt;}
@list l0:level4
{mso-level-tab-stop:none;
mso-level-number-position:left;
text-indent:-.25in;}
@list l0:level5
{mso-level-number-format:alpha-lower;
mso-level-tab-stop:none;
mso-level-number-position:left;
text-indent:-.25in;}
@list l0:level6
{mso-level-number-format:roman-lower;
mso-level-tab-stop:none;
mso-level-number-position:right;
text-indent:-9.0pt;}
@list l0:level7
{mso-level-tab-stop:none;
mso-level-number-position:left;
text-indent:-.25in;}
@list l0:level8
{mso-level-number-format:alpha-lower;
mso-level-tab-stop:none;
mso-level-number-position:left;
text-indent:-.25in;}
@list l0:level9
{mso-level-number-format:roman-lower;
mso-level-tab-stop:none;
mso-level-number-position:right;
text-indent:-9.0pt;}
@list l1
{mso-list-id:1930235544;
mso-list-type:hybrid;
mso-list-template-ids:-1734300642 67698705 67698713 67698715 67698703 67698713 67698715 67698703 67698713 67698715;}
@list l1:level1
{mso-level-text:"%1\)";
mso-level-tab-stop:none;
mso-level-number-position:left;
text-indent:-.25in;}
@list l1:level2
{mso-level-number-format:alpha-lower;
mso-level-tab-stop:none;
mso-level-number-position:left;
text-indent:-.25in;}
@list l1:level3
{mso-level-number-format:roman-lower;
mso-level-tab-stop:none;
mso-level-number-position:right;
text-indent:-9.0pt;}
@list l1:level4
{mso-level-tab-stop:none;
mso-level-number-position:left;
text-indent:-.25in;}
@list l1:level5
{mso-level-number-format:alpha-lower;
mso-level-tab-stop:none;
mso-level-number-position:left;
text-indent:-.25in;}
@list l1:level6
{mso-level-number-format:roman-lower;
mso-level-tab-stop:none;
mso-level-number-position:right;
text-indent:-9.0pt;}
@list l1:level7
{mso-level-tab-stop:none;
mso-level-number-position:left;
text-indent:-.25in;}
@list l1:level8
{mso-level-number-format:alpha-lower;
mso-level-tab-stop:none;
mso-level-number-position:left;
text-indent:-.25in;}
@list l1:level9
{mso-level-number-format:roman-lower;
mso-level-tab-stop:none;
mso-level-number-position:right;
text-indent:-9.0pt;}
ol
{margin-bottom:0in;}
ul
{margin-bottom:0in;}
--></style><!--[if gte mso 9]><xml>
<o:shapedefaults v:ext="edit" spidmax="1026" />
</xml><![endif]--><!--[if gte mso 9]><xml>
<o:shapelayout v:ext="edit">
<o:idmap v:ext="edit" data="1" />
</o:shapelayout></xml><![endif]-->
</head>
<body lang="EN-US" link="blue" vlink="purple">
<div class="WordSection1">
<p class="MsoNormal">Hi numpy list!<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal">I am trying to do some image processing on a number of images, 72 to be specific. I am seeing the python memory usage continually increase. I will say this is being done also using arcpy so that COULD be an issue, but want to rule out scipy
first since I am only using arcpy to grab a set of polygons and convert them to a binary image, and to save images to disk. I pass the raster to my numpy/scipy processing. The code is organized as:<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoListParagraph" style="text-indent:-.25in;mso-list:l0 level1 lfo2"><![if !supportLists]><span style="mso-list:Ignore">1)<span style="font:7.0pt "Times New Roman"">
</span></span><![endif]>Main loop iterating over features<o:p></o:p></p>
<p class="MsoListParagraph" style="text-indent:-.25in;mso-list:l0 level1 lfo2"><![if !supportLists]><span style="mso-list:Ignore">2)<span style="font:7.0pt "Times New Roman"">
</span></span><![endif]>Function which does image processing using scipy<o:p></o:p></p>
<p class="MsoListParagraph" style="text-indent:-.25in;mso-list:l0 level1 lfo2"><![if !supportLists]><span style="mso-list:Ignore">3)<span style="font:7.0pt "Times New Roman"">
</span></span><![endif]>Class which implements some of the image processing features<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal">I placed 2 sleep calls in the code.. one at the end of the image processing function and one at the end of each iteration in the main loop. I am seeing a 20MB memory increase each iteration while looking at my processes in Task Manager:<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal">Iteration 1: 204,312K and 209,908K, respectively<o:p></o:p></p>
<p class="MsoNormal">Iteration 2: 233,728K and 230,192K, respectively<o:p></o:p></p>
<p class="MsoNormal">Iteration 3: 255,676K and 250,600K, respectively<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal">Any ideas? <o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal">Thanks so much!<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal">Joe<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoListParagraph" style="text-indent:-.25in;mso-list:l1 level1 lfo4"><![if !supportLists]><span style="color:#00B0F0"><span style="mso-list:Ignore">1)<span style="font:7.0pt "Times New Roman"">
</span></span></span><![endif]><span style="color:#00B0F0">Main loop<o:p></o:p></span></p>
<p class="MsoNormal">#for each feature class, convert to raster by "Value" and reclass to 1 and 0<o:p></o:p></p>
<p class="MsoNormal">temp_raster = os.path.join(scratchWorkspace,"temp_1.tif")<o:p></o:p></p>
<p class="MsoNormal">temp_raster_elev = os.path.join(scratchWorkspace,"temp_elev.tif")<o:p></o:p></p>
<p class="MsoNormal">temp_reclass = os.path.join(scratchWorkspace,"temp_reclass.tif")<o:p></o:p></p>
<p class="MsoNormal">remap = arcpy.sa.RemapValue([[1,1], ["NoData",0]])<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal">for i,fc in enumerate(potWaterFeatNames):<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal"> #try to delete the temp rasters in case they exist<o:p></o:p></p>
<p class="MsoNormal"> try:<o:p></o:p></p>
<p class="MsoNormal"> #delete temp rasters<o:p></o:p></p>
<p class="MsoNormal"> arcpy.Delete_management(temp_raster)<o:p></o:p></p>
<p class="MsoNormal"> arcpy.Delete_management(temp_raster_elev)<o:p></o:p></p>
<p class="MsoNormal"> arcpy.Delete_management(temp_reclass)<o:p></o:p></p>
<p class="MsoNormal"> except e:<o:p></o:p></p>
<p class="MsoNormal"> arcpy.AddMessage(e.message())<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal"> #clear in memory workspace"<o:p></o:p></p>
<p class="MsoNormal"> arcpy.Delete_management("in_memory/temp_table")<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal"> #print "on feature {} of {}".format(i+1, len(potWaterFeatNames))<o:p></o:p></p>
<p class="MsoNormal"> arcpy.AddMessage("on feature {} of {}".format(i+1, len(potWaterFeatNames)))<o:p></o:p></p>
<p class="MsoNormal"> arcpy.AddMessage("fc = {}".format(fc))<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal"> arcpy.PolygonToRaster_conversion(fc, "Value", temp_raster, "", "", 1)<o:p></o:p></p>
<p class="MsoNormal"> outReclass = arcpy.sa.Reclassify(temp_raster,"Value",remap)<o:p></o:p></p>
<p class="MsoNormal"> outReclass.save(temp_reclass)<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal"> del outReclass<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal"> #call function to process connected components<o:p></o:p></p>
<p class="MsoNormal"> try:<o:p></o:p></p>
<p class="MsoNormal"> proc_im = processBinaryImage(temp_reclass, scratchWorkspace, sm_pixels)<o:p></o:p></p>
<p class="MsoNormal"> except Exception as e:<o:p></o:p></p>
<p class="MsoNormal"> print "FAILED! scipy image processing"<o:p></o:p></p>
<p class="MsoNormal"> arcpy.AddMessage(e.message())<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal"> #convert processed raster to polygons<o:p></o:p></p>
<p class="MsoNormal"> try:<o:p></o:p></p>
<p class="MsoNormal"> out_fc = os.path.join(outWorkspace,fc + "_cleaned")<o:p></o:p></p>
<p class="MsoNormal"> arcpy.RasterToPolygon_conversion(proc_im, out_fc)<o:p></o:p></p>
<p class="MsoNormal"> arcpy.Delete_management(proc_im)<o:p></o:p></p>
<p class="MsoNormal"> except Exception as e:<o:p></o:p></p>
<p class="MsoNormal"> print "FAILED! converting cleaned binary raster to polygon"<o:p></o:p></p>
<p class="MsoNormal"> arcpy.AddMessage(e.message())<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal"> #delete zero value polyons, gridcode = 0<o:p></o:p></p>
<p class="MsoNormal"> try:<o:p></o:p></p>
<p class="MsoNormal"> uc = arcpy.UpdateCursor(out_fc, "gridcode=0",arcpy.Describe(out_fc).spatialReference)<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal"> #loop through 0 rows and delete<o:p></o:p></p>
<p class="MsoNormal"> for row in uc:<o:p></o:p></p>
<p class="MsoNormal"> uc.deleteRow(row)<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal"> del uc<o:p></o:p></p>
<p class="MsoNormal"> except Exception as e:<o:p></o:p></p>
<p class="MsoNormal"> print "FAILED! deleting rows from table"<o:p></o:p></p>
<p class="MsoNormal"> arcpy.AddMessage(e.message())<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal"> #check that number of polygons with gridcode = 1 is greater than 0<o:p></o:p></p>
<p class="MsoNormal"> count = arcpy.GetCount_management(out_fc)<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal"> #add elevation field back in<o:p></o:p></p>
<p class="MsoNormal"> if (count > 0):<o:p></o:p></p>
<p class="MsoNormal"> arcpy.PolygonToRaster_conversion(fc, "z_val", temp_raster_elev, "", "", 1)<o:p></o:p></p>
<p class="MsoNormal"> arcpy.sa.ZonalStatisticsAsTable(out_fc, "Id", temp_raster_elev, "in_memory/temp_table","#","MEAN")<o:p></o:p></p>
<p class="MsoNormal"> arcpy.JoinField_management(out_fc, "Id", "in_memory/temp_table", "Id", ["MEAN"])<o:p></o:p></p>
<p class="MsoNormal"> else:<o:p></o:p></p>
<p class="MsoNormal"> arcpy.Delete_management(out_fc)<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal"> #delete temp rasters<o:p></o:p></p>
<p class="MsoNormal"> arcpy.Delete_management(temp_raster)<o:p></o:p></p>
<p class="MsoNormal"> arcpy.Delete_management(temp_raster_elev)<o:p></o:p></p>
<p class="MsoNormal"> arcpy.Delete_management(temp_reclass)<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal"> #python garbage collection<o:p></o:p></p>
<p class="MsoNormal"> #collected = gc.collect()<o:p></o:p></p>
<p class="MsoNormal"> #print "Garbage collector: collected %d objects." % (collected)<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal"> print "sleeping for 10 seconds"<o:p></o:p></p>
<p class="MsoNormal"> time.sleep(10)<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoListParagraph" style="text-indent:-.25in;mso-list:l1 level1 lfo4"><![if !supportLists]><span style="mso-list:Ignore">2)<span style="font:7.0pt "Times New Roman"">
</span></span><![endif]><span style="color:#00B0F0">Image processing function</span><o:p></o:p></p>
<p class="MsoNormal">def processBinaryImage(imageName, save_dir, sm_pixels):<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal"> fname = os.path.basename(imageName)<o:p></o:p></p>
<p class="MsoNormal"> imRaster = arcpy.Raster(imageName)<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal"> #Grab AGS info from image<o:p></o:p></p>
<p class="MsoNormal"> #use describe module to grab info<o:p></o:p></p>
<p class="MsoNormal"> descData = arcpy.Describe(imRaster)<o:p></o:p></p>
<p class="MsoNormal"> cellSize = descData.meanCellHeight<o:p></o:p></p>
<p class="MsoNormal"> extent = descData.Extent<o:p></o:p></p>
<p class="MsoNormal"> spatialReference = descData.spatialReference<o:p></o:p></p>
<p class="MsoNormal"> pnt = arcpy.Point(extent.XMin, extent.YMin)<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal"> del imRaster<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal"> converter = ConvertRasterNumpy();<o:p></o:p></p>
<p class="MsoNormal"> imArray = converter.rasterToNumpyArray(imageName)<o:p></o:p></p>
<p class="MsoNormal"> imMin = np.min(imArray)<o:p></o:p></p>
<p class="MsoNormal"> imMax = np.max(imArray)<o:p></o:p></p>
<p class="MsoNormal"> print imMin, imMax<o:p></o:p></p>
<p class="MsoNormal"> arcpy.AddMessage("image min: " + str(imMin))<o:p></o:p></p>
<p class="MsoNormal"> arcpy.AddMessage("image max: " + str(imMax));<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal"> <o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal"> #other flags<o:p></o:p></p>
<p class="MsoNormal"> show_image = False #verbose<o:p></o:p></p>
<p class="MsoNormal"> save_flag = False #under verbose<o:p></o:p></p>
<p class="MsoNormal"> gray = False #threshold but keep gray levels<o:p></o:p></p>
<p class="MsoNormal"> make_binary = False #flag if grayscale is true, to binarize final image<o:p></o:p></p>
<p class="MsoNormal"> save_AGS_rasters = True<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal"> #create filter object<o:p></o:p></p>
<p class="MsoNormal"> filter = thresholdMedianFilter()<o:p></o:p></p>
<p class="MsoNormal"> lowValue = 0<o:p></o:p></p>
<p class="MsoNormal"> tImage = filter.thresholdImage(imArray, lowValue, gray)<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal"> #median filter image<o:p></o:p></p>
<p class="MsoNormal"> filtImage1 = filter.medianFiltImage(tImage,1)<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal"> #create structuring element for morphological operations<o:p></o:p></p>
<p class="MsoNormal"> sElem = np.array([[0., 1., 0.],<o:p></o:p></p>
<p class="MsoNormal"> [1., 1., 1.],<o:p></o:p></p>
<p class="MsoNormal"> [0., 1., 0.]])<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal"> #open filtered image<o:p></o:p></p>
<p class="MsoNormal"> gray = False<o:p></o:p></p>
<p class="MsoNormal"> mImage = filter.mOpenImage(filtImage1, sElem, gray)<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal"> #set list to store info change<o:p></o:p></p>
<p class="MsoNormal"> num_its = 100<o:p></o:p></p>
<p class="MsoNormal"> the_change = np.zeros(num_its)<o:p></o:p></p>
<p class="MsoNormal"> for it in range(100):<o:p></o:p></p>
<p class="MsoNormal"> prev = mImage<o:p></o:p></p>
<p class="MsoNormal"> filtImage = filter.medianFiltImage(mImage,1)<o:p></o:p></p>
<p class="MsoNormal"> mImage = filter.mOpenImage(filtImage, sElem, gray)<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal"> #calculate difference<o:p></o:p></p>
<p class="MsoNormal"> (m_info, z_info) = filter.calcInfoChange(prev, mImage)<o:p></o:p></p>
<p class="MsoNormal"> the_change[it] = z_info<o:p></o:p></p>
<p class="MsoNormal"> del filtImage<o:p></o:p></p>
<p class="MsoNormal"> del prev<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal"> #if the change is less than 5% of the initial change, exit<o:p></o:p></p>
<p class="MsoNormal"> cur_perc = the_change[it]/the_change[0]*100<o:p></o:p></p>
<p class="MsoNormal"> if cur_perc < 5.0:<o:p></o:p></p>
<p class="MsoNormal"> print "exiting filter on iteration " + str(it)<o:p></o:p></p>
<p class="MsoNormal"> print "change is less than 5% (this iteration: " + str(cur_perc) + ")"<o:p></o:p></p>
<p class="MsoNormal"> break<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal"> ############################################################################<o:p></o:p></p>
<p class="MsoNormal"> # now we have a binary mask. Let's find and label the connected components #<o:p></o:p></p>
<p class="MsoNormal"> ############################################################################<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal"> #clear some space<o:p></o:p></p>
<p class="MsoNormal"> del tImage<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal"> del filtImage1<o:p></o:p></p>
<p class="MsoNormal"> del m_info<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal"> label_im_init, nb_labels = ndimage.label(mImage)<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal"> #Compute size, mean_value, etc. of each region:<o:p></o:p></p>
<p class="MsoNormal"> label_im = label_im_init.copy()<o:p></o:p></p>
<p class="MsoNormal"> del label_im_init<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal"> sizes = ndimage.sum(mImage, label_im, range(nb_labels + 1))<o:p></o:p></p>
<p class="MsoNormal"> mean_vals = ndimage.sum(imArray, label_im, range(1, nb_labels + 1))<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal"> #clean up small components<o:p></o:p></p>
<p class="MsoNormal"> mask_size = sizes < sm_pixels<o:p></o:p></p>
<p class="MsoNormal"> remove_pixel = mask_size[label_im]<o:p></o:p></p>
<p class="MsoNormal"> label_im[remove_pixel] = 0;<o:p></o:p></p>
<p class="MsoNormal"> labels = np.unique(label_im)<o:p></o:p></p>
<p class="MsoNormal"> label_im = np.searchsorted(labels, label_im)<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal"> #make label image to a binary image, and convert to arcgis raster<o:p></o:p></p>
<p class="MsoNormal"> label_im[label_im > 0] = 1<o:p></o:p></p>
<p class="MsoNormal"> label_im = np.array(label_im, dtype = 'float32')<o:p></o:p></p>
<p class="MsoNormal"> print label_im.dtype<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal"> saveit = False<o:p></o:p></p>
<p class="MsoNormal"> if ~saveit:<o:p></o:p></p>
<p class="MsoNormal"> outRaster = save_dir + "\\" + fname[:-4] + "_CC_" + str(sm_pixels) + ".tif"<o:p></o:p></p>
<p class="MsoNormal"> temp = arcpy.NumPyArrayToRaster(label_im,pnt,cellSize,cellSize)<o:p></o:p></p>
<p class="MsoNormal"> arcpy.DefineProjection_management(temp, spatialReference)<o:p></o:p></p>
<p class="MsoNormal"> arcpy.CopyRaster_management(temp, outRaster, "DEFAULTS","0","","","","8_BIT_UNSIGNED")<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal"> #clear more space<o:p></o:p></p>
<p class="MsoNormal"> del mImage<o:p></o:p></p>
<p class="MsoNormal"> del nb_labels<o:p></o:p></p>
<p class="MsoNormal"> del sizes<o:p></o:p></p>
<p class="MsoNormal"> del mean_vals<o:p></o:p></p>
<p class="MsoNormal"> del mask_size<o:p></o:p></p>
<p class="MsoNormal"> del remove_pixel<o:p></o:p></p>
<p class="MsoNormal"> del label_im<o:p></o:p></p>
<p class="MsoNormal"> del labels<o:p></o:p></p>
<p class="MsoNormal"> del temp<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal"> del the_change<o:p></o:p></p>
<p class="MsoNormal"> del sElem<o:p></o:p></p>
<p class="MsoNormal"> del filter<o:p></o:p></p>
<p class="MsoNormal"> del imArray<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal"> print 'sleeping'<o:p></o:p></p>
<p class="MsoNormal"> time.sleep(20)<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal"> return outRaster<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoListParagraph" style="text-indent:-.25in;mso-list:l1 level1 lfo4"><![if !supportLists]><span style="color:#00B0F0"><span style="mso-list:Ignore">3)<span style="font:7.0pt "Times New Roman"">
</span></span></span><![endif]><span style="color:#00B0F0">Image processing class<o:p></o:p></span></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal">class thresholdMedianFilter():<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal"> def thresholdImage(self, imArray, thresh, binary = True):<o:p></o:p></p>
<p class="MsoNormal"> """ threshold the image. values equal or below thresh are 0, above are 1"""<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal"> tImage = imArray<o:p></o:p></p>
<p class="MsoNormal"> tImage[tImage <= thresh] = 0<o:p></o:p></p>
<p class="MsoNormal"> if binary:<o:p></o:p></p>
<p class="MsoNormal"> tImage[tImage > thresh] = 1<o:p></o:p></p>
<p class="MsoNormal"> return tImage<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal"> def medianFiltImage(self,imArray,n=1):<o:p></o:p></p>
<p class="MsoNormal"> """median filter the image n amount of times. a single time is the default"""<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal"> for n in range(n):<o:p></o:p></p>
<p class="MsoNormal"> prev = imArray<o:p></o:p></p>
<p class="MsoNormal"> imArray = signal.medfilt2d(prev)<o:p></o:p></p>
<p class="MsoNormal"> del prev<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal"> return imArray<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal"> def mOpenImage(self,imArray, sElem, gray = False):<o:p></o:p></p>
<p class="MsoNormal"> """ morphological opening """<o:p></o:p></p>
<p class="MsoNormal"> #Mnp = np.array(morph.binary_erosion(imArray, sElem))<o:p></o:p></p>
<p class="MsoNormal"> #Mnp = np.array(morph.binary_dilation(Mnp, sElem))<o:p></o:p></p>
<p class="MsoNormal"> if gray:<o:p></o:p></p>
<p class="MsoNormal"> imArray1 = np.array(morph.grey_dilation(imArray, structure = sElem))<o:p></o:p></p>
<p class="MsoNormal"> imArray2 = np.array(morph.grey_erosion(imArray1, structure = sElem), dtype = 'float32')<o:p></o:p></p>
<p class="MsoNormal"> del imArray1<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal"> return imArray2<o:p></o:p></p>
<p class="MsoNormal"> else:<o:p></o:p></p>
<p class="MsoNormal"> Mnp1 = np.array(morph.binary_dilation(imArray, sElem))<o:p></o:p></p>
<p class="MsoNormal"> Mnp2 = np.array(morph.binary_erosion(Mnp1, sElem), dtype = 'float32')<o:p></o:p></p>
<p class="MsoNormal"> del Mnp1<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal"> return Mnp2<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal"> def calcInfoChange(self,imArray1, imArray2):<o:p></o:p></p>
<p class="MsoNormal"> """calculate entropy of an image"""<o:p></o:p></p>
<p class="MsoNormal"> diff = imArray1 - imArray2<o:p></o:p></p>
<p class="MsoNormal"> m_norm = sum(abs(diff)) #Manhattan Norm<o:p></o:p></p>
<p class="MsoNormal"> z_norm = norm(diff.ravel(), 0) #Zero norm<o:p></o:p></p>
<p class="MsoNormal"> del diff<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal"> return (m_norm, z_norm)<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal"><span style="color:black">Joe McGlinchy <b>|</b> </span><b><span style="color:#00B050">Imagery Scientist</span></b><span style="color:#00B050"><o:p></o:p></span></p>
<p class="MsoNormal"><b><span style="color:black">Database Services – 3D and Imagery Team<o:p></o:p></span></b></p>
<p class="MsoNormal"><b><i><span style="color:#00B050">ESRI</span></i></b><span style="color:black">
<b>||</b> </span><b><i><span style="color:#00B050">380 New York St.</span></i></b><span style="color:black">
<b>||</b><i> </i></span><b><i><span style="color:#00B050">Redlands, CA</span></i></b><i><span style="color:#00B050">
<b>92373</b></span></i><span style="color:black"> <b>||</b> </span><b><i><span style="color:#00B050">USA</span></i></b><span style="color:black"><o:p></o:p></span></p>
<p class="MsoNormal"><i><span style="color:black">T 909-793-2853, ext. 4783<o:p></o:p></span></i></p>
<p class="MsoNormal"><o:p> </o:p></p>
</div>
</body>
</html>