schism-dev / pyschism

Python interface for handling the SCHISM model.
https://schism-dev.github.io/schism/master/getting-started/pre-processing-with-pyschism/overview.html
Apache License 2.0
25 stars 21 forks source link

Errors reading hgrid.gr3 #99

Open calquigs opened 1 year ago

calquigs commented 1 year ago

I am new to SCHISM and pyschism, so apologies if there are some basic mistakes I am making here:

I am trying to read in and plot an hgrid.gr3 file, but I get the following errors:

hgrid = Hgrid.open('Desktop/hgrid.gr3')

File "/<home>/miniforge3/envs/opendrift/lib/python3.11/site-packages/requests/models.py", line 439, in prepare_url
    raise MissingSchema(
requests.exceptions.MissingSchema: Invalid URL 'Desktop/hgrid.gr3': No scheme supplied. Perhaps you meant https://Desktop/hgrid.gr3?

During handling of the above exception, another exception occurred:

File "/<home>/miniforge3/envs/opendrift/lib/python3.11/site-packages/pyschism/mesh/parsers/grd.py", line 64, in buffer_to_dict
    npts, ibtype = map(int, buf.readline().split()[:2])
    ^^^^^^^^^^^^
ValueError: invalid literal for int() with base 10: '='

Any idea where I may be going wrong?

calquigs commented 1 year ago

Update, I went into pyschism/mesh/parsers/grd.py to add some debugging print statements, and I still get the same error message with the same line numbers, but the line numbers are now referring to a diferrent part of the code because I have added new lines above. So the error message now reads:

File "/<home>/miniforge3/envs/opendrift/lib/python3.11/site-packages/pyschism/mesh/parsers/grd.py", line 64, in buffer_to_dict
    NBOU = int(buf.readline().split()[0])
        ^^^^^^^^^^^^
ValueError: invalid literal for int() with base 10: '='

A bit confusing.

cuill commented 1 year ago

This is something wrong with the boundary information. Is it possible to share the hgrid.gr3 file?

calquigs commented 1 year ago

Yes. The file is over 25mb so I wasn't able to attach it here, but it is available at https://github.com/calquigs/superreefs/blob/master/hgrid.gr3.

cuill commented 1 year ago

The hgrid is missing information of exterior land boundary/island for land boundaries:

16 = number of land boundaries (including islands)
1743 = Total number of land boundary nodes
753 0 = Number of nodes for land boundary 1 ('0' means the exterior land boundary)
30381 ! first node on this segment
.......................................
1 !last node on this segment
741 0 ! Number of nodes for land boundary 2 ('0' means the exterior boundary)
.
.
.
10 1 = Number of nodes for island boundary 1 ('1' means island)
29448 ! first node on this island
.
.
.
29449 !last node on this island (note this is different from the first node ‘29448’ above)

See hgrid.gr3 example here: https://schism-dev.github.io/schism/master/input-output/hgrid.html

calquigs commented 1 year ago

My colleague who is developing the model took a look at your notes and made some edits (new file here) to "close" the island boundaries, however I'm still getting a similar error:

>>> hgrid = Hgrid.open('Desktop/hgrid_islfix.gr3')
hgrid_islfix.gr3
<_io.TextIOWrapper name='Desktop/hgrid_islfix.gr3' mode='r' encoding='UTF-8'>
Traceback (most recent call last):
  File "/<env>/lib/python3.11/site-packages/pyschism/mesh/hgrid.py", line 31, in open
    response = requests.get(path)
               ^^^^^^^^^^^^^^^^^^
  File "/<env>/lib/python3.11/site-packages/requests/api.py", line 73, in get
    return request("get", url, params=params, **kwargs)
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/<env>/lib/python3.11/site-packages/requests/api.py", line 59, in request
    return session.request(method=method, url=url, **kwargs)
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/<env>/lib/python3.11/site-packages/requests/sessions.py", line 575, in request
    prep = self.prepare_request(req)
           ^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/<env>/lib/python3.11/site-packages/requests/sessions.py", line 486, in prepare_request
    p.prepare(
  File "/<env>/lib/python3.11/site-packages/requests/models.py", line 368, in prepare
    self.prepare_url(url, params)
  File "/<env>/lib/python3.11/site-packages/requests/models.py", line 439, in prepare_url
    raise MissingSchema(
requests.exceptions.MissingSchema: Invalid URL 'Desktop/hgrid_islfix.gr3': No scheme supplied. Perhaps you meant https://Desktop/hgrid_islfix.gr3?

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/<env>/lib/python3.11/site-packages/pyschism/mesh/hgrid.py", line 38, in open
    _grd = grd.read(path, crs=crs)
           ^^^^^^^^^^^^^^^^^^^^^^^
  File "/<env>/lib/python3.11/site-packages/pyschism/mesh/parsers/grd.py", line 183, in read
    grd = buffer_to_dict(stream)
          ^^^^^^^^^^^^^^^^^^^^^^
  File "/<env>/lib/python3.11/site-packages/pyschism/mesh/parsers/grd.py", line 68, in buffer_to_dict
    npts, ibtype = map(int, buf.readline().split()[:2])
    ^^^^^^^^^^^^
ValueError: invalid literal for int() with base 10: '='

(posting w/ full traceback this time, don't know if that helps at all)

Thanks for your help with this!

cuill commented 1 year ago

@calquigs This line npts, ibtype = map(int, buf.readline().split()[:2]) takes header information for each land or island boundary, for example, line 848638 in your hgrid.gr3, and assigns the first two strings (separated by space) to npts (number of nodes), ibtype (boundary type). In your hgrid, line 848638 is: 333 = Number of nodes for open boundary 1 which is missing boundary type (0/1) information.

The correct format should be: 333 0 = Number of nodes for open boundary 1

calquigs commented 1 year ago

Ok great! I've gone through and in lines specifying Number of nodes for open/land boundary I've added 0 for open, 1 for land.

I'm still getting the exception requests.exceptions.MissingSchema: Invalid URL 'hgrid_islfix.gr3': No scheme supplied. Perhaps you meant https://hgrid_islfix.gr3? which I'm not sure what that is referring to or if I should be worried about it? And below that I now get:

Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/<env>/lib/python3.11/site-packages/pyschism/mesh/hgrid.py", line 38, in open
    _grd = grd.read(path, crs=crs)
           ^^^^^^^^^^^^^^^^^^^^^^^
  File "/<env>/lib/python3.11/site-packages/pyschism/mesh/parsers/grd.py", line 183, in read
    grd = buffer_to_dict(stream)
          ^^^^^^^^^^^^^^^^^^^^^^
  File "/<env>/lib/python3.11/site-packages/pyschism/mesh/parsers/grd.py", line 68, in buffer_to_dict
    npts, ibtype = map(int, buf.readline().split()[:2])
    ^^^^^^^^^^^^
ValueError: not enough values to unpack (expected 2, got 1)

Thanks again for your help and patience. For errors caused by input file formatting, it would be helpful if the error message had a way of showing the line number in the file where the error occurred, that way I would have a better chance of being able to debug by myself!

cuill commented 1 year ago

@calquigs In your hgrid.gr3, no need to add an extra node to make it closed, which caused inconsistent between the number of nodes read from header information and the actual lines:

2898 0 = Number of nodes for land boundary
1
2
....
150
67
22
21
1
1
5061 0 = Number of nodes for land boundary 2

For example, for the first land boudanry in your hgrid.gr3, you don't need to add an extra node 1 in the end.