KeyError upon Instantiating an instance of the DEMProcessor class --Projection

[tgmueller]

I am running PyDEM (0.2.0) on Python (2.7.14) and trying to follow the instructions on the readme.md file. Below I have listed the commands I used, the output, and the input elevation file. I am getting the KeyError with the third command. It looks like it has something to do with the projection files. Has anyone else had this issue and if so, how did you solve it. I tried to search this error with Google but did not find anything useful.

Here are the commands I typed in the Python Console:

from pydem.dem_processing import DEMProcessor
boy = 'elevationUTM.tif'
dem_proc = DEMProcessor(boy)

Then I receive the messages including the key error at the end:

Traceback (most recent call last):
  File "<input>", line 1, in <module>
  File "C:\Python27\lib\site-packages\pydem\dem_processing.py", line 565, in __init__
    elev, = elev_file.raster_layers
  File "C:\Python27\lib\site-packages\traits\has_traits.py", line 895, in decorator
    self.__dict__[ name ] = result = function( self )
  File "C:\Python27\lib\site-packages\pydem\reader\gdal_reader.py", line 195, in _get_raster_layers
    for raster_band in raster_bands]
  File "C:\Python27\lib\site-packages\pydem\reader\gdal_reader.py", line 143, in _raster_layer_from_raster_band
    layer.grid_coordinates = self.grid_coordinates
  File "C:\Python27\lib\site-packages\traits\has_traits.py", line 895, in decorator
    self.__dict__[ name ] = result = function( self )
  File "C:\Python27\lib\site-packages\pydem\reader\gdal_reader.py", line 130, in _get_grid_coordinates
    wkt = d_wkt_to_name[wkt_]
KeyError: 'PROJCS["WGS_1984_UTM_Zone_16N",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-87],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AUTHORITY["EPSG","32616"]]'

Here is the input elevation model I used (this file was created with ArcGISPro using the project raster command and the projection was WGS_1984_UTM_Zone_16N):
elevationUTM.tif.zip

I have tried the solution recommended by mpu-creare but this did not work. Here is the link for that GitHub Issue. https://github.com/creare-com/pydem/issues/1

[mpu-creare]

Seems like I need to add something to the readme file! pydem only really supports geotiffs using EPSG:4326. It's only been tested on files of that type.

You can probably still do some calculation by first reading your data into a numpy array. Looking at the docstring of DEMProcessor's constructor:

Parameters
        -----------
        file_name : str, np.ndarray, tuple
            If isinstance(str): filename of of elevation data
            If isinstance(np.ndarray): a numpy array containing elevation data
            If isinstance(tuple): (elev, lat, lon), three numpy arrays containing
                elevation data, the latitude (rows), and longitude (columns) of the array.
                Note: only the max and min values of the latitude and longitude inputs are used.
                       and a square N-S E-W aligned array is assumed.
        dx_dy_from_file : bool, optional
            Default True. If true, will extract coordinates from geotiff file
            and use those to calculate the magnitude/direction of slopes.
            Otherwise assumes rectangular (uniform) coordinates. This flag is not
            used if file_name is not a string.
        plotflag : bool, optional
            Default False: If True, will plot debug image. For a large
            file this is not advised.
        """

So in your case you can try:

>>> boy = [some_geotiff_reader_like_rasterio]('elevationUTM.tif')
>>> type(boy)
np.ndarray
>>> dem_proc = DEMProcessor(boy)

Note, this assumes a uniform grid in the x-y directions, and calculates the spacing as follows:

            dX = np.ones(self.data.shape[0] - 1) / self.data.shape[1]  #dX only changes in latitude
            dY = np.ones(self.data.shape[0] - 1) / self.data.shape[0]

where data = boy

[tgmueller]

Hi, Thank you for responding. I will try out your suggestion and I will post my results here.

First, I hope you could clarify something. I believe that TauDEM calculates slope as rise over run and therefore requires Cartesian coordinates. But isn't EPSG:4326 a geographic coordinate system? If so, then I assume PyDEM must be projecting the tif for the analysis? Could you shine some light on this?

Thanks, Tom

[mpu-creare]

You are correct about TauDEM -- in fact that's one of the motivations for PyDEM.

So yes, PyDEM uses geopy to calculate the distances between grid cells, code for that is here: https://github.com/creare-com/pydem/blob/master/pydem/utils.py#L145

[tgmueller]

I took a different approach to this issue. I downloaded elevation data from the national map which were by default in the EPSG: 4269 coordinate system. I was able to project that data to with ArcGIS Pro to the EPSG: 4326 coordinate system. I ran the same pyDEM commands and this eliminated teh KeyError problem. Unfortunately there was a new error: "ImportError: No module named _gdal_array". I pasted screen shots of these steps and posted them below. Any suggestions on how to address this Import Error would be appreciated. Also, I'm wondering if there was a way for pyDEM to be modified to take in the EPSG: 4269 coordinate system so the data won't have to be projected when downloading from the National Map.

I sent to the national map

This is what the data looked like:

I opened up the data in ArcGIS Pro and the data were the EPSG: 4269 as shown below:

Then I took the data into ArcGIS and projected it

The new dataset was in in EPSG: 4326 as required:

These were the files that were created:

It fixed the KeyError! But the bad news is now I have an import Error. I found some help with Google on gis.stackexchange.com but I didn't really understand how to solve this. Let me know if you have suggestions. Thank you!

[mpu-creare]

Hurray, progress!

It looks like it's failing when trying to import gdal , so you're gdal installation is likely broken. It can be a bear getting that to work on windows, sometimes.

Have a look at the Install.md for some hints.

Recently I've been having better luck with the conda-forge anaconda channel... something like:

conda install -c conda-forge gdal

... if you're using conda.

But I can't help much beyond this on a gdal installation.

On you're second question, feel free to put out another issue and I'll mark it as a feature request. No idea when someone will be able to get to it. I'm not that familiar with that part of the pydem, but epsg:4269 might actually work -- you could give it a try once you get the gdal matters sorted out.

[tgmueller]

I was finally able to get get pydem installed and working with Anaconda. Thank you for the suggestions! I"m having some trouble interpreting the results i'm getting. If you have time to help troubleshoot, I have copied the code from a Jupyter Notebook, the details about the elevation file from ArcGIS, and a list of the packages I have installed obtained with the command "pip freeze". Please let me know if you have suggestions.

The elevation model I'm creating has the correct projection. You can see the screen shot below that indicate the details.

I believe I have all the packages installed correctly. See below.

I was also trying to configure TauDEM to work with PyDEM according to the instructions found here: and as shown below:

You can see how I put the executable in taudem_Windows directly as shown below. Can you confirm that this is correct?

[mpu-creare]

@tgmueller

Great! Looks like you've got it working. The output is as expected -- it's pretty verbose and probably makes sense to no one but myself. It basically tells you what part of the file (chunk) it's working on, and what stage of the calculation it's at.

Also, it looks like you put taudem in the correct location. You should now be able to run the tests that come with pydem.

[tgmueller]

@mpu-creare Thank you so much. I was able to get look at a DEM in ArcGIS but it was very clumsy and It made a 6 band image rather than a 5 band image. But you can see from the images below that the wetness index flowpaths match very well with the aerial photo where dark lines indicate drainage areas. What is the best way to create geo-tif's in your experience? Thanks, Tom.

[mpu-creare]

We actually have a convenience function for that: https://github.com/creare-com/pydem/blob/master/pydem/dem_processing.py#L691

so, you can try: dem_proc.save_twi('filepath')

[tgmueller]

Thank you so much @mpu-creare! I am sorry but I'm just learning how to program. I received errors when I tried this in a number of different ways. The last time I tried still gave me errors: dem_proc.save_twi(twi,'E:\TerrainAnalysis\TWI_new.tif'). Please see below. What am I doing wrong?

best, Tom

[mpu-creare]

No problem. Try:

dem_proc.save_twi('E:\TerrainAnalysis')

[tgmueller]

@mpu-creare .... thanks, I tried that but am still getting an error. BTW where can I find examples in the documentation?

[mpu-creare]

Check the pydem/examples folder for examples.

I remember what the problem was with the save_twi command, sorry, this is not obvious at all. Try making the directory E:\TerrainAnalysis\twi, and then use the same code again. pydem expects that the twi directory already exists when using that function. The save_twi function is really a helper function for the processing_manager module.

[tgmueller]

@mpu-creare Thanks.... It worked. UCA also works. I have a problem when I export slope (see below). Also, I wasn't able to find the examples subdirectory on my hardrive or on GitHub. Using Google, I found this subeditor: https://pypkg.com/pypi/pydem/f/pydem/examples/. One other thing, can I assume that it is using D-infinity rather than D8. I'm looking for an example of how to create a depressionless DEM. I had previously done this with ArcGIS. Thanks again. Were you involved in developing PyDEM?

[tgmueller]

Oh Now I found the examples: https://github.com/creare-com/pydem/tree/master/pydem/examples Thanks!

[mpu-creare]

The slope saving probably has the same problem as the twi. Make sure you have the following directories created: ang mag twi uca elev and then I think you'll be all set saving out geotiffs.

In terms of who developed pydem, you can always check the insights on github: https://github.com/creare-com/pydem/graphs/contributors

[mpu-creare]

Closing due to inactivity.