Top 5 Recent Articles

ARTICLES CATEGORIES

Computing Cardinal Directions to nearby geometries

I was asked the other day if I could help derive the direction from a selected land parcel to its neighbouring parcels (request slightly modified):

Below is a screenshot showing an example of what I mean. I need to get the directions of parcel no (12) , it should give me (15) as north, (14) as south, (1,8) as west.

Function

I have chosen to implement the direction from one point to another in a generic function that can accept and process geodetic (long/lat) data and planar data. The function does not determine if a supplied SDO_GEOMETRY p_start_point’s sdo_srid is geodetic or not (this can, however, be done easily), rather it exposes a parameter called p_projected which the user must set to 1 (or any positive number) if the SDO_SRID is planar and 0 if geodetic (eg 8307). Failure to set the right value will generate strange results.

In addition, the created function computes a bearing from p_start_point to p_end_point and then turns this into one of 16 cardinal directions based on the compass rose. If you only need the 4 main points, modify the function accordingly.

 CREATE OR REPLACE
 FUNCTION CardinalDirection(p_start_point IN sdo_geometry,
                            p_end_point   IN sdo_geometry,
                            p_tolerance   IN NUMBER DEFAULT 0.05,
                            p_projected   IN INTEGER DEFAULT 1)
 RETURN varchar2 deterministic
 AS
   c_PI CONSTANT NUMBER := acos(-1);
   v_bearing NUMBER;
   v_tilt    NUMBER;
   v_vertex1 mdsys.vertex_type;
   v_vertex2 mdsys.vertex_type;
   FUNCTION ST_Azimuth(p_dE1 IN NUMBER, p_dN1 IN NUMBER,
                       p_dE2 IN NUMBER, p_dN2 IN NUMBER)
   RETURN NUMBER
   IS
       dBearing NUMBER;
       dEast    NUMBER;
       dNorth   NUMBER;
   BEGIN
       IF (p_dE1 IS NULL OR p_dN1 IS NULL OR
           p_dE2 IS NULL OR p_dE2 IS NULL ) THEN
          RETURN NULL;
       END IF;
       IF ( (p_dE1 = p_dE2) AND (p_dN1 = p_dN2) ) THEN
          RETURN NULL;
       END IF;
       dEast  := p_dE2 - p_dE1;
       dNorth := p_dN2 - p_dN1;
       IF ( dEast = 0 ) THEN
           dBearing := CASE WHEN ( dNorth < 0 )
                            THEN c_PI
                            ELSE 0
                        END;
       ELSE
           dBearing := -aTan(dNorth / dEast) + c_PI / 2;
       END IF;
       RETURN CASE WHEN ( dEast < 0 )
                   THEN dBearing + c_PI
                   ELSE dBearing
               END;
   END ST_Azimuth;
 BEGIN
   IF (p_start_point IS NULL OR p_end_point IS NULL) THEN
      RETURN NULL;
   END IF;
   -- If Geodetic....
   IF (NVL(ABS(p_projected),1) = 0 ) THEN
     SDO_UTIL.BEARING_TILT_FOR_POINTS(
              p_start_point,
              p_end_point,
              p_tolerance,v_bearing,v_tilt);
   ELSE
     v_vertex1 := sdo_util.getVertices(p_start_point)(1);
     v_vertex2 := sdo_util.getVertices(p_end_point)(1);
     v_bearing := ST_Azimuth(v_vertex1.X, v_vertex1.Y,
                             v_vertex2.X, v_vertex2.Y);
   END IF;
   v_bearing := round(v_bearing * (180.0 / acos(-1)),2);
   RETURN CASE
          WHEN v_bearing IS NULL THEN 'NULL'
          WHEN /* Cardinal Point 0     */ v_bearing BETWEEN 348.75 AND   0.0
                                       OR v_bearing BETWEEN   0.0  AND  11.25 THEN 'N'
          WHEN /* Cardinal Point 22.5  */ v_bearing BETWEEN  11.25 AND  33.75 THEN 'NNE'
          WHEN /* Cardinal Point 45    */ v_bearing BETWEEN  33.75 AND  56.25 THEN 'NE'
          WHEN /* Cardinal Point 67.5  */ v_bearing BETWEEN  56.25 AND  78.75 THEN 'ENE'
          WHEN /* Cardinal Point 90    */ v_bearing BETWEEN  78.75 AND 101.25 THEN 'E'
          WHEN /* Cardinal Point 112.5 */ v_bearing BETWEEN 101.25 AND 123.75 THEN 'ESE'
          WHEN /* Cardinal Point 135   */ v_bearing BETWEEN 123.75 AND 146.25 THEN 'SE'
          WHEN /* Cardinal Point 157.5 */ v_bearing BETWEEN 146.25 AND 168.75 THEN 'SSE'
          WHEN /* Cardinal Point 180   */ v_bearing BETWEEN 168.75 AND 191.25 THEN 'S'
          WHEN /* Cardinal Point 202.5 */ v_bearing BETWEEN 191.25 AND 213.75 THEN 'SSW'
          WHEN /* Cardinal Point 225   */ v_bearing BETWEEN 213.75 AND 236.25 THEN 'SW'
          WHEN /* Cardinal Point 247.5 */ v_bearing BETWEEN 236.25 AND 258.75 THEN 'WSW'
          WHEN /* Cardinal Point 270   */ v_bearing BETWEEN 258.75 AND 281.25 THEN 'W'
          WHEN /* Cardinal Point 292.5 */ v_bearing BETWEEN 281.25 AND 303.75 THEN 'WNW'
          WHEN /* Cardinal Point 315   */ v_bearing BETWEEN 303.75 AND 326.25 THEN 'NW'
          WHEN /* Cardinal Point 337.5 */ v_bearing BETWEEN 326.25 AND 348.75 THEN 'NNW'
          ELSE to_char(v_bearing)
        END;
 END CardinalDirection;
 /
 SHOW errors

To use the function, we can execute the following.

 DROP   TABLE ADJOINING_PARCELS;
 CREATE TABLE ADJOINING_PARCELS (
   objectid INTEGER,
   shape    sdo_geometry);
 INSERT INTO ADJOINING_PARCELS (OBJECTID,SHAPE) VALUES (1,SDO_GEOMETRY(2003,82362,NULL,SDO_ELEM_INFO_ARRAY(1,1003,1),SDO_ORDINATE_ARRAY(359154.918,2892270.678,359132.697,2892276.612,359129.913,2892266.19,359125.472,2892249.561,359147.693,2892243.626,359153.37,2892264.881,359154.918,2892270.678)));
 INSERT INTO ADJOINING_PARCELS (OBJECTID,SHAPE) VALUES (4,SDO_GEOMETRY(2003,82362,NULL,SDO_ELEM_INFO_ARRAY(1,1003,1),SDO_ORDINATE_ARRAY(359170.4,2892328.646,359148.179,2892334.581,359145.869,2892325.931,359140.438,2892305.596,359162.659,2892299.662,359164.723,2892307.391,359170.4,2892328.646)));
 INSERT INTO ADJOINING_PARCELS (OBJECTID,SHAPE) VALUES (8,SDO_GEOMETRY(2003,82362,NULL,SDO_ELEM_INFO_ARRAY(1,1003,1),SDO_ORDINATE_ARRAY(359162.659,2892299.662,359140.438,2892305.596,359132.697,2892276.612,359154.918,2892270.678,359159.047,2892286.136,359162.659,2892299.662)));
 INSERT INTO ADJOINING_PARCELS (OBJECTID,SHAPE) VALUES (10,SDO_GEOMETRY(2003,82362,NULL,SDO_ELEM_INFO_ARRAY(1,1003,1),SDO_ORDINATE_ARRAY(359200.35,2892320.647,359170.4,2892328.646,359164.723,2892307.391,359194.674,2892299.392,359200.35,2892320.647)));
 INSERT INTO ADJOINING_PARCELS (OBJECTID,SHAPE) VALUES (12,SDO_GEOMETRY(2003,82362,NULL,SDO_ELEM_INFO_ARRAY(1,1003,1),SDO_ORDINATE_ARRAY(359188.997,2892278.137,359159.047,2892286.136,359154.918,2892270.678,359153.37,2892264.881,359183.32,2892256.882,359185.647,2892265.593,359188.997,2892278.137)));
 INSERT INTO ADJOINING_PARCELS (OBJECTID,SHAPE) VALUES (14,SDO_GEOMETRY(2003,82362,NULL,SDO_ELEM_INFO_ARRAY(1,1003,1),SDO_ORDINATE_ARRAY(359183.32,2892256.882,359153.37,2892264.881,359147.693,2892243.626,359158.297,2892240.794,359167.326,2892238.382,359177.643,2892235.627,359180.063,2892244.687,359183.32,2892256.882)));
 INSERT INTO ADJOINING_PARCELS (OBJECTID,SHAPE) VALUES (15,SDO_GEOMETRY(2003,82362,NULL,SDO_ELEM_INFO_ARRAY(1,1003,1),SDO_ORDINATE_ARRAY(359194.674,2892299.392,359164.723,2892307.391,359162.659,2892299.662,359159.047,2892286.136,359188.997,2892278.137,359193.687,2892295.698,359194.674,2892299.392)));
 commit;
 DELETE FROM user_sdo_geom_metadata WHERE TABLE_NAME = 'ADJOINING_PARCELS';
 commit;
 INSERT INTO USER_SDO_GEOM_METADATA (TABLE_NAME,COLUMN_NAME,DIMINFO,SRID)
      VALUES ('ADJOINING_PARCELS','SHAPE',MDSYS.SDO_DIM_ARRAY(MDSYS.SDO_DIM_ELEMENT('X',358925.345,359200.351,1),MDSYS.SDO_DIM_ELEMENT('Y',2892235.627,2892414.914,1)),82362);
 commit;
 CREATE INDEX ADJOINING_PARCELS_SHAPE
           ON ADJOINING_PARCELS(SHAPE)
        IndexType IS mdsys.spatial_index
        Parameters('sdo_indx_dims=2,layer_gtype=POLYGON');
 -- Test
 SELECT p1.objectid,
        p2.objectid,
        CardinalDirection(
           sdo_geom.sdo_centroid(p1.shape,0.05),
           sdo_geom.sdo_centroid(p2.shape,0.05),
           0.05,
           1 /* Projected */) AS cardnlDirn,
        sdo_geom.sdo_centroid(p2.shape,0.05) target
   FROM ADJOINING_PARCELS p1,
        ADJOINING_PARCELS p2
  WHERE p1.objectid = 12
    AND sdo_anyinteract(p2.shape,p1.shape) = 'TRUE'
    AND p2.objectid <> p1.objectid
    AND sdo_geom.sdo_intersection(p1.shape,p2.shape,0.005).sdo_gtype <> 2001;
 -- Results
 --
 OBJECTID OBJECTID CARDNLDIRN
 -------- -------- ----------
       12        1 WSW
       12        8 NW
       12       14 SSW
       12       15 NNE

This results in the following visually.

I hope this is useful to someone, because my interlocutor thought so:

Thank you very much the code is working as I accepted and much more . God bless you that you are publishing your knowledge to others , you are the man!