Time series extraction, output points outside the mesh

EDIT:

Recap of the discussion for this issue

Some of our points are ignored, because they lie outside the domain.

I'm proposing a new feature to enable this capacity. There are 2 ways of implementing this:

1 - Correct the existing NN (nearest neighbour) implementation in out_history.f90:

• the INTERP_METHOD is only accessible from the out_history.f90 and should be a recognized keyword in the TELEMAC Cas file. @alexander.breugem , from the way it has been implemented, I guess it was your first idea.
• the INTERP_METHOD = 0 does not work for locations outside the domain mesh. IMO it should with any kind of point location.

2 - change the way we read the time series coordinates file

• This method (a bit easier to implement) consists in reading directly the nodes ids in the time series coordinates files. Columns should be: node_id ikle_id name, instead of x y ikle_id name.

I provided a modified version of out_history.f90 that works for option (2) but that is not well implemented at all, because it's using INTERP_METHOD = 0, and it is mixing the 2 different concepts I listed above.

I think I would rather need your feedback on how you picture the implementation of this new feature (and what is the best option (1/2 or both?).

---

Original message

The issue arises when using time series output file.

The implementation of time export is recent in TELEMAC, and work only with spherical coordinates (ESPG 3857) at the moment. (I opened an issue on that matter)

For some reason, during the execution of the fortran code, some observation points get ignored because

 WARNING!            2  OUTPUT POINTS ARE OUTSIDE THE MODEL DOMAIN.
 NO HISTORY FILE WILL BE GENERATED FOR THESE POINTS.

Whereas I made sure the points where right on the closest nodes. It must be because because of some rounding error in the transformation (from WGS84 to mercator).

Looking at the code:

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!                         INTERPOLATION ROUTINES
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

!                   *******************************
                    SUBROUTINE INTERP_POINT_PREPARE &
!                   *******************************
!
     (XINT,YINT,NPINT,INTVAR,INTERP_METHOD,MESH2D)
!
!***********************************************************************
! BIEF
!***********************************************************************
!
!>@brief    Determines weight factors for interpolation
!
!           Note: May be slow for large number of output points
!
!history WA BREUGEM (IMDC)
!+        1/09/2020
!+
!+
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!>@param[in]      XINT   X coordinates to interpolate
!>@param[in]      YINT   Y coordinates to interpolate
!>@param[in]      NPINT  Number of points to interpolate
!>@param[in]      XINT   Structure with coordinates to use
!>@param[in,out]  INTVAR Structure with data for interpolation
!>@param[in]      INTERP_METHOD 1 if linear interpolation, 0 if nearest
!                               neighbour interpolation 
!>@param[in]      MESH2D Structure with the 2d mesh
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
...
! Loop over all points that need to interpolated

...
          CALL IN_TRIANGLE(XP,YP,XTRI,YTRI,INTRI)

          IF(INTRI) THEN
            ! Store points
...

            IF(INTERP_METHOD.EQ.1) THEN
!             Linear interpolation
              INTVAR%WEIGHTS (I,1) = W1/W
              INTVAR%WEIGHTS (I,2) = W2/W
              INTVAR%WEIGHTS (I,3) = W3/W
            ELSE
              ! Nearest neighbour
              IF((W1.GE.W2).AND.(W1.GE.W3)) THEN
                INTVAR%WEIGHTS (I,1) = 1.D0
              ENDIF
              IF((W2.GE.W1).AND.(W2.GE.W3)) THEN
                INTVAR%WEIGHTS (I,2) = 1.D0
              ENDIF
              IF((W3.GE.W2).AND.(W3.GE.W1)) THEN
                INTVAR%WEIGHTS (I,3) = 1.D0
              ENDIF
            ENDIF
            ! Continue to next point of the list
            EXIT
          ENDIF
        ENDDO
      ENDDO

I tried the interpolation 'nearest neightbour' INTERP_METHOD=0.

It also failed because the method IN_TRIANGLE does not find the points inside the triangles. I think there are 2 problems here:

• because of the rounding error I mention above (transformation (from WGS84 to mercator), I think the method IN_TRIANGLE should be a little more more forgiving for points really close to the outskirts of the mesh.
• I don't think IN_TRIANGLE should be use for the Nearest Neighbour method. This method should assign the nearest node to point, whether the point is in or out the mesh.