POST
|
Hi Dan, For the final output, I am would like to have a raster with information on SLOPE, INTERCEPT, and Maximum rainfall threshold that could trigger a flood. So if I have rainfall forecast in a raster, then I could use above mentioned data to calculate the probability if each pixel will flooded or not. Now, I am still preparing all the data requirement for that, in the question above I have attached example data (rainfall and percent of water occurrence). As the pixel size is 0.05 deg ~ to 5.6km/pixels and by doing 3x3 or 5x5 group of cells, in my opinion is enough to see the correlation between both data. Most of extreme rainfall followed by flood event are occurred locally, then 3x3 ~ 282.24 sqkm or 5x5 ~ 784 sqkm is a good start to analyse it. Do you have any script that ready to use? I also review some script related to regression written by you and/or Xander Bakker in other thread, but doesn't consider 3x3 or 5x5 neighbour value. I found similar problem that already answered using R approach: https://gis.stackexchange.com/questions/278979/linear-regression-between-every-3×3-pixels-between-two-rasters-using-r
... View more
08-07-2020
03:20 AM
|
0
|
1
|
938
|
POST
|
I have two raster with same resolution and dimension: rainfall (mm) and water-occurrence (%). I want to create a new raster of SLOPE (a) and INTERCEPT (b) by performing linear regression between every 3×3 pixels of the two rasters (rainfall and water-occurrence), such that each pixel of the SLOPE and INTERCEPT will hold the regression slope and intercept value obtained from linear regression of the corresponding 3×3 pixels in rainfall and water-occurrence that surround that pixel. Is there any example solution using ArcGIS or python/arcpy related above problem? Example data in the attachment.
... View more
08-07-2020
12:16 AM
|
0
|
3
|
1045
|
POST
|
Dear Xander Bakker Sorry to open this thread again. Let say all the input file has date information in filename ie. r001_NPP_YYYY.MM.tif From your script, is it possible to make the time dimension enabled based on date in the filename when saving into netCDF?
... View more
07-05-2020
08:33 AM
|
0
|
1
|
561
|
POST
|
Yes, the previous script to calculate consecutive dry days works well using Daily data with similar naming convention chirps-v2.0.1981.01.01.tif chirps-v2.0.1981.01.02.tif chirps-v2.0.1981.01.03.tif chirps-v2.0.1981.01.04.tif .. chirps-v2.0.2020.05.31.tif punctuation in filenames already eliminated in line 80. I think the problem because the script read input file with date YYYYMMDD information, and the input is just YYYYMMD => the D is not Date but Dekad with value 1, 2 and 3. I dont know how to modified that
... View more
06-02-2020
06:45 PM
|
0
|
0
|
1670
|
POST
|
Yes correct, because I am not using Daily data, but Dekad or 10 days data with naming convention chirps-v2.0.1981.01.1.tif chirps-v2.0.1981.01.2.tif chirps-v2.0.1981.01.3.tif chirps-v2.0.1981.02.1.tif chirps-v2.0.1981.02.2.tif .. .. chirps-v2.0.2020.04.3.tif Dekad 1 in chirps-v2.0.1981.01.1.tif means rainfall accumulation from Date 1 to 10, and Dekad 2 from 11 to 20, Dekad 3 from 21 to 28/29/30/31
... View more
06-02-2020
05:45 PM
|
0
|
2
|
1670
|
POST
|
There is print statement line 48. Or do you mean another print statement?
... View more
06-02-2020
05:13 PM
|
0
|
4
|
1670
|
POST
|
Hi Dan, thanks. I have updated my question and include a screenshot from cmd window.
... View more
06-02-2020
04:51 PM
|
0
|
6
|
1670
|
POST
|
Hi, with limited knowledge on python and arcpy, I tried to modified script from former staff in my office with no success. The existing script is: Try to calculate number of consecutive days without rainfall (rainfall < 1mm) based on daily CHIRPS data with naming convention "chirps-v2.0.YYYY.MM.DD.tif" The calculation required number of days without rain from the previous day. If this data not available, it will skipped to next data. Let say, the data start from 1 Jan 1981, then the script will skip this date. And continue the process for the next day (2 Jan 1981) with following process: IF(chirps-v2.0.1981.01.01.tif <1, 1, 0) and save as cli_chirps_dslr_19810102.tif IF(chirps-v2.0.1981.01.02.tif <1, cli_chirps_dslr_19810102.tif + 1, cli_chirps_dslr_19810102.tif = 0) and save as cli_chirps_dslr_19810103.tif and so on. And now I would like to use the script and modified some part to calculate similar approach with rainfall > 50mm, but using Dekad CHIRPS data with naming convention chirps-v2.0.1981.01.1.tif chirps-v2.0.1981.01.2.tif chirps-v2.0.1981.01.3.tif chirps-v2.0.1981.02.1.tif chirps-v2.0.1981.02.2.tif .. .. chirps-v2.0.2020.04.3.tif I got an error because the script try to read based on Date information, after finished process chirps-v2.0.1981.01.03.tif the scipt will continue to find chirps-v2.0.1981.01.04.tif and couldn't find the data then error. Any reference or example on how to make it work? Below is the error, script and few example of the data (subset from area with smaller size). .... found chirps-v2.0.2001.01.2.tif in the dekad rainfall folder found chirps-v2.0.2016.07.2.tif in the dekad rainfall folder found chirps-v2.0.1994.11.2.tif in the dekad rainfall folder found chirps-v2.0.2000.01.2.tif in the dekad rainfall folder found chirps-v2.0.2017.07.2.tif in the dekad rainfall folder found chirps-v2.0.1995.11.2.tif in the dekad rainfall folder execute first rainy data from date 1981011 creating precipgt50 data 1981-01-02 from dekad rainfall data from 1981-01-01 1981-01-02 Processing chirps_precipgt50_1981012.tif file chirps_precipgt50_1981012.tif is created start reading existing PRECIPgt50 Dataset.... looking for file with naming chirps_precipgt50_YYYYMMD found chirps_precipgt50_1981012.tif in the PRECIPgt50 folder PRECIPgt50 Data found. Looking for latest PRECIPgt50 Data... next dekad data is available... start processing next PRECIPgt50... Processing PRECIPgt50 from Z:\Temp\DryWetSeason\Output\chirps_precipgt50_1981012.tif and Z:\Temp\DryWetSeason\Dekad\chirps-v2.0.1981.01.2.tif PRECIPgt50 File chirps_precipgt50_1981013.tif is created next dekad data is available... start processing next PRECIPgt50... Processing PRECIPgt50 from Z:\Temp\DryWetSeason\Output\chirps_precipgt50_1981013.tif and Z:\Temp\DryWetSeason\Dekad\chirps-v2.0.1981.01.3.tif PRECIPgt50 File chirps_precipgt50_1981014.tif is created next dekad data is not available. Exit... next dekad data is not available. Exit... next dekad data is not available. Exit... next dekad data is not available. Exit... next dekad data is not available. Exit... next dekad data is not available. Exit... next dekad data is available... start processing next PRECIPgt50... Processing PRECIPgt50 from Z:\Temp\DryWetSeason\Output\chirps_precipgt50_19810110.tif and Z:\Temp\DryWetSeason\Dekad\chirps-v2.0.1981.01.1.tif Traceback (most recent call last): File "PRECIPgt50.py", line 118, in <module> list = create_PRECIPgt50('Z:\\Temp\\DryWetSeason\\Output','Z:\\Temp\\DryWetSeason\\Dekad', 50) File "PRECIPgt50.py", line 107, in create_PRECIPgt50 execute_PRECIPgt50(last_date, _tiffolder, _PRECIPgt50_folder, threshold) File "PRECIPgt50.py", line 50, in execute_PRECIPgt50 outPRECIPgt50Con = Con(Raster(next_dekaddata) < int(threshold), Raster(last_precipgt50file)+1, 0) RuntimeError: ERROR 000732: Input Raster: Dataset Z:\Temp\DryWetSeason\Output\chirps_precipgt50_19810110.tif does not exist or is not supported Z:\Temp\DryWetSeason> import os
import arcpy
from arcpy.sa import *
from datetime import date, timedelta
# to execute first PRECIPgt50 data
def execute_first_PRECIPgt50(_tiffolder, _PRECIPgt50Folder, threshold):
sr = arcpy.SpatialReference(4326)
print("looking at the first dekad rainfall data in tif folder...")
dekad_list = create_dekad_List(_tiffolder)
first_date = min(dekad_list)
print("execute first rainy data from date "+first_date)
first_data_name = 'chirps-v2.0.{0}.{1}.{2}.tif'.format(first_date[0:4], first_date[4:6], first_date[6:7])
first_dekad_data = os.path.join(_tiffolder, first_data_name)
dekad_Date = date(int(first_date[0:4]), int(first_date[4:6]), int(first_date[6:7]))
precipgt50_date = dekad_Date + timedelta(days=1)
print("creating precipgt50 data "+str(precipgt50_date)+ " from dekad rainfall data from "+str(dekad_Date))
PRECIPgt50Year = str(precipgt50_date.year)
PRECIPgt50month = str(precipgt50_date.month)
PRECIPgt50day = str(precipgt50_date.day)
print(str(precipgt50_date))
PRECIPgt50Filename = 'chirps_precipgt50_{0}{1}{2}.tif'.format(PRECIPgt50Year.zfill(4), PRECIPgt50month.zfill(2), PRECIPgt50day.zfill(1))
print("Processing "+PRECIPgt50Filename)
arcpy.CheckOutExtension("spatial")
outCon = Con(Raster(first_dekad_data) < int(threshold), 1, 0)
outCon.save(os.path.join(_PRECIPgt50Folder, PRECIPgt50Filename))
arcpy.DefineProjection_management(os.path.join(_PRECIPgt50Folder, PRECIPgt50Filename), sr)
arcpy.CheckInExtension("spatial")
print("file " + PRECIPgt50Filename + " is created")
# to execute next PRECIPgt50 data
def execute_PRECIPgt50(_lastdate, _tiffolder, _PRECIPgt50_folder, threshold):
sr = arcpy.SpatialReference(4326)
date_formatted = date(int(_lastdate[0:4]), int(_lastdate[4:6]), int(_lastdate[6:7]))
last_precipgt50name = 'chirps_precipgt50_{0}'.format(_lastdate)
last_precipgt50file = os.path.join(_PRECIPgt50_folder, last_precipgt50name)
next_dekadname = 'chirps-v2.0.{0}.{1}.{2}.tif'.format(_lastdate[0:4], _lastdate[4:6], _lastdate[6:7])
next_dekaddata = os.path.join(_tiffolder, next_dekadname)
if arcpy.Exists(next_dekaddata):
print("next dekad data is available...")
print("start processing next PRECIPgt50...")
new_precipgt50_date = date_formatted + timedelta(days=1)
PRECIPgt50Year1 = str(new_precipgt50_date.year)
PRECIPgt50month1 = str(new_precipgt50_date.month)
PRECIPgt50day1 = str(new_precipgt50_date.day)
new_precipgt50_name = 'chirps_precipgt50_{0}{1}{2}.tif'.format(PRECIPgt50Year1.zfill(4), PRECIPgt50month1.zfill(2), PRECIPgt50day1.zfill(1))
print("Processing PRECIPgt50 from "+last_precipgt50file+" and "+next_dekaddata)
arcpy.CheckOutExtension("spatial")
outPRECIPgt50Con = Con(Raster(next_dekaddata) < int(threshold), Raster(last_precipgt50file)+1, 0)
outPRECIPgt50Con.save(os.path.join(_PRECIPgt50_folder, new_precipgt50_name))
arcpy.DefineProjection_management(os.path.join(_PRECIPgt50_folder, new_precipgt50_name), sr)
arcpy.CheckInExtension("spatial")
print("PRECIPgt50 File "+new_precipgt50_name+" is created")
else:
print("next dekad data is not available. Exit...")
# to check if there is PRECIPgt50 data in output folder
def create_PRECIPgt50_List(_PRECIPgt50_folder):
print("start reading existing PRECIPgt50 Dataset....")
print("looking for file with naming chirps_precipgt50_YYYYMMD")
PRECIPgt50_Date_List=[]
for PRECIPgt50_Data in os.listdir(_PRECIPgt50_folder):
if PRECIPgt50_Data.endswith(".tif") or PRECIPgt50_Data.endswith(".tiff"):
print("found " + PRECIPgt50_Data + " in the PRECIPgt50 folder")
parse_String = PRECIPgt50_Data.split('_')
PRECIPgt50_Data_Date = parse_String[2]
PRECIPgt50_Date_List.append(PRECIPgt50_Data_Date)
return PRECIPgt50_Date_List
# to check input data
def create_dekad_List(_tif_folder):
print("start reading list of dekad rainfall data....")
print("looking for file with naming chirps-v2.0.YYYY.MM.DD")
Dekad_Date_List=[]
for Dekad_Data in os.listdir(_tif_folder):
if Dekad_Data.endswith(".tif") or Dekad_Data.endswith(".tiff"):
print("found " + Dekad_Data+ " in the dekad rainfall folder")
parse_String = Dekad_Data.split('.')
Dekad_Data_Date = parse_String[2]+parse_String[3]+parse_String[4]
Dekad_Date_List.append(Dekad_Data_Date)
return Dekad_Date_List
# to run the script
def create_PRECIPgt50(_PRECIPgt50_folder, _tiffolder, threshold):
PRECIPgt50_Date_List = create_PRECIPgt50_List(_PRECIPgt50_folder)
Dekad_list = create_dekad_List(_tiffolder)
# if there is no PRECIPgt50 data, creating new PRECIPgt50 data
if len(PRECIPgt50_Date_List)==0:
print("No PRECIPgt50 Data found...")
print("creating first PRECIPgt50 data...")
execute_first_PRECIPgt50(_tiffolder, _PRECIPgt50_folder, threshold)
PRECIPgt50_Date_List = create_PRECIPgt50_List(_PRECIPgt50_folder)
# if there is PRECIPgt50 data
print("PRECIPgt50 Data found. Looking for latest PRECIPgt50 Data...")
#Check last PRECIPgt50 available
last_date = max(PRECIPgt50_Date_List)
#Check last dekad data availabke
max_dekad_date = max(Dekad_list)
last_PRECIPgt50_date = date(int(last_date[0:4]), int(last_date[4:6]), int(last_date[6:7]))
last_dekad_date = date(int(max_dekad_date[0:4]), int(max_dekad_date[4:6]), int(max_dekad_date[6:7]))
# process PRECIPgt50 to every dekad data available after last PRECIPgt50 data
while last_dekad_date + timedelta(days=1) > last_PRECIPgt50_date:
execute_PRECIPgt50(last_date, _tiffolder, _PRECIPgt50_folder, threshold)
last_PRECIPgt50_date=last_PRECIPgt50_date+timedelta(days=1)
PRECIPgt50Year2 = str(last_PRECIPgt50_date.year)
PRECIPgt50month2 = str(last_PRECIPgt50_date.month)
PRECIPgt50day2 = str(last_PRECIPgt50_date.day)
last_date='{0}{1}{2}.tif'.format(PRECIPgt50Year2.zfill(4), PRECIPgt50month2.zfill(2), PRECIPgt50day2.zfill(1))
print("All PRECIPgt50 data is available")
# run the function (output folder, input folder, threshold)
list = create_PRECIPgt50('Z:\\Temp\\DryWetSeason\\Output','Z:\\Temp\\DryWetSeason\\Dekad', 50)
... View more
06-02-2020
07:54 AM
|
0
|
9
|
1742
|
POST
|
Thanks for the suggestion Dan, I also found this from GitHub - bkomaki/pydrought: Standardized Precipitation Index (SPI) in Python and looks promising. The input is CSV, we'll see if I can find an example on how to modify the script to be able read raster as input.
... View more
12-05-2019
06:18 AM
|
0
|
0
|
2381
|
POST
|
Hi Dan, thanks for the answer. Currently I convert all the raster precipitation data into csv, and put it together into spreadsheet, so far its work. But its easy for small areas, and make a big problem if calculate for the whole country and for the whole period of data (1981-now).
... View more
12-05-2019
04:03 AM
|
0
|
2
|
2381
|
POST
|
Hello, I usually calculate monthly or 3-monthly Standardized Precipitation Index (SPI) using spreadsheet, based on daily rainfall from weather station. Attached is the spreadsheet example. And now I would like to use satellite precipitation estimate (raster data) to calculate the SPI. The problem is in the spreadsheet I use some function like: FREQUENCY, GAMA.DIST, NORMSINV. My question, is there any similar above function in ArcGIS? If not, I am open to any suggestion.
... View more
12-05-2019
02:54 AM
|
0
|
4
|
2701
|
POST
|
Hello, I have daily rainfall data in milimeters in GeoTIFF format with naming convention chirps_YYYYMMDD.tif (example chirps_20100101.tif), and I also have 1 raster of dry-spell with name dslr_chirps_20091231.tif and both raster have same spatial resolution and extent. Then I would like to calculate dry-spell for 1 Jan 2010, IF(rainfall>1,dslr=0,dslr+1). Using raster calculator I can use this formula: Con("chirps_20100101.tif" > 1, "dslr_chirps_20091231.tif" == 0, "dslr_chirps_20101231.tif" + 1) The raster output will be dslr_chirps_20180101.tif Problem: I have 10 years daily rainfall data and would like to calculate daily dry-spell for all the available data period. How to loop the above calculation using model builder, when the output for each calculation will use as input for the next calculation?
... View more
12-04-2019
07:23 PM
|
0
|
2
|
916
|
POST
|
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?
... View more
12-19-2018
12:20 AM
|
0
|
2
|
1791
|
POST
|
Hi Xander Bakker, Sorry for late reply, I tried to prepared the data and example to better understand the process and output that I want to achieve. So this is two different data and the correlation is not based on the distance. I will explain more details below. I would like to see general sensitivity of rainfall in a country to sea surface temperature (SST) changes of NINO-3.4 region in Pacific. Why Pacific? Because this region is optimal for monitoring El Nino-Southern Oscillation and its impacts in Southeast Asia. Simple regression equation is applied to indicate the correlation between rainfall anomaly in each area to anomaly of SST in Pacific Ocean which represent ENSO signals. Y = a + bX, where: Y = Rainfall anomaly a = Y intercept b = Slope X = SST anomaly From this approach (in this case, I used Timor-Leste as an example) and the map, we can say Timor-Leste is heavily affected by El Niño which is associated with a rise in SST, negatively impacting rainfall in much of the country. To get the map, I tried manual process using spreadsheet - precip_sst_regression.xlsx , I convert monthly rainfall anomaly raster data - tls_precip_anomaly.zip to csv and calculate the regression between rainfall anomaly and SST anomaly - sst_anom.csv I hope this help and how to use your script above in my case?
... View more
11-12-2018
04:22 AM
|
0
|
4
|
1791
|
Online Status |
Offline
|
Date Last Visited |
04-07-2022
04:30 AM
|