Latest Contributions by ChrisBater

Rather than use a set, I wonder if you can use numpy's "unique" function, with "return_index=True" if necessary. Check this out: http://docs.scipy.org/doc/numpy/reference/generated/numpy.unique.html Chris ... View more

Try something like this..... import os draftFileName = r'd:\temp\draft.txt' draftFile = open(draftFileName, 'w') #do stuff draftFile.close() if os.path.exists(draftFileName):     os.remove(draftFileName) ... View more

As this is an esri forum this may not be terribly helpful, but if you have access to envi/idl you should check out the modis conversion toolkit, which you can find here: http://www.exelisvis.com/UserCommunity/CodeLibrary.aspx Chris ... View more

I'd try (if you haven't already), removing EACH variable from memory until the lock file (hopefully) disappears. I'd start with typing "del row" and "del rows" into the console. Those lines should error since it looks like they are being properly removed in the script. Keep removing variables (like "del pnt"). "dir()" will give you the full list of all variable in memory. Let me know if you find your culprit. Good Luck! Try "del arcpy" too. ... View more

You may want to check out matplotlib, which is a python library for plotting. http://matplotlib.sourceforge.net/gallery.html The gallery includes code with the examples. Chris ... View more

None should not be in quotes. ... View more

When I import a module I've written, I also like to reload the module to update any changes I may have made to it, so: sys.path.append(r'modulePath') import module reload(module) ... View more

If you ever find yourself using big arrays you may find you need to delete them to free up memory. ... View more

You can convert the raster to a numpy array, flatten and sort the array, then determine the percentiles based on the flattened array's size. It would look something like this: def firstPercentile(array):     #remove no data values, in this case zeroes, returning a flattened array     print 'removing no data values and flattening array.....'     flatArray = array[array != 0.0]     print 'sorting array....'     flatArray = np.sort(flatArray)         #report some summary stats     print 'min = ', np.min(flatArray)       print 'median = ', flatArray[int(np.size(flatArray) * 0.50)]     print 'max = ', np.max(flatArray)      firstPercentile = flatArray[int(np.size(flatArray) * 0.01)]     print '1st percentile = ', firstPercentile Chris ... View more

Probably the easiest way forward is to use image differencing, where you calculate a vegetation index for each image (e.g. NDVI), difference the images (e.g. TT1-T2), then apply threhsolds to the change images to classify disturbance. A google search on change detection and image differencing might help with more specifics. ... View more

This might help get you started with gdal.... Based on code I found here: http://www.gis.usu.edu/~chrisg/python/ ''' --------------------------------------------------------------------------- --------------------------------------------------------------------------- program: TOA to NDVI.py Chris Bater GIS and Remote Sensing Technologist Wildfire Resource Information Unit Wildfire Management Branch Forestry Division Alberta Sustainable Resource Development 9th floor Great West Life Building 9920 - 108 Street Edmonton, AB T5K 2M4 Ph: 780-422-0669 Cell: 780-446-1899 Email: chris.bater@gov.ab.ca Date: September 2011 Purpose: Converts Landsat-5 TOA reflectance data to NDVI --------------------------------------------------------------------------- ---------------------------------------------------------------------------- ''' #---------------------------------------------------------------------- #Identify the directory containing the landsat data workspace = r"C:\_LOCALdata\temp" #---------------------------------------------------------------------- #---------------------------------------------------------------------- #Import and load modules #---------------------------------------------------------------------- # Import system modules import os, sys from osgeo import gdal import numpy as np import numpy.ma as ma #register the GDAL drivers gdal.AllRegister() driver = gdal.GetDriverByName('GTiff') #not sure if this is necessary driver.Register()#not sure if this is necessary print "" print "All systems go" print "Hit CTRL + c at any time to stop processing" print "" #---------------------------------------------------------------------- #Iterate through the landsat scenes and convert to TOA #---------------------------------------------------------------------- try:         files = os.listdir(workspace)     for f in files:         if not f.endswith('tif'):             continue         name = f[:-4]                 #open the landsat TOA file         TOA = gdal.Open(workspace + "/" + f)         #get the image size         rows = TOA.RasterYSize         cols = TOA.RasterXSize         bands = TOA.RasterCount                 if bands != 6:             continue         print "processing " + f                 #create an empty tif file         outData = driver.Create(os.path.join(workspace, name + "_NDVI.tif"), cols, rows, 1, gdal.GDT_Float32)         outBand = outData.GetRasterBand(1)         #assign the bands to variables         b3 = TOA.GetRasterBand(3)         b4 = TOA.GetRasterBand(4)                         #set the initial blocksize for analysis         xBlockSize = 5000         yBlockSize = 5000         for i in range(0, rows, yBlockSize):             if i + yBlockSize < rows:                 numRows = yBlockSize             else:                 numRows = rows - i             for j in range(0, cols, xBlockSize):                 if j + xBlockSize < cols:                     numCols = xBlockSize                 else:                     numCols = cols - j                          #read bands into arrays                 d3 = b3.ReadAsArray(j, i, numCols, numRows).astype(np.float32)                 d4 = b4.ReadAsArray(j, i, numCols, numRows).astype(np.float32)                                 #mask no data values (no data = 0)                 m3 = ma.masked_equal(d3, 0)                 m4 = ma.masked_equal(d4, 0)                              del d3, d4                                 ndvi = (m4 - m3) /(m4 + m3)                 del m3, m4                         outBand.WriteArray(ndvi, j, i)                 # flush data to disk, set the NoData value and calculate stats         outBand.FlushCache()         outBand.SetNoDataValue(0)         stats = outBand.GetStatistics(0, 1)         # georeference the image and set the projection         outData.SetGeoTransform(TOA.GetGeoTransform())         outData.SetProjection(TOA.GetProjection())         #memory management         del TOA, outData, outBand                     print 'script complete'                         except Exception, e:     # If an error occurred, print line number and error message     tb = sys.exc_info()[2]     print "Line %i" % tb.tb_lineno     print e.message ... View more

I've been struggling with this as well. My work-around has been to use gdal. However, if I'm forced to use arcpy, then I crank out the array as a raster with no spatial reference, then use "rescale," "shift," and "define projection" to put the final surface back to where it should be. I wrote a function to do this, but it's definitely not optimal. There needs to be a better, more in depth usage example included in the numpyarraytoraster tool's help doc. ... View more

I imagine because calculating standard deviation on a bunch of integers still requires a floating point number for maximum precision. For example the standard deviation of [1,2,2,2] is 0.5.  You could always run the integer (INT) tool on the output standard deviation raster to convert back. Chris ... View more

Assuming you know "general python coding" (loops, if/then, while, etc.) you can easily put together a meaningful script. Just in case you don't, this site is pretty good: http://www.python-course.eu/index.php ... View more

Can you just run your regression on a subset of your data? Take a systematic sample of your rasters and run your anayses on those. ... View more