mdbartos / pysheds

:earth_americas: Simple and fast watershed delineation in python.
GNU General Public License v3.0
708 stars 191 forks source link

grid.catchment ..... TypeError: '<' not supported between instances of 'AxesSubplot' and 'float' #174

Closed fluviotect closed 2 years ago

fluviotect commented 2 years ago
# Specify outlet
x, y = 1824187, 5456997

# specify crs
crs = 'EPSG:2193'

# new_crs = pyproj.Proj('+init=epsg:2193') # pyproj syntax has changed
crs_proj = pyproj.CRS(crs)

# Resolve flats in DEM
inflated_dem = grid.resolve_flats(dem)

# compute flow directions...as_crs kwarg added per https://github.com/mdbartos/pysheds/issues/20#issuecomment-416018819
fdir = grid.flowdir(inflated_dem, as_crs=crs_proj)

# Delineate a catchment
catch = grid.catchment(x=x, y=y, fdir=fdir, xytype='coordinate')

TypeError Traceback (most recent call last) ~\AppData\Local\Temp/ipykernel_6164/1784417402.py in 1 # Delineate a catchment ----> 2 catch = grid.catchment(x=x, y=y, fdir=fdir, xytype='coordinate') 3 4 # Clip the view to the catchment 5 grid.clip_to(catch)

~\miniconda3\envs\CatchDelin\lib\site-packages\pysheds\sgrid.py in catchment(self, x, y, fdir, pour_value, dirmap, nodata_out, xytype, routing, snap, **kwargs) 685 xmin, ymin, xmax, ymax = fdir.bbox 686 if xytype in {'label', 'coordinate'}: --> 687 if (x < xmin) or (x > xmax) or (y < ymin) or (y > ymax): 688 raise ValueError('Pour point ({}, {}) is out of bounds for dataset with bbox {}.' 689 .format(x, y, (xmin, ymin, xmax, ymax)))

TypeError: '<' not supported between instances of 'AxesSubplot' and 'float'

Based on the plotting: image

the coordinate values seem to be within the graphical space.

mdbartos commented 2 years ago

The way that projections are handled was changed in v0.3. See subheading "Converting the Raster coordinate reference system" in the docs: https://mattbartos.com/pysheds/raster.html

The as_crs keyword argument was removed from all functions. It is now probably being swallowed by **kwargs.

For your case, consider something like this:

# Specify outlet
x, y = 1824187., 5456997.

# specify crs
crs = 'EPSG:2193'

# new_crs = pyproj.Proj('+init=epsg:2193') # pyproj syntax has changed
crs_proj = pyproj.CRS(crs)

# Resolve flats in DEM
inflated_dem = grid.resolve_flats(dem)

# Project DEM to new CRS
proj_dem = inflated_dem.to_crs(crs_proj)

# Set grid's viewfinder to projected viewfinder
grid.viewfinder = proj_dem.viewfinder

# Compute flow directions
fdir = grid.flowdir(proj_dem)

# Delineate a catchment
catch = grid.catchment(x=x, y=y, fdir=fdir, xytype='coordinate')

The error message you are getting makes me think that a matplotlib plot was passed in as one of the arguments though.

fluviotect commented 2 years ago

Cheers!

I was getting a type error (TypeError: crs must be a pyproj.Proj object.) on the reproj and had to edit one bit of the above:

crs_proj = pyproj.CRS(crs)

needs to be:

crs_proj = pyproj.Proj(crs)

your link to the docs got me sorted. Cheers

fluviotect commented 2 years ago

and good thought on the matplotlib

If I run this block of code:

import numpy as np
import matplotlib.pyplot as plt
from matplotlib import colors

fig, ax = plt.subplots(figsize=(8,6))
fig.patch.set_alpha(0)

ax = plt.subplot(1,2,1)
plt.imshow(dem, extent=grid.extent, cmap='terrain', zorder=1)
plt.colorbar(label='Elevation (m)')
plt.grid(zorder=0)
plt.title('Digital elevation map', size=14)
plt.xlabel('NZTM - Easting')
plt.ylabel('NZTM - Northing')

x = plt.subplot(1,2,2)
plt.imshow(pit_filled_dem, extent=grid.extent, cmap='terrain', zorder=1)
plt.colorbar(label='Elevation (m)')
plt.grid(zorder=0)
plt.title('Flow Direction Map', size=14)
plt.xlabel('NZTM - Easting')
plt.ylabel('NZTM - Northing') 

plt.tight_layout()

I get the error. If I skip it, then grid.catchment runs.

This is made even more strange by the fact that I don't even call fdir in the block.

Maybe it's a jupyter issue?

mdbartos commented 2 years ago

It's because x is being set to plt.subplot

fluviotect commented 2 years ago

classic...good catch