Go to content Go to navigation and search

Home

Current Oracle Spatial Blogs

Search

RSS / Atom

Email me

textpattern

New To_3D Function

· Aug 12, 10:51 pm by Simon Greener

I often have need to convert 2D geometry objects from 2D to 3D and have wrote a function that was part of my free PL/SQL packages that did this. That function depended on another package.

Just this week I had need of a standalone version, so I went back to my old function and removed all dependencies on any other package.

This function, called TO_3D, actually does a number of things:

Here is the function.

   /*
   --  function To_3D
   --  precis   Converts a 2D or 4D geometry to a 3D geometry
   --  version  1.0
   --  usage    v_3D_geom := geom.To_3D(MDSYS.SDO_Geometry(2001,....),50)
   --  history  Simon Greener,   May 2007 Original coding
   --  history  Simon Greener,   Aug 2009 Added support for interpolating Z values 
   */
  Function To_3D( p_geom      IN MDSYS.SDO_Geometry,
                  p_start_z   IN NUMBER,
                  p_end_z     IN NUMBER := NULL,
                  p_tolerance IN NUMBER := 0.05)
     Return MDSYS.SDO_Geometry Deterministic
  Is
    v_sign              PLS_INTEGER := SIGN(p_end_z - p_start_z);
    v_isMeasured        BOOLEAN;
    v_gtype             INTEGER;     -- geometry type (single digit)
    v_dim               INTEGER;
    v_npoints           INTEGER;
    v_i                 PLS_INTEGER;
    v_j                 PLS_INTEGER;
    v_offset            PLS_INTEGER;
    v_length            NUMBER := 0;
    v_cumulative_length NUMBER := 0;
    v_round_factor      NUMBER := case when p_tolerance is null
                                       then null
                                       else round(log(10,(1/p_tolerance)/2))
                                   end;
    v_3D_geom           MDSYS.SDO_Geometry;

  Function isCompoundElement(p_elem_type in number)
    return boolean
  Is
  Begin
    Return ( p_elem_type in (4,5,1005,2005) );
  End isCompoundElement;

  Function isMeasured( p_gtype in number )
    return boolean
  is
  Begin
    Return CASE WHEN MOD(trunc(p_gtype/100),10) = 0
                THEN False
                ELSE True
             END;
  End isMeasured;

  Function GetDimensions( p_gtype in number )
    return integer
  Is
  Begin
    return TRUNC(p_gtype/1000,0);
  End GetDimensions;

  Begin
    -- If the input geometry is null, just return null
    IF ( p_geom IS NULL ) THEN
      RETURN (NULL);
    END IF;

    -- Get the number of dimensions and the gtype
    v_dim   := GetDimensions(p_geom.sdo_gtype);
    v_gtype := MOD(p_geom.sdo_gtype,10); -- Short gtype
    v_isMeasured := isMeasured(p_geom.sdo_gtype);
   
    IF ( v_dim = 3 ) THEN
      -- Nothing to do, p_geom is already 3D
      Return (p_geom);
    END IF;

    IF ( p_start_z is not null And p_end_z is not null ) THEN
       v_length := mdsys.sdo_geom.SDO_Length( p_geom, p_tolerance );
    END IF;

    -- Compute number of points
    v_npoints := mdsys.sdo_util.GetNumVertices(p_geom);

    -- Construct output object ...
    v_3D_geom           := MDSYS.SDO_GEOMETRY (3000 + v_gtype,
                             p_geom.sdo_srid,
                             p_geom.sdo_point,
                             p_geom.sdo_elem_info,
                             MDSYS.sdo_ordinate_array ()
                           );

    -- Does geometry have a valid sdo_point?
    IF ( v_3D_geom.sdo_point is not null ) Then
      -- It's a point, there's not much to it...
      v_3D_geom.sdo_point.Z := p_start_z;
    END IF;

    -- If Single Point, all done, else...    
    IF ( v_gtype != 1 ) THEN
      -- It's not a single point ...
      -- Process the geometry's ordinate array

      -- Create space in ordinate array for 3D coordinates
      v_3D_geom.sdo_ordinates.EXTEND ( v_npoints * 3 );

      -- Copy the ordinates array
      v_i := p_geom.sdo_ordinates.FIRST;      -- index into input ordinate array
      v_j := 1;                               -- index into output ordinate array
      v_cumulative_length := 0;
      FOR i IN 1 .. v_npoints LOOP
        v_3D_geom.sdo_ordinates (v_j)     := p_geom.sdo_ordinates (v_i);      -- copy X
        v_3D_geom.sdo_ordinates (v_j + 1) := p_geom.sdo_ordinates (v_i + 1);  -- copy Y
        -- compute and assign Z only if doesn't have existing value
        If ( ( v_dim = 4 ) Or ( v_dim = 3 And Not v_isMeasured ) ) Then
           v_3D_geom.sdo_ordinates (v_j + 2) := p_geom.sdo_ordinates (v_i + 2);  -- copy Z
        Else
           -- Compute new Z
          If ( i = 1 ) Then
            v_3D_geom.sdo_ordinates (v_j + 2) := p_start_z; 
          Else
            If ( v_length <> 0 ) Then
              v_cumulative_length := v_cumulative_length + 
                 round(sdo_geom.sdo_distance(mdsys.sdo_geometry(2001,
                                                                p_geom.sdo_srid,
                                                                sdo_point_type(p_geom.sdo_ordinates(v_i - v_dim),
                                                                               p_geom.sdo_ordinates(v_i - v_dim + 1),
                                                                               null),
                                                                null,null),
                                             mdsys.sdo_geometry(2001,
                                                                p_geom.sdo_srid,
                                                                sdo_point_type(p_geom.sdo_ordinates (v_i),
                                                                               p_geom.sdo_ordinates (v_i + 1),null),
                                                                null,null),
                                             p_tolerance),
                       v_round_factor);
            End If;
            If ( i = v_npoints ) Then
              v_3D_geom.sdo_ordinates (v_j + 2) := case when v_length = 0 then p_start_z else p_end_z end; 
            Else
              v_3D_geom.sdo_ordinates (v_j + 2) := case when p_end_z is null then p_start_z
                                                        when v_length != 0   then p_start_z + v_sign * round( ( ( p_end_z - p_start_z ) / v_length ) * v_cumulative_length,v_round_factor)
                                                        else p_start_z
                                                    end; 
            End If;
          End If;
        End If;
        v_i := v_i + v_dim;
        v_j := v_j + 3;
      END LOOP;

      -- Process the element info array
      -- by adjust the offsets
      v_i := v_3D_geom.sdo_elem_info.FIRST;
      WHILE v_i < v_3D_geom.sdo_elem_info.LAST LOOP
            If Not ( isCompoundElement(v_3D_geom.sdo_elem_info (v_i + 1) ) ) Then
              -- Adjust Elem Info offsets
              v_offset := v_3D_geom.sdo_elem_info (v_i);
              v_3D_geom.sdo_elem_info(v_i) := ( v_offset - 1 ) / v_dim * 3 + 1;
            End If;
            v_i := v_i + 2;
      END LOOP;

    END IF;
    RETURN v_3D_geom;
  END To_3D;

Now let’s test it by building a 3D geometry from a 2D one with full Z interpolation.

select geom.to_3d( a.original_geom, -1, -200, 0.05) as threed_geom
  from ( select mdsys.sdo_geometry(2002,null,null,sdo_elem_info_array(1,2,1),sdo_ordinate_array(0,0,5,5,10,10)) as original_geom
           from dual 
       ) a;

THREED_GEOM
-----------
MDSYS.SDO_GEOMETRY(3002,null,null,MDSYS.SDO_ELEM_INFO_ARRAY(1,2,1),MDSYS.SDO_ORDINATE_ARRAY(0,0,-1,5,5,98.9,10,10,-200))

1 rows selected

Now, take the above and push it into a function to add measure information to the new 3D geometry.

select MDSYS.SDO_LRS.CONVERT_TO_LRS_GEOM( threed_geom, 1, 10) as lrs_geom, threed_geom
  from (select geom.to_3d( a.original_geom, -1, -200, 0.05) as threed_geom
          from ( select mdsys.sdo_geometry(2002,null,null,sdo_elem_info_array(1,2,1),sdo_ordinate_array(0,0,5,5,10,10)) as original_geom
                   from dual 
               ) a
       ) b;

LRS_GEOM                                                                                                                          threed_geom
--------------------------------------------------------------------------------------------------------------------------------- ------------------------------------------------------------------------------------------------------------------------
MDSYS.SDO_GEOMETRY(4402,null,null,MDSYS.SDO_ELEM_INFO_ARRAY(1,2,1),MDSYS.SDO_ORDINATE_ARRAY(0,0,-1,1,5,5,98.9,5.5,10,10,-200,10)) MDSYS.SDO_GEOMETRY(3002,null,null,MDSYS.SDO_ELEM_INFO_ARRAY(1,2,1),MDSYS.SDO_ORDINATE_ARRAY(0,0,-1,5,5,98.9,10,10,-200))

1 rows selected

Now let’s strip the measure information off the geometry to return our 3D geometry of above.

select geom.to_3d(c.lrs_geom,null,null,0.05) as geom_3d, c.threed_geom
          from (select MDSYS.SDO_LRS.CONVERT_TO_LRS_GEOM( threed_geom, 1, 10) as lrs_geom, threed_geom
                  from (select geom.to_3d( a.original_geom, -1, -200, 0.05) as threed_geom
                          from ( select mdsys.sdo_geometry(2002,null,null,sdo_elem_info_array(1,2,1),sdo_ordinate_array(0,0,5,5,10,10)) as original_geom
                                   from dual 
                               ) a
                        ) b
               ) c;

GEOM_3D                                                                                                                  threed_geom
------------------------------------------------------------------------------------------------------------------------ ------------------------------------------------------------------------------------------------------------------------
MDSYS.SDO_GEOMETRY(3002,null,null,MDSYS.SDO_ELEM_INFO_ARRAY(1,2,1),MDSYS.SDO_ORDINATE_ARRAY(0,0,-1,5,5,98.9,10,10,-200)) MDSYS.SDO_GEOMETRY(3002,null,null,MDSYS.SDO_ELEM_INFO_ARRAY(1,2,1),MDSYS.SDO_ORDINATE_ARRAY(0,0,-1,5,5,98.9,10,10,-200))

1 rows selected

Finally, let’s check the validity of the final geometry which should be the same as the original.

select sdo_geom.validate_geometry(d.geom_3d,0.05)
  from (select geom.to_3d(c.lrs_geom,null,null,0.05) as geom_3d, c.threed_geom
          from (select MDSYS.SDO_LRS.CONVERT_TO_LRS_GEOM( threed_geom, 1, 10) as lrs_geom, threed_geom
                  from (select geom.to_3d( a.original_geom, -1, -200, 0.05) as threed_geom
                          from ( select mdsys.sdo_geometry(2002,null,null,sdo_elem_info_array(1,2,1),sdo_ordinate_array(0,0,5,5,10,10)) as original_geom
                                   from dual 
                               ) a
                        ) b
               ) c
        ) d;

SDO_GEOM.VALIDATE_GEOMETRY(D.GEOM_3D,0.05)
------------------------------------------
TRUE

1 rows selected

Finally, let’s reduce the 3D geometry back to a 2D geometry using a TO_2D function to see if we get the same geometry.

select sdo_geom.relate(geom.to_2D(d.geom_3d),'DETERMINE',mdsys.sdo_geometry(2002,null,null,sdo_elem_info_array(1,2,1),sdo_ordinate_array(0,0,5,5,10,10)),0.05) as relate
  from (select geom.to_3d(c.lrs_geom,null,null,0.05) as geom_3d, c.threed_geom
          from (select MDSYS.SDO_LRS.CONVERT_TO_LRS_GEOM( threed_geom, 1, 10) as lrs_geom, threed_geom
                  from (select geom.to_3d( a.original_geom, -1, -200, 0.05) as threed_geom
                          from ( select mdsys.sdo_geometry(2002,null,null,sdo_elem_info_array(1,2,1),sdo_ordinate_array(0,0,5,5,10,10)) as original_geom
                                   from dual 
                               ) a
                        ) b
               ) c
        ) d;

RELATE
---------
EQUAL

1 rows selected

I hope this function is of use to someone as it certainly is useful to me!

post this at del.icio.uspost this at Diggpost this at Technoratipost this at Newsvinepost this at Ma.gnoliapost this at Furlpost this at Blinklistpost this at Spurlpost this at Wistspost this at Simpypost this at Redditpost this at Farkpost this at Blogmarkspost this at Yahoo! my webpost this at Mr. Wongpost this at Windows Livepost this at Google Bookmarkspost this to Twitter
  Textile Help