2

I have 80 .tif files, and each pixel value represents NDVI. The whole stack of images is a time-series. I want to run a regression for each pixel using NDVI and time. I have tried to approach this two ways. The first way is to combine all of my rasters, using the combine tool, and then run the regression by manipulating the dbf table in python, and the second is to convert each pixel to a point, intersect the points, and then run the regression, again just by organizing the data and regressing in python.

The problem is that certain images have clouds masked out as No Data. When the rasters are combined, or when the points are intersected, this results in the final layer having no No Data assigned in all of the locations where any individual input layer had No Data. To get around this I assigned No Data for all of the input rasters as -999, and then I use the same process described above. When I combine all of the rasters though the original data that is not -999 is not retained in the output combine, and when I intersect using the points method every single point is given a value of -999, so again the original data that is not -999 is not retained. This makes no sense to me why this would be happening and was hoping to get some insight.

The entire process for the combine method I use can be seen with this code:

import arcpy 
from arcpy.sa import *
import os
arcpy.env.overwriteOutput = True
arcpy.CheckOutExtension('Spatial')

#assign no data to -999
folders=r'F:\Sheyenne\Final_Imagery\Sheyenne'
out=r'F:\Sheyenne\Pixel_Regression\no_data'
for folder in os.listdir(folders):
    if not os.path.isdir(os.path.join(out,folder.split('_')[0])):
        os.mkdir(os.path.join(out,folder.split('_')[0]))
    arcpy.env.workspace=os.path.join(folders,folder)
    rasters=arcpy.ListRasters('*.tif')
    for raster in rasters:
        ras = Raster(raster)
        outraster=Con(IsNull(raster),-999, raster)
        outraster.save(os.path.join(out,folder.split('_')[0], raster))

#combine rasters for the month of May
folders=r'F:\Sheyenne\Pixel_Regression\no_data'                                                                                                                                     '
out=r'F:\Sheyenne\Pixel_Regression\Sheyenne_combined'
for folder in os.listdir(folders):
    if not os.path.isdir(os.path.join(out,folder.split('_')[0])):
        os.mkdir(os.path.join(out,folder.split('_')[0]))
    arcpy.env.workspace=os.path.join(folders,folder)
    rasters=arcpy.ListRasters('*.tif')
    rasters2=[]
    for raster in rasters: 
        if raster[13:15] == '05':
            rasters2.append(raster)
            rasters=[os.path.join(os.path.join(folders,folder), i) for i in rasters2] 
    outcombine=Combine(rasters)
    outcombine.save(os.path.join(out, folder.split('_')[0], '05_Combine.tif'))
print "Done"

If someone wants to see the code for the point method just let me know and I will post it.

EDIT:

I got it so that the original data points are retained when converting to points. But now when I intersect the output is empty. I see it temporarily populate in my folder but then it disappears. In Desktop mode I get the error output is empty. The code I am using for this method is as follows:

folders=r'F:\Sheyenne\Pixel_Regression\no_data'   
out=r'F:\Sheyenne\Pixel_Regression\raster_to_point\Sheyenne'
for folder in os.listdir(folders):
    if not os.path.isdir(os.path.join(out,folder.split('_')[0])):
        os.mkdir(os.path.join(out,folder.split('_')[0]))
    arcpy.env.workspace=os.path.join(folders,folder)
    outDir=os.path.join(out,folder.split('_')[0])
    rasters=arcpy.ListRasters('*.tif')
    for raster in rasters:
        rastername=raster[:-4]
        arcpy.RasterToPoint_conversion(raster, os.path.join(outDir,rastername),"Value")
print "Done converting to point"
#    

folders=r'F:\Sheyenne\Pixel_Regression\raster_to_point\Sheyenne'
out=r'F:\Sheyenne\Pixel_Regression\Intersected\Sheyenne'
for folder in os.listdir(folders):
    if not os.path.isdir(os.path.join(out,folder)):
        os.mkdir(os.path.join(out,folder))
    arcpy.env.workspace=os.path.join(folders,folder)
    outDir=os.path.join(out,folder)  
    shapefiles=arcpy.ListFeatureClasses('*.shp')
    shapefiles2=[]
    for shape in shapefiles:
         if shape[13:15] == '05':
            shapefiles2.append(shape)
            shapefiles=[os.path.join(os.path.join(folders,folder), i) for i in shapefiles2]  
    arcpy.Intersect_analysis(shapefiles, os.path.join(outDir, '05_Instersect'),"ALL","#","INPUT")

EDIT:

The method I used, which works, is to convert the rasters to points (after assigning all nodata to -999) then combine all the 80 resulting text files together in python pandas. From there just refill -999 as NaN, do a little bit of data pivoting and run the regression using scipy and a groupby. Can't believe this isn't just a tool at this point but this works to run pixel level regression.

  • This post,gis.stackexchange.com/questions/65787/…, shows something similar to what you are trying to achieve but I do not think will remedy -999 problem that you have. In any case if you do not treat NoData as no data or null, your results will be biased since -999 has a value (ratio). More on the latter matter is somewhat touched in stats.stackexchange.com/questions/225175/…. – fatih_dur Jul 27 '16 at 1:20
  • I will just turn -999 back to Null when I am organizing my data later, it won't actually be used in the regression. I am also going to need RMSE, r value and p-value so I don't think just the raster calculator will work for me. Thank you though and I will take a look – Stefano Potter Jul 27 '16 at 2:30
  • Brute force approach is to convert clouds free raster to points an d sample rasters. How big is the raster? Will it fit to excel ncols*nrows – FelixIP Jul 27 '16 at 5:20
  • That is exactly what I am attempting to do. Each raster is only 3.5 kb and I have 80 of them. When I convert to point every single point is given the value of -999 though. It's really not that much data, I just need to intersect them to retain the spatial locations across the rasters, unless there is another way to do that. – Stefano Potter Jul 27 '16 at 16:38

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.