mheberger / delineator

Fast, accurate watershed delineation using hybrid vector- and raster-based methods and data from MERIT-Hydro
MIT License
84 stars 20 forks source link

Catchment shapefile: no CRS definied #3

Closed AlexDo1 closed 11 months ago

AlexDo1 commented 11 months ago

Hi Matthew!

I currently try to delineate catchments in Germany (Basin 23). I followed your instructions from the front page to download the data for Germany and set up everything in config.py. When I run python delineate.py in the terminal, I get the following output and error:

$ python delineate.py 
/home/alexander/Github/camels/camelsp/merit_hydro/delineator/venv/lib/python3.9/site-packages/geopandas/_compat.py:123: UserWarning: The Shapely GEOS version (3.10.3-CAPI-1.16.1) is incompatible with the GEOS version PyGEOS was compiled with (3.10.4-CAPI-1.16.2). Conversions between both will be slow.
  warnings.warn(
/home/alexander/Github/camels/camelsp/merit_hydro/delineator/venv/lib/python3.9/site-packages/geopandas/_compat.py:123: UserWarning: The Shapely GEOS version (3.10.3-CAPI-1.16.1) is incompatible with the GEOS version PyGEOS was compiled with (3.10.4-CAPI-1.16.2). Conversions between both will be slow.
  warnings.warn(
Reading your outlets data in: /home/alexander/Github/camels/camelsp/merit_hydro/data/outlets.csv
Finding out which Level 2 megabasin(s) your points are in
Your watershed outlets are in 1 basin(s)

Beginning delineation for 10 outlet point(s) in Level 2 Basin #23.
Reading catchment geodata in /home/alexander/Github/camels/camelsp/merit_hydro/data/shp/merit_catchments/cat_pfaf_23_MERIT_Hydro_v07_Basins_v01.shp
Traceback (most recent call last):
  File "/home/alexander/Github/camels/camelsp/merit_hydro/delineator/delineate.py", line 701, in <module>
    delineate()
  File "/home/alexander/Github/camels/camelsp/merit_hydro/delineator/delineate.py", line 449, in delineate
    catchments_gdf.to_crs(crs, inplace=True)
  File "/home/alexander/Github/camels/camelsp/merit_hydro/delineator/venv/lib/python3.9/site-packages/geopandas/geodataframe.py", line 1364, in to_crs
    geom = df.geometry.to_crs(crs=crs, epsg=epsg)
  File "/home/alexander/Github/camels/camelsp/merit_hydro/delineator/venv/lib/python3.9/site-packages/geopandas/geoseries.py", line 1124, in to_crs
    self.values.to_crs(crs=crs, epsg=epsg), index=self.index, name=self.name
  File "/home/alexander/Github/camels/camelsp/merit_hydro/delineator/venv/lib/python3.9/site-packages/geopandas/array.py", line 762, in to_crs
    raise ValueError(
ValueError: Cannot transform naive geometries.  Please set a crs on the object first.

I think this error occurs because there is no .prj file for the merit catchments on Google Drive (linked in your instruction): https://drive.google.com/drive/folders/1owkvZQBMZbvRv3V4Ff3xQPEgmAC48vJo

Your sample data for Iceland includes the file cat_pfaf_27_MERIT_Hydro_v07_Basins_v01.prj which I cannot find for my area.

There are also some incompatibility warnings at the beginning of the output but I do not think that those are causing the problem.

I tried to run python delineate.py in a fresh conda and a fresh venv environemnt, both produce the same error.

Let me know if you need anything else and thank you!

mheberger commented 11 months ago

That error is GeoPandas saying it does not know what projection the data is in. It looks like the authors at Reach Hydro failed to include a .prj file with their updated "bugfix" shapefiles. (I used the older files, pre-bug fix.) I suggest you just add your own .prj file. Just create a text file, name it same as the other shapefile components (.shp, .dbf, ...) but give it the extension ".prj". Then put the definition for the in the file:

GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]] .

By the way, if you're just doing a few watersheds, it is probably faster to use the web app. It is very fast, and lets you immediately check the results on the map. You can download the results in a few different geodata formats.

On Mon, Sep 25, 2023 at 2:41 PM Alexander Dolich @.***> wrote:

Hi Matthew!

I currently try to delineate catchments in Germany (Basin 23). I followed your instructions https://github.com/mheberger/delineator#detailed-instructions from the front page to download the data for Germany and set up everything in config.py. When I run python delineate.py in the terminal, I get the following output and error:

$ python delineate.py /home/alexander/Github/camels/camelsp/merit_hydro/delineator/venv/lib/python3.9/site-packages/geopandas/_compat.py:123: UserWarning: The Shapely GEOS version (3.10.3-CAPI-1.16.1) is incompatible with the GEOS version PyGEOS was compiled with (3.10.4-CAPI-1.16.2). Conversions between both will be slow. warnings.warn( /home/alexander/Github/camels/camelsp/merit_hydro/delineator/venv/lib/python3.9/site-packages/geopandas/_compat.py:123: UserWarning: The Shapely GEOS version (3.10.3-CAPI-1.16.1) is incompatible with the GEOS version PyGEOS was compiled with (3.10.4-CAPI-1.16.2). Conversions between both will be slow. warnings.warn( Reading your outlets data in: /home/alexander/Github/camels/camelsp/merit_hydro/data/outlets.csv Finding out which Level 2 megabasin(s) your points are in Your watershed outlets are in 1 basin(s)

Beginning delineation for 10 outlet point(s) in Level 2 Basin #23. Reading catchment geodata in /home/alexander/Github/camels/camelsp/merit_hydro/data/shp/merit_catchments/cat_pfaf_23_MERIT_Hydro_v07_Basins_v01.shp Traceback (most recent call last): File "/home/alexander/Github/camels/camelsp/merit_hydro/delineator/delineate.py", line 701, in delineate() File "/home/alexander/Github/camels/camelsp/merit_hydro/delineator/delineate.py", line 449, in delineate catchments_gdf.to_crs(crs, inplace=True) File "/home/alexander/Github/camels/camelsp/merit_hydro/delineator/venv/lib/python3.9/site-packages/geopandas/geodataframe.py", line 1364, in to_crs geom = df.geometry.to_crs(crs=crs, epsg=epsg) File "/home/alexander/Github/camels/camelsp/merit_hydro/delineator/venv/lib/python3.9/site-packages/geopandas/geoseries.py", line 1124, in to_crs self.values.to_crs(crs=crs, epsg=epsg), index=self.index, name=self.name File "/home/alexander/Github/camels/camelsp/merit_hydro/delineator/venv/lib/python3.9/site-packages/geopandas/array.py", line 762, in to_crs raise ValueError( ValueError: Cannot transform naive geometries. Please set a crs on the object first.

I think this error occurs because there is no .prj file for the merit catchments on Google Drive (linked in your instruction): https://drive.google.com/drive/folders/1owkvZQBMZbvRv3V4Ff3xQPEgmAC48vJo

Your sample data for Iceland includes the file cat_pfaf_27_MERIT_Hydro_v07_Basins_v01.prj which I cannot find for my area.

There are also some incompatibility warnings at the beginning of the output but I do not think that those are causing the problem.

I tried to run python delineate.py in a fresh conda and a fresh venv environemnt, both produce the same error.

Let me know if you need anything else and thank you!

— Reply to this email directly, view it on GitHub https://github.com/mheberger/delineator/issues/3, or unsubscribe https://github.com/notifications/unsubscribe-auth/ADYEMR6YX2Z4E7QU6EWPKBTX4F3WRANCNFSM6AAAAAA5GA4XMU . You are receiving this because you are subscribed to this thread.Message ID: @.***>

AlexDo1 commented 11 months ago

Thank you very much, creating the .prj file manually did work.

But I also downloaded the "old" shapefiles without the bugfixes, the .prj file is not there: https://drive.google.com/drive/folders/1PxHZB-vO5O7EUQJSxYS3ptV-tuI39uKD So in theory, everyone who follows your (great) Instructions on the Github front page should stumble upon this error.

I tried the delineation first with the bugfix version, but delineate.py failed, as the catchment file must have the name of the "old" version, without bugfixes: https://github.com/mheberger/delineator/blob/ae8575bf18dc49d109fbf77cbc2d740d0edd9585/delineate.py#L448

mheberger commented 11 months ago

Good to know, thanks. I must have added those files myself and forgotten about it. I could add a line to the code to set the projection, since it is known.

mheberger commented 11 months ago

I see what the issue was. Line 455 was not doing what I thought it was. I need to change it from:

catchments_gdf.to_crs(crs, inplace=True)

to

catchments_gdf.set_crs(crs, inplace=True)

AlexDo1 commented 11 months ago

Thank you, that fixed the issue!