bug in the classical PSHA calculator

[Vitor Silva]

A bug occurs when the attached demo is ran:

Traceback (most recent call last):
 File "bin/openquake", line 198, in <module>
   log_file=args.log_file)
 File "/Users/acerisara/Projects/MyOpenQuake/openquake/engine.py", line 815, in run_job
   _launch_job(job_ctxt, sections)
 File "/Users/acerisara/Projects/MyOpenQuake/openquake/engine.py", line 888, in _launch_job
   calculator.execute()
 File "/Users/acerisara/Projects/MyOpenQuake/openquake/calculators/risk/classical/core.py", line 240, in execute
   task.wait()
 File "/Library/Python/2.7/site-packages/celery/result.py", line 99, in wait
   return self.get(*args, **kwargs)
 File "/Library/Python/2.7/site-packages/celery/result.py", line 95, in get
   interval=interval)
 File "/Library/Python/2.7/site-packages/celery/backends/amqp.py", line 147, in wait_for
   raise self.exception_to_python(meta["result"])
django.core.exceptions.ObjectDoesNotExist: HazardCurveData matching query does not exist.

[Andrea Cerisara (acerisara)]

From what I understood looking at the data, the problem could arise when we have region borders that are not an exact multiple of the grid spacing. For example, having a region like this,

region = shapes.Region.from_coordinates([(1.0, 2.0), (2.0, 2.0), (2.0, 1.0), (1.0, 1.0)])
grid_spacing = 0.5

and iterating over it, we generate these sites (that are supposed to be the center of the cell):

[Site(1.0, 1.0), Site(1.5, 1.0), Site(2.0, 1.0), Site(1.0, 1.5), Site(1.5, 1.5), Site(2.0, 1.5), Site(1.0, 2.0), Site(1.5, 2.0), Site(2.0, 2.0)]

and these points:

Point(row=0, column=0)
Point(row=0, column=1)
Point(row=0, column=2)
Point(row=1, column=0)
Point(row=1, column=1)
Point(row=1, column=2)
Point(row=2, column=0)
Point(row=2, column=1)
Point(row=2, column=2)

In this case whatever site I use inside the region the component is able to correctly translate to the grid point, for example:

point = region.grid.point_at(shapes.Site(1.9, 2.0))
site = region.grid.site_at(point)

`point` is Point(row=2, column=2) and `site` is Site(2.0, 2.0) (the center of the cell)

Let's now take a similar region:

region = shapes.Region.from_coordinates([(1.0, 2.0), (2.3, 2.0), (2.3, 1.0), (1.0, 1.0)])
grid_spacing = 0.5

In this case the region borders are not an exact multiple of the grid spacing. We always have the same sites:

[Site(1.0, 1.0), Site(1.5, 1.0), Site(2.0, 1.0), Site(1.0, 1.5), Site(1.5, 1.5), Site(2.0, 1.5), Site(1.0, 2.0), Site(1.5, 2.0), Site(2.0, 2.0)]

and I think it's correct because the next generated points/sites would be outside the region (Site(2.5, 2.0), for example). We also have the same points:

Point(row=0, column=0)
Point(row=0, column=1)
Point(row=0, column=2)
Point(row=1, column=0)
Point(row=1, column=1)
Point(row=1, column=2)
Point(row=2, column=0)
Point(row=2, column=1)
Point(row=2, column=2)

Here is the problem. If I take, for example, Site(2.25, 2.0) that is inside the region (note: risk sites are taken from the exposure file that is defined by the user, so every site inside the region is valid, also really closed to the borders) and ask for the point, I get:

point = region.grid.point_at(shapes.Site(2.25, 2.0))
site = region.grid.site_at(point)

`point` is Point(row=2, column=3) and `site` is Site(2.5, 2.0) (the center of the cell)

Site(2.5, 2.0) is not in the set of the sites computed by the region. Since we compute the hazard by iterating over these sites this means we can't find the pre-computed hazard input for the risk site we are interested in.

A solution I see to address this problem is to add an additional row and column in the grid we generate, to be sure we handle correctly sites really closed the region border.

[Lars Butler (lars-butler)]

This should solve the problem in the short term, but I think we should step back and consider the logic surrounding Risk and Hazard geometry. This is especially important if we ever want to be able to run Risk calculations on pre-computed Hazard curves.

The problem is, in my opinion, that we're treating Risk grid cells first and foremost as points. However, since we have a 'cell size', we're really dealing with polygons (location + width + height = polygon). Thus, I think it would be much more stable to start treating risk cells as polygons; it would also remove the Region/Grid abstraction layer, which, at least for me, is difficult to wrap my head around. I think there is a simpler solution.

It would be important, however, to go back to the Risk science crew (Vitor) and make sure we really understand what it is the Classical Risk calculator is trying to do. We should definitely clarify this issue if it becomes a priority to calculate Risk on pre-computed results.

There are more details I've thought about which I haven't written here. Maybe we have a chat soon and do some brainstorming.