POST
|
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
01-14-2014
10:07 AM
|
0
|
0
|
182
|
POST
|
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
07-25-2012
09:39 AM
|
0
|
0
|
358
|
POST
|
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
07-14-2012
07:15 AM
|
0
|
0
|
320
|
POST
|
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
06-19-2012
11:18 AM
|
0
|
0
|
496
|
POST
|
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
06-15-2012
06:40 AM
|
0
|
0
|
244
|
POST
|
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
05-30-2012
07:57 AM
|
0
|
0
|
1158
|
POST
|
If you ever find yourself using big arrays you may find you need to delete them to free up memory.
... View more
04-20-2012
09:04 AM
|
0
|
0
|
224
|
POST
|
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
04-13-2012
06:56 AM
|
0
|
0
|
208
|
POST
|
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
04-05-2012
12:29 PM
|
0
|
0
|
215
|
POST
|
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
01-27-2012
08:25 AM
|
0
|
0
|
466
|
POST
|
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
01-25-2012
05:57 AM
|
0
|
0
|
466
|
POST
|
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
01-23-2012
05:40 AM
|
0
|
0
|
263
|
POST
|
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
01-20-2012
06:15 AM
|
1
|
0
|
4769
|
POST
|
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
10-28-2011
10:14 AM
|
0
|
0
|
233
|
Title | Kudos | Posted |
---|---|---|
1 | 01-20-2012 06:15 AM | |
1 | 07-08-2011 12:53 PM |
Online Status |
Offline
|
Date Last Visited |
11-11-2020
02:23 AM
|