UMEP-dev / UMEP

Urban Multi-scale Environmental Predictor
https://umep-docs.readthedocs.io/
61 stars 15 forks source link

URock fails after generating invalid geometry in step "Creates the 2D Röckle zones" #589

Closed jlegewie closed 7 months ago

jlegewie commented 7 months ago

Describe the bug During the step "Creates the 2D Röckle zones", URock generates an invalid geometry and fails with the error message below. Note that all the input geometries are valid. I think it happens during the part labeled "Identify vegetation zones being in building wake zones" of the function vegetationZones in Zones.py. I am including the building and vegetation vector files. I tried to crop the vector files but the error didn't happen when I restrict the file to the part in the upper right corner.

This error happens in almost half (24 out of 58) of the areas for which I am running URock. I am just focussing on one as an example here.

Vector files vectors.zip

To Reproduce Steps to reproduce the behavior:

  1. Open the linked building and vegetation vector files
  2. Open URock v2023a
  3. Run with following settings:
    {
    "inputs": {
    "ATTENUATION_FIELD": null,
    "BUILDINGS": LOCATION OF BUILDING VECTOR FILE,
    "HEIGHT_FIELD_BUILD": "HEIGHTROOF",
    "HORIZONTAL_RESOLUTION": 10,
    "INPUT_PROFILE_FILE": null,
    "INPUT_PROFILE_TYPE": 0,
    "INPUT_WIND_DIRECTION": 315.0,
    "INPUT_WIND_HEIGHT": 10.0,
    "INPUT_WIND_SPEED": 10.8,
    "LOAD_OUTPUT": false,
    "OUTPUT_FILENAME": "urock_output",
    "RASTER_OUTPUT": null,
    "SAVE_NETCDF": false,
    "SAVE_RASTER": true,
    "SAVE_VECTOR": true,
    "UROCK_OUTPUT": "TEMPORARY_OUTPUT",
    "VEGETATION": LOCATION OF VEGETATION VECTOR FILE,
    "VEGETATION_CROWN_BASE_HEIGHT": null,
    "VEGETATION_CROWN_TOP_HEIGHT": "height",
    "VERTICAL_RESOLUTION": 10,
    "WIND_HEIGHT": "1.5"
    }
    }
  4. See error

Desktop (please complete the following information):

Screenshots The first screenshot shows the invalid geometry (MULTIPOLYGON) referenced in the error message below. The second screenshot shows a larger area with an arrow pointing to the invalid geometry (although not the full vector files).

Screenshot 2024-02-22 at 10 46 00 AM Screenshot 2024-02-22 at 10 46 19 AM

Data If applicable, add data files to help us reproduce your problem.

Error meesage

----------------
Inputs
----------------

ATTENUATION_FIELD:  
BUILDINGS:  tile-24.shp
HEIGHT_FIELD_BUILD: HEIGHTROOF
HORIZONTAL_RESOLUTION:  10
INPUT_PROFILE_FILE: 
INPUT_PROFILE_TYPE: 0
INPUT_WIND_DIRECTION:   315
INPUT_WIND_HEIGHT:  10
INPUT_WIND_SPEED:   10.8
LOAD_OUTPUT:    false
OUTPUT_FILENAME:    urock_output
RASTER_OUTPUT:  template-raster-10m.tif
SAVE_NETCDF:    false
SAVE_RASTER:    true
SAVE_VECTOR:    true
UROCK_OUTPUT:   tile-24
VEGETATION: tile-24.shp
VEGETATION_CROWN_BASE_HEIGHT:   
VEGETATION_CROWN_TOP_HEIGHT:    height
VERTICAL_RESOLUTION:    10
WIND_HEIGHT:    2

Writing settings for this model run to specified output folder (Filename: RunInfoURock_YYYY_DOY_HHMM.txt)
Initiating algorithm
Creates an H2GIS Instance and load data
Creates the stacked blocks used as obstacles
Rotates obstacles to the right direction and calculates geometry properties
Creates the 2D Röckle zones
Traceback (most recent call last):
  File "JdbcPreparedStatement.java", line 237, in org.h2.jdbc.JdbcPreparedStatement.execute
org.locationtech.jts.geom.org.locationtech.jts.geom.TopologyException: org.locationtech.jts.geom.TopologyException: found non-noded intersection between LINESTRING ( 296800.41809878143 49923.29916504909, 296800.41809878143 49923.29916504906 ) and LINESTRING ( 296800.41809878143 49923.29916504909, 296800.41809878143 49923.299165049066 ) [ (296800.41809878143, 49923.29916504909, NaN) ]

The above exception was the direct cause of the following exception:

Traceback (most recent call last):
  File "JdbcPreparedStatement.java", line 237, in org.h2.jdbc.JdbcPreparedStatement.execute
Exception: Java Exception

The above exception was the direct cause of the following exception:

Traceback (most recent call last):
  File "/n/home11/jlegewie/.local/lib/python3.10/site-packages/jaydebeapi/__init__.py", line 534, in execute
    is_rs = self._prep.execute()
org.h2.jdbc.org.h2.jdbc.JdbcSQLNonTransientException: org.h2.jdbc.JdbcSQLNonTransientException: Exception calling user-defined function: "union(MULTIPOLYGON (((296801.59327265644 49922.12399117407, 296800.17905909405 49920.709777611715, 296798.76484553167 49922.1239911741, 296800.17905909405 49923.538204736455, 296801.59327265644 49922.12399117407)), ((296801.59327265644 49922.12399117407, 296800.17905909405 49920.709777611715, 296798.76484553167 49922.1239911741, 296800.17905909405 49923.538204736455, 296801.59327265644 49922.12399117407)), ((296801.59327265644 49922.12399117407, 296800.3024651866 49920.833183704315, 296800.09086883924 49923.45001448166, 296800.17905909405 49923.538204736455, 296801.59327265644 49922.12399117407)), ((296801.59327265644 49922.12399117407, 296800.17905909405 49920.709777611715, 296798.76484553167 49922.1239911741, 296800.17905909405 49923.538204736455, 296801.59327265644 49922.12399117407)), ((296800.8329505903 49921.363669107974, 296800.4210386577 49920.95175717536, 296800.3913030351 49923.325960795446, 296800.41809878143 49923.29916504909, 296800.71 49921.94, 296800.8329505903 49921.363669107974)), ((296800.8329505903 49921.363669107974, 296800.4210386577 49920.95175717536, 296800.3913030351 49923.325960795446, 296800.41809878143 49923.29916504909, 296800.71 49921.94, 296800.8329505903 49921.363669107974)), ((296800.8329505903 49921.363669107974, 296800.4210386577 49920.95175717536, 296800.3913030351 49923.325960795446, 296800.41809878143 49923.29916504909, 296800.71 49921.94, 296800.8329505903 49921.363669107974)))): found non-noded intersection between LINESTRING ( 296800.41809878143 49923.29916504909, 296800.41809878143 49923.29916504906 ) and LINESTRING ( 296800.41809878143 49923.29916504909, 296800.41809878143 49923.299165049066 ) [ (296800.41809878143, 49923.29916504909, NaN) ]"; SQL statement:
(SELECT    ST_UNION(ST_ACCUM(THE_GEOM)) AS THE_GEOM,
                                        MIN(MIN_HEIGHT) AS MIN_HEIGHT,
                                        MIN(MAX_HEIGHT) AS MAX_HEIGHT,
                                        MIN(ATTENUATIO) AS ATTENUATIO,
                                        ID_VEG
                            FROM temporary_built_vegetation_20240221220500
                            GROUP BY ID_VEG) [90105-200]
j3r3m1 commented 7 months ago

Thanks for the issue. It is under investigation.

A comment out of the error scope: I suppose you know that but in case other persons read it better remind it:

j3r3m1 commented 7 months ago

OK it should now be solved with the new UMEP version (I have updated the H2GIS version to the newest which seems to solve the geometry union there was before). The code modifications are diff here.

Please try the latest UMEP-Processing SNAPSHOT downloading the zip file here and importing it into QGIS plugin.

However, this fix is not sufficient for your case. The complex vegetation was also creating an issue that was also met by other users (cf. https://github.com/UMEP-dev/UMEP-processing/issues/63). This is now solved by a new version of the code.

Last (hopefully), the resolution used lead to new errors in the code. I am currently investigating the problem in more detail. The log is the following:

** jaydebeapi.DatabaseError: org.h2.jdbc.JdbcSQLDataException: Division par zéro: "-9.0 Division by zero: "-9.0"; SQL statement: CREATE TABLE STREET_CANYON_INIT_POINTS AS SELECT b.ID_POINT, (SIN(2(a.THETA_WIND-PI()/2))(0.5+(a.Y_WALL-b.Y_POINT) (a.D_0S-(a.Y_WALL-b.Y_POINT))/ (0.5POWER(a.D_0S,2)))) (1+(0.6a.CANYON_DELTAH)/(a.Hc+0.6a.CANYON_DELTAH)) (POWER(a.D_0S/a.Hc,2)/(1+POWER(a.D_0S/a.Hc,2))) AS U, (1-POWER(COS(a.THETA_WIND-PI()/2),2)(1+(a.Y_WALL-b.Y_POINT) (a.D_0S-(a.Y_WALL-b.Y_POINT))/(POWER(0.5a.D_0S,2)))) (1+(0.6a.CANYON_DELTAH)/(a.Hc+0.6a.CANYON_DELTAH)) (POWER(a.D_0S/a.Hc,2)/(1+POWER(a.D_0S/a.Hc,2))) AS V, (-ABS(0.5(1-(a.Y_WALL-b.Y_POINT)/(0.5a.D_0S))) (1-(a.D_0S-(a.Y_WALL-b.Y_POINT))/(0.5a.D_0S))) (1+(0.6a.CANYON_DELTAH)/(a.Hc+0.6a.CANYON_DELTAH)) (POWER(a.D_0S/a.Hc,2)/(1+POWER(a.D_0S/a.Hc,2))) AS W, a.ID_CANYON, a.UPSTREAM_HEIGHT, a.BASE_HEIGHT, a.ID_X, CAST(a.Y_WALL AS INTEGER) AS Y_WALL, a.UPWIND_FACADE_ID, a.Hc, a.ID_STACKED_BLOCK, b.ID_Y, a.DOWNWIND_FACADE_ID, b.Y_POINT, a.D_0S, a.UPSTREAM_HEIGHT*SQRT(1-POWER((a.Y_WALL-b.Y_POINT)/a.D_0C, 2)) AS Z_UP_THRESH FROM ZONE_LIMITS_STREET_CANYON_20240226160057 AS a RIGHT JOIN TEMPO_STREET_CANYON_INIT_POINTS_20240226160057 AS b ON a.ID_CANYON = b.ID_CANYON AND a.ID_X = b.ID_X [22012-224] at org.h2.message.DbException.getJdbcSQLException(DbException.java:518) at org.h2.message.DbException.getJdbcSQLException(DbException.java:489) at org.h2.message.DbException.get(DbException.java:223) at org.h2.message.DbException.get(DbException.java:199) at org.h2.value.ValueNumeric.divide(ValueNumeric.java:105) at org.h2.expression.BinaryOperation.getValue(BinaryOperation.java:123) at org.h2.expression.BinaryOperation.getValue(BinaryOperation.java:99) at org.h2.expression.BinaryOperation.getValue(BinaryOperation.java:99) at org.h2.expression.BinaryOperation.getValue(BinaryOperation.java:98) at org.h2.expression.Alias.getValue(Alias.java:37) at org.h2.command.query.Select$LazyResultQueryFlat.fetchNextRow(Select.java:1851) at org.h2.result.LazyResult.hasNext(LazyResult.java:78) at org.h2.result.FetchedResult.next(FetchedResult.java:34) at org.h2.command.query.Select.queryFlat(Select.java:728) at org.h2.command.query.Select.queryWithoutCache(Select.java:833) at org.h2.command.query.Query.queryWithoutCacheLazyCheck(Query.java:197) at org.h2.command.query.Query.query(Query.java:520) at org.h2.command.dml.Insert.insertRows(Insert.java:197) at org.h2.command.dml.Insert.update(Insert.java:135) at org.h2.command.dml.DataChangeStatement.update(DataChangeStatement.java:74) at org.h2.command.ddl.CreateTable.insertAsData(CreateTable.java:199) at org.h2.command.ddl.CreateTable.update(CreateTable.java:157) at org.h2.command.CommandContainer.update(CommandContainer.java:169) at org.h2.command.Command.executeUpdate(Command.java:256) at org.h2.command.CommandList.update(CommandList.java:65) at org.h2.command.CommandList.executeRemaining(CommandList.java:58) at org.h2.command.CommandList.update(CommandList.java:66) at org.h2.command.CommandList.executeRemaining(CommandList.java:58) at org.h2.command.CommandList.update(CommandList.java:66) at org.h2.command.CommandList.executeRemaining(CommandList.java:58) at org.h2.command.CommandList.update(CommandList.java:66) at org.h2.command.CommandList.executeRemaining(CommandList.java:58) at org.h2.command.CommandList.update(CommandList.java:66) at org.h2.command.CommandList.executeRemaining(CommandList.java:58) at org.h2.command.CommandList.update(CommandList.java:66) at org.h2.command.CommandList.executeRemaining(CommandList.java:58) at org.h2.command.CommandList.update(CommandList.java:66) at org.h2.command.CommandList.executeRemaining(CommandList.java:58) at org.h2.command.CommandList.update(CommandList.java:66) at org.h2.command.CommandList.executeRemaining(CommandList.java:58) at org.h2.command.CommandList.update(CommandList.java:66) at org.h2.command.CommandList.executeRemaining(CommandList.java:58) at org.h2.command.CommandList.update(CommandList.java:66) at org.h2.command.CommandList.executeRemaining(CommandList.java:58) at org.h2.command.CommandList.update(CommandList.java:66) at org.h2.command.CommandList.executeRemaining(CommandList.java:58) at org.h2.command.CommandList.update(CommandList.java:66) at org.h2.command.CommandList.executeRemaining(CommandList.java:58) at org.h2.command.CommandList.update(CommandList.java:66) at org.h2.command.CommandList.executeRemaining(CommandList.java:58) at org.h2.command.CommandList.update(CommandList.java:66) at org.h2.command.CommandList.executeRemaining(CommandList.java:58) at org.h2.command.CommandList.update(CommandList.java:66) at org.h2.command.CommandList.executeRemaining(CommandList.java:58) at org.h2.command.CommandList.update(CommandList.java:66) at org.h2.command.CommandList.executeRemaining(CommandList.java:58) at org.h2.command.CommandList.update(CommandList.java:66) at org.h2.command.CommandList.executeRemaining(CommandList.java:58) at org.h2.command.CommandList.update(CommandList.java:66) at org.h2.command.CommandList.executeRemaining(CommandList.java:58) at org.h2.command.CommandList.update(CommandList.java:66) at org.h2.command.CommandList.executeRemaining(CommandList.java:58) at org.h2.command.CommandList.update(CommandList.java:66) at org.h2.command.CommandList.executeRemaining(CommandList.java:58) at org.h2.command.CommandList.update(CommandList.java:66) at org.h2.command.CommandList.executeRemaining(CommandList.java:58) at org.h2.command.CommandList.update(CommandList.java:66) at org.h2.command.CommandList.executeRemaining(CommandList.java:58) at org.h2.command.CommandList.update(CommandList.java:66) at org.h2.command.CommandList.executeRemaining(CommandList.java:58) at org.h2.command.CommandList.update(CommandList.java:66) at org.h2.command.CommandList.executeRemaining(CommandList.java:58) at org.h2.command.CommandList.update(CommandList.java:66) at org.h2.command.CommandList.executeRemaining(CommandList.java:58) at org.h2.command.CommandList.update(CommandList.java:66) at org.h2.command.CommandList.executeRemaining(CommandList.java:58) at org.h2.command.CommandList.update(CommandList.java:66) at org.h2.command.CommandList.executeRemaining(CommandList.java:58) at org.h2.command.CommandList.update(CommandList.java:66) at org.h2.command.CommandList.executeRemaining(CommandList.java:58) at org.h2.command.CommandList.update(CommandList.java:66) at org.h2.command.CommandList.executeRemaining(CommandList.java:58) at org.h2.command.CommandList.update(CommandList.java:66) at org.h2.command.Command.executeUpdate(Command.java:256) at org.h2.server.TcpServerThread.process(TcpServerThread.java:413) at org.h2.server.TcpServerThread.run(TcpServerThread.java:191) at java.base/java.lang.Thread.run(Thread.java:834)

j3r3m1 commented 7 months ago

The problem was actually an error in the equation proposed in the issue https://github.com/UMEP-dev/UMEP-processing/issues/45.

The last version of the code also fix this. If the initial error is solved, please close. Open a new one if a new one raises. Calculations are very long for this case, it did not reach the end of my computer so I cannot say...

jlegewie commented 7 months ago

Thanks! I will try the new version sometime this week and open a new issue if necessary.

And thanks for the context as well. We are also planing to test SimScale's Pedestrian Wind Comfort Simulation (or some other CFD), which is supposed to work with high-rise buildings. 10m resolution is just a first try with the idea to increase resolution later. In any case, I will report back but it will all take time.