Latest Contributions by Leo_KrisPalao

My intention is to stack my 2D array created from my MODIS rasters using numpy dstack.  After that, I hope to implement savitzky-golay filtering algorithm from scipy signal processing so I think it is better to stack my 2D array to access the time series values of a particular pixel. Right now, I am just testing for 10 MODIS rasters (but later on will use my whole dataset -- approx 800 rasters (15 years of data)) only to test how should I create my script. So my initial workflow is to: 1) Access my raster files 2) For each file I need to convert it to a 2D array using RasterToNumPyArray 3) Append all those created 2D arrays to a list, and 4) Implement numpy.dstack() to stack my 2D arrays and hopefully I can begin some analysis. I am getting an error in the 4th step (it says memory error) -- converting my list of 2D arrays to a stack (please see script 1 below) but if I revise my script and load in dstack using the index of each element in my 2D array it runs okay. Could you guide how to approach my problem. I am just building my script so for now its not so much. please see below. Thanks. Script 1 import arcpy from arcpy.sa import * import os, sys import numpy as np from scipy.signal import savgol_filter arcpy.CheckOutExtension("Spatial") rasters = [] ws = 'G:/Test/Raster' for folder, subs, files in os.walk(ws):   for filename in files:   aSrc = arcpy.RasterToNumPyArray(os.path.join(folder,filename))   rasters.append(aSrc) stackRast = np.dstack(rasters) Script 2 import arcpy from arcpy.sa import * import os, sys import numpy as np from scipy.signal import savgol_filter arcpy.CheckOutExtension("Spatial") rasters = [] ws = 'G:/Test/Raster' for folder, subs, files in os.walk(ws):   for filename in files:   aSrc = arcpy.RasterToNumPyArray(os.path.join(folder,filename))   rasters.append(aSrc) stack = np.dstack((rasters[0], rasters[1], rasters[2], rasters[3], rasters[4], rasters[5],rasters[6], rasters[7], rasters[8], rasters[9])) ... View more

Dear ArcPy/Python Users, I am working on four MODIS datasets on the same tile. I reclassified the map into boolean (1 and 0) based on some criteria. My objective is to have a reclassified map which contains some values based on the combination of values for each pixel from my four MODIS dataset. For instance, each pixel of my four rasters would contain one of the possible combinations, i.e., (1,1,0,0) the values are based on the sequence of my rasters -- (pixel value raster 1, pixel value raster 2, ...). How can I convert this to 1100? Another example would be: >>> a = np.array([1,2,3,4]) >>> b = np.array([5,6,7,8]) >>> c = np.array([9,10,11,12]) >>> d = np.array([4,7,8,9]) >>> x      array([[[ 1,  5,  9,  4],           [ 2,  6, 10,  7],           [ 3,  7, 11,  8],           [ 4,  8, 12,  9]]]) convert the items to 1594, 26107, 37118, 48129 then for instance 1594 will be reclassified to say 2120 The reason is that, I have a dictionary variable (see below) which will reclassify the values for each possible combination. So in the end I will have an array of, i.e., from 1100 to 2120. The corresponding values are representative on the number of occurrence and the sequence of that event. For instance the 1100 will be reclassified into 2120 -- this means the event of interest occurred twice in week 1 and week 2. {1000:1100,   0100:1200,   0010:0130,   0001:1400,   1001:2140,   1100:2120,   0110:2230,   0011:2340,   1011:3134,   1101:3124,   1110:3123,   0111:3234,   1111:4111,   0000:1 255,255,255,255 = 0} Here is the code that I am trying to complete and run but I am stuck now: # import arcpy module of arcgis import arcpy, os, sys from arcpy import env # import spatial extension from arcpy.sa import * arcpy.env.overwriteOutput = True # activate spatial extension arcpy.CheckOutExtension("Spatial") # environment settings arcpy.env.extent = "MAXOF" # workspace arcpy.env.workspace = r'raster\directory' rasters = arcpy.ListRasters("*", "tif") # output workspace out_ws = 'output/workspace/' rastersCount = len(rasters) counter = 0 # starting index pixReclass = {1000:1100,               0100:1200,               0010:0130,               0001:1400,               1001:2140,               1100:2120,               0110:2230,               0011:2340,               1011:3134,               1101:3124,               1110:3123,               0111:3234,               1111:4111,               0000:0}    # compute dekad rasters while counter < 4: # I am processing/testing only the first four rasters in my list   rast = []   rastMon = []   for i in range(counter,counter+4):        #print rasters        rast.append(rasters)        print "rasters:{}".format(rast)        print len(rast)      for r in rast:        name, ext = r.split(".")        crit_Out = out_ws + name + "_reclass.tif"           criteria = Raster(r) >= 0001        criteria.save(crit_Out)        print criteria           vRast = arcpy.RasterToNumPyArray(crit_Out)        rastMon.append(vRast)        print "number of array:{}".format(len(rastMon))      zipped = np.dstack((rastMon[0], rastMon[1], rastMon[2], rastMon[3]))      counter += 4 Hope that I made my query clear. Thanks for any advice and suggestion. -Leo Message was edited by: Leo Kris Palao ... View more

Dear Arcpy Users, I am reclassifying my raster data but I cannot get the output that I want. The case is, I have a floating point raster, in which, I want to reclassify it into categorical values. It is quite weird that some of my values are reclassified outside the range of values that I set. For instance, output categorical raster          original raster 1                                                 1.24941356103 1                                                 0.944287159553 1                                                 0.70514940818 0                                                 0.0177492300117 0                                                 2.72532446104e-07 250                                             -6.0889773748 This is the reclassification method that I want to do. -50 to -2.58 = 6 -2.58 to -1.96 = 5 -1.96 to -1.65 = 4 -1.65 to -1.29 = 3 -1.29 to 2.58 = 2 2.58 to 50 = 1 other pixels = No data Here is the code that I am running. import arcpy, os, sys from arcpy import env from arcpy.sa import * arcpy.env.overwriteOutput = True arcpy.CheckOutExtension("Spatial") arcpy.env.extent = "MAXOF" arcpy.env.workspace = r'directory\of\rasters\test' rasters = arcpy.ListRasters("*","TIF") out_ws = 'output/directory/Reclass/' if not os.path.exists(out_ws):   os.makedirs(out_ws) for raster in rasters:   name, ext = raster.split(".")   fileout = out_ws + name + "_reclass" + ".tif"   outReclass = Reclassify(raster, "Value", RemapRange([[-50.0,-2.58,6],[-2.58,-1.96,5],\   [-1.96,-1.65,4],[-1.65,-1.29,3],[-1.29,2.58,2],[2.58,50.0,1]]),"NODATA")      outReclass.save(fileout) ... View more

Hi ArcPy users, I have a bit of problem in my script. I am working with MODIS data, and I already have a stack of raster layers. I want to rescale my MODIS data (stack layers) but whenever I use the raster calculator or the divide function in Math the output is converted into single band. I have 10 bands(rasters) stacked in each of my TIFF files. out_Div = out_ws + "MOD_" + date_parse + "_clip_Geo_rescale" + ".tif" outDivide = Divide(bands_WGS, 10000.00) outDivide.save(out_Div) bands_WGS is a stack of 10 rasters. When I execute this the output is converted into a single band. Not a stack of rasters divided by 10000.00 Hope my query is clear. Thanks in advance for any help. -Leo ... View more

Hi ALL, I am doing a Temperature-Vegetation Dryness Index (TDVI) analysis using Normalized Difference Vegetation Index (NDVI) and Land Surface Temperature (LST) for drought assessment. I need to create a scatterplot of LST and NDVI. Can you help me do this in Python? or at least guide me on how to do this. Thanks, -Leo ... View more

Hi Python Users, I got stuck with my code trying to merge modis rasters (in tiff format) using Mosaic to New Raster. Any help will be much appreciated. I have  time series NDVI rasters from MODIS data for 2014 for two tiles located in separate directory. What I would like to achieve is to automate the merging of the times series NDVI rasters using the Mosaic to New Raster tool in arcpy. MOD_h26v06_2014001_ndvi.tif merged with MOD_h27v06_2014001.tif . . . MOD_h26v06_2014297_ndvi.tif merged with MOD_h27v06_2014297.tif Any help/suggestion is very much appreciated. Below is the code that I am working on. Thanks, -Leo import arcpy, os, sys from arcpy import env from arcpy.sa import * arcpy.env.overwriteOutput = True arcpy.CheckOutExtension("Spatial") raster1 = 'path/toRaster1' arcpy.env.workspace = raster1 tile1 = [os.path.join(raster1, l) for l in arcpy.ListRasters("*2014*", "TIF")] raster2 = 'path/toRaster2' arcpy.env.workspace = raster2 tile2 = [os.path.join(raster2, r) for r in arcpy.ListRasters("*2014*", "TIF")] out_ws = 'path/Out_rast' # (input rasters,outputloc,name w/ ext, coor, pixeltype, cellsize, #of bands # mosaic method, mosaic color map) for h in range(0,46):     parse = l[18:25]     filename = "modis_evi_" + parse + ".tif"     mosaic = arcpy.MosaicToNewRaster_management(Raster(tile1 );Raster(tile2 ), out_ws, filename, "", "16_BIT_SIGNED", "", "1", "LAST","FIRST")     print mosaic ... View more