KMarkert / servir-vic-training

Training materials for the VIC hydrologic model setup developed for the @SERVIR program
GNU General Public License v3.0
36 stars 32 forks source link

max(nbands) in format_snow_params returning empty sequence #29

Open Saadi4469 opened 3 years ago

Saadi4469 commented 3 years ago

Hi Dr. Markert,

I have a VIC mask grid raster and 30 meter DEM, and I am trying to create the VIC 4.2.d elevation band file by using the following code Now, I am getting an empty sequence error at line 188 of the original "format_snow_params.py"script, what can I do to fix this issue please? My data can be accessed from here https://app.box.com/s/8qg5e8qif9r8hqr7han440zdw61mlxfn

On Conda prompt I wrote:

python format_snow_params.py ~/VIC_GRID_Rstr1.tif ~/DEM.tif ~/snow_params.txt 492 Where VIC_GRID_Rstr1.tif is the lower resolution grid raster (basinMask), DEM is a high resolution raster(elvHiRes), snow_params.txt (outsnow) is the output file to generated and lastly 492 (interval) is the elevation interval in meters.

While running the script I get the following error:

F``` ile "format_snow_params.py", line 218, in main() File "format_snow_params.py", line 212, in main format_snow_params(sys.argv[1],sys.argv[2],sys.argv[3],sys.argv[4]) File "format_snow_params.py", line 70, in format_snow_params elvhires[np.where(elvhires<0)] = np.nan ValueError: cannot convert float NaN to integer


I exported (using ArcGIS) the DEM (DEM_0-n2) again in .tif format and set no data to 0, but the Float Nan error still remains. And if I remove the line

`elvhires[np.where(elvhires<0)] = np.nan
`then I get the error at line 188, which is:

`ValueError: max() arg is an empty sequence
`                                        
  Even when I view the DEM values in array form in Python, the actual values of the DEM start at column 910-1134. The rest of the columns are all zeros. I am attaching an screenshot of the array just for reference. The script can also be downloaded from the same link. I also tried the script both in Python 2.7 and 3.6, but the error still remains the same. Please let me know if there are data access issues.
![DEM_Array](https://user-images.githubusercontent.com/55733842/97129997-775c3280-16fd-11eb-9c68-6a7644607550.PNG)
KMarkert commented 3 years ago

The point where the program threw an error was where it was trying to mask elevation values less that 0 as no data. This is done to prevent no data values from skewing the elevation band information. This error basically states a nan value cannot be used within a integer array

The first thing I would suggest trying is to convert your elevation data as float. You can do this either in a GIS application (ArcMap or QGIS) or you can add a line in the script where the elevation data is read to convert it to float (elvhires = elvhires.astype(np.float))

Please let me know if that solves your error.

Saadi4469 commented 3 years ago

Even after converting the DEM into float, I am getting the same error of empty sequence. I have also attached a screenshot showing the description of the DEM (in float) array, and as you can see the array is in Float 32. I have even uploaded the Float DEM (Float_DEM_0-n.tif) on the share box link in case you want to test the data. Also, the DEM does not have any values less than 0, because upon exporting the DEM in .tif format, I set no-data values to 0. Float_DEM_Array.

Saadi4469 commented 3 years ago

Dr. Markert, I have an update, the problem was with my Grid raster which had Unique Grid IDs stored as pixel values, hence the condition in line 113 was coming out to be false, therefore I was getting an empty array. I solved the problem by creating a binary grid. This did solve the problem. However, some of the fractional areas are coming out to be less than 1, why is that so? It is counting the right amount of pixels (which are 256 in number with a cell size of 0.0625*0.0625 degrees). I am attaching the output snow parameter file for further reference. I have uploaded the correct grid raster on the link that I shared in case you test the data. [outsnow(with binary valuesi.e grid cell ID is not there).txt](https://github.com/KMarkert/servir-vic-training/files/5442711/outsnow.with.binary.values(grid.cell.ID.not.stored as pixel value).txt)

out_snow

KMarkert commented 3 years ago

Glad you were able to solve the problem! Yes, the grid raster is expected to only be a value of 0 or 1.

I imagine you are getting a fraction less than one because not all of the pixels are being counted in the statistics. Are you using a DEM that is clipped to the basin? If so, then that may be causing the issue. The program is looking for pixels that have been classified as an elevation band between the min and max values (https://github.com/KMarkert/servir-vic-training/blob/master/scripts/format_snow_params.py#L149) and counts only those bands up to maxbands. So, if there are really large values in the data on the edges that can be throwing it off.

This actually is less of a problem though because the VIC model provide a warning and normalize these fractions for you by dividing each band fraction by the sum of band fractions for the grid cell when executing the model. So the fractions for grid cell 1 will be 0.6009 0.3991 on execution. Or if you would like piece of mind, you can add that normalization in the script.

Saadi4469 commented 3 years ago

Hi Dr. Markert,

The script works fine. Now, when I create elevation bands in ArcGIS with the same interval, the distribution of bands across the basin don't match with the grid cell elevation band based distribution created by the script (although, the number of bands is same). In ArcGIS, I am using the Contour tool to create elevation bands. Please advise, that if I want to work with the same elevation bands distribution (as created by the script) in ArcGIS, which tool should I use? Basically, the reason for doing this step in ArcGIS is that, this would tell me where exactly these elevation bands lie in a grid cell. Elevation_Bands

Saadi4469 commented 3 years ago

This is because the VIC model does not explicitly tell how these bands are spatially distributed in a grid cell i.e. the model does not tell spatially tell the fractional elevation band coverages in a grid cell. So I want to reconstruct the elevation band information on sub-grid level in ArcGIS.