Correlation between rainfall and temperature anomaly in Python

Question

This question has been branched from this question: https://community.esri.com/thread/198969-trend-analysis-through-time-series-of-raster-data to give it more visibility

Dear Xander Bakker‌ and Dan Patterson thank you for the awesome discussion.

I have different case. Let say I have raster timeseries for monthly rainfall anomaly data from 1981-2017 for 1 province. And then I have sea surface temperature anomaly data for the same periods and from 1 location only (text data).

I want to see the correlation between both data (rainfall and temperature anomaly) for each pixel in raster data. Is it possible to do this using above script? I want to have sets of raster (slope, intercept, r value, p value and std error) as an output.

From this output I hope I can see if the temperature increase in the Pacific, which area will experience more/less rainfall.

Thanks, Benny

XanderBakker · Answer

Hi Benny Istanto ,

I gave this some more thought and still come to the conclusion that having 444 regression rasters will make it hard to define a conclusion of the analysis. If the sequence of rasters is indeed a sequence of values through time, maybe it would be better to loop through the 444 raster, on a pixel level determine the difference between the rainanom raster and the corresponding sstanom raster and determine the regression values for the sequence of 444 differences on a pixel level. Does that make sense? With some modifications to the script this should be possible to achieve.

XanderBakker · Answer

Hi Benny Istanto ,

Thanks for sharing the simplified data. I haven't had the time to look into it, since it is a "little" hectic here at the end of the year, but I haven't forgotten about it. I dowwloaded the data and will start soon to give it a try. Thanks for your patience.

I'm not sure if a regression based on two values will yield a good result, but I would have to look into this to see what the result will be.

As you can see in the "other" thread (https://community.esri.com/thread/198969-trend-analysis-through-time-series-of-raster-data#comment-711148 ) you can access the raster(s) with arcpy.RasterToNumPyArray(ras_path), get the number of rows and columns using the shape of the array like this:

rows = ras_np.shape[0]
 cols = ras_np.shape[1]

Based on this you can create a loop through the rows and columns to access each pixel in multiple rasters (two in your case) and do something with these values and write the  result to a new numpy array which will be converted to a raster at the end. I will have to look into the data to see if a regression based on two values makes sense.

BennyIstanto · Answer

Hi Xander Bakker‌

I have tried to simplify the input data. See attached: rainfall anomaly and sst anomaly

Both folder (rainfall and sst) have 444 geotiff. For SST, I have created constant raster, and its have same dimension with rainfall anomaly, the value is following data in sst_anom.csv in my previous answer.

Using your script in the different thread, how to do pixel-wise regression between two sets of raster timeseries?