Sample Oracle Spatial Queries

From NGDCWiki

Jump to: navigation, search

Contents

point-in-polygon, user-specified area, geodetic data

  • in 10g, optimized geodetic rectangles are supported and are the preferred way to do an envelope query
select count(*) from JCC.CITIES point where
   SDO_FILTER(
      point.SHAPE,
      sdo_geometry(
         2003,
         8307,
         null,
         mdsys.sdo_elem_info_array(1,1003,3),
         mdsys.sdo_ordinate_array(-180.0,-90.0, 180.0,90.0)
      )
   ) = 'TRUE'

access the components of a SDO_GEOMETRY instance directly

  • polygon, sdo_point field is normally null
select 
   a.objectid, 
   a.shape.sdo_point, 
   a.shape.sdo_elem_info, 
   a.shape.sdo_ordinates, 
   a.shape.sdo_gtype 
from polygon_table a
  • get a point's X,Y values
select p.SHAPE.SDO_POINT.X, p.SHAPE.SDO_POINT.Y from randompoints p

useful methods on SDO_GEOMETRY

returns 0 if invalid, 1 if valid

select s.shape.st_isvalid() from poly_tst s

returns the Well Known Text representation of the geometry

select s.shape.get_wkt() from poly_tst s

return the GType of the geometry

select s.shape.get_gtype() from poly_tst s

GTypes:

1Point
2Line
3Polygon
5Multipoint
6Multiline
7Multipolygon

Create a multipart polygon

use SDO_UTIL.APPEND to create a multipart poly from two simple polygons. Note that optimized rectangles are expanded in the output geometry

insert into poly_tst values (
   3,
   sdo_util.append(
      sdo_geometry(
         2003,
         8307,
         null,
         sdo_elem_info_array(1,1003,3),
        sdo_ordinate_array(10,10, 20,20)
      ),
      sdo_geometry(
         2003,8307,null,
         sdo_elem_info_array(1, 1003, 3),
        sdo_ordinate_array(30,30, 40,40)
      )
   )
)

find features that cross the antemeridian or poles

select
   s.objectid, s.shape
from
   poly_tst s
where
   sdo_overlaps(
      s.shape,
      sdo_geometry(
         2003,
         8307,
         null,
         sdo_elem_info_array(1,1003,3),
         sdo_ordinate_array(-180,-90, 180,90)
      )
   ) = 'TRUE'

Use Oracle to re-project data

convert point in geographic (NAD83) coords to UTM zone14:

select sdo_cs.transform(SDO_GEOMETRY(2001,8265,
     MDSYS.SDO_POINT_TYPE(-102, 24, NULL),
     NULL,
     NULL
  ),82232) from dual

Validate a Layer

validate table RANDOMPOINTS and store results of validation in RANDOM_POINT_RESULTS. The results tablename should be uppercase

  • create a table to hold validation results
CREATE TABLE random_point_results (sdo_rowid ROWID, result VARCHAR2(2000));
  • validate the geometries
EXECUTE SDO_GEOM.VALIDATE_LAYER_WITH_CONTEXT('RANDOMPOINTS','SHAPE','RANDOM_POINT_RESULTS',10000);
  • join the geometry table with the validation results
select 
  A.OBJECTID, B.RESULT 
from 
  POLY0 A, RESULTS B 
where 
  A.ROWID = B.SDO_ROWID;

List of spatial indexes and their types

select INDEX_NAME,TABLE_NAME,SDO_INDEX_TYPE,SDO_INDEX_STATUS from user_sdo_index_info;

note that a VALID SDO_INDEX_STATUS does not always mean that everything is OK!

Find points within fixed distance of feature

select cities within 250km of Denver

select p.city_name from cities p where
   sdo_within_distance(p.shape,
                       (select shape from cities where city_name = 'Denver'),
                       'distance=250 unit=km'
   ) = 'TRUE';

Find the 3 cities closest to Denver

note the result count is set to 4 because the origin is included

select p.city_name from cities p where 
   sdo_nn(p.shape, 
          (select shape from cities where city_name = 'Denver'),
          'sdo_num_res=4'
   ) = 'TRUE';

Split geometries at the antemeridian

provided by Daniel Geringer

REM
REM  If east_west_or_multi = EAST, function returns the portion left of the 180 meridian
REM     which would display on the EAST side of a flat representation of the world.
REM
REM  If east_west_or_multi = WEST, function returns the portion right of the 180 meridian.
REM     which would display on the WEST side of a flat representation of the world.
REM
REM  If east_west_or_multi = MULTI, function returns a multi polygon with with both portions.
REM
REM  Note: The multipolygon can be used for ANYINTERACT queries, and for DISPLAY, but 
REM        INSIDE queries will not return objects that interact with the 180 meridian
REM        due to the seam in the multipolygon geometry. Also, the multipolgon geometry
REM        will not pass validation because it will violate the OGC rule
REM        that states outer rings in a multipolygon can only intersect at a point.
REM
CREATE OR REPLACE FUNCTION split_geom (in_geom            SDO_GEOMETRY,
                                       east_west_or_multi VARCHAR2) 
  RETURN SDO_GEOMETRY DETERMINISTIC AS

  -- West and east boxes each cover 1/4 of the world
  --   Rotation must be counterclockwise when defining these polygons.
  west_box  SDO_GEOMETRY := sdo_geometry (2003,8265,null,
                              sdo_elem_info_array (1,1003,1),
                              sdo_ordinate_array (-180,90,-180,0,-180,-90,-90,0,-90,90));
  east_box  SDO_GEOMETRY := sdo_geometry (2003,8265,null,
                              sdo_elem_info_array (1,1003,1),
                              sdo_ordinate_array (0, 90, 90, 0, 0, -90, 180,0, 180,90));
BEGIN
  IF upper (east_west_or_multi) = 'EAST'
  THEN
    return sdo_geom.sdo_intersection (in_geom, east_box, .05);
  ELSIF upper (east_west_or_multi) = 'WEST'
  THEN
    return sdo_geom.sdo_intersection (in_geom, west_box, .05);
  ELSE
    return sdo_util.append (sdo_geom.sdo_intersection (in_geom, east_box, .05),
                            sdo_geom.sdo_intersection (in_geom, west_box, .05));
  END IF;
END;
/
REM
REM Example on how to call split_geom
REM
set pages 1000
SELECT split_geom (shape, 'EAST') from DR_SURVEY_CROSS180 where objectid = 2097;
SELECT split_geom (shape, 'WEST') from DR_SURVEY_CROSS180 where objectid = 2097;
SELECT split_geom (shape, 'MULTI') from DR_SURVEY_CROSS180 where objectid = 2097;
REM
REM Example on how to find all the geometries that cross the 180 meridian,
REM   and then call split_geom.
REM
DECLARE
  west_geom  SDO_GEOMETRY;
  east_geom  SDO_GEOMETRY;
  multi_geom SDO_GEOMETRY;
BEGIN
  FOR r IN (select objectid, shape
            from DR_SURVEY_CROSS180
            where sdo_anyinteract (shape, 
                                   sdo_geometry (2002,8265,null,sdo_elem_info_array (1,2,1),
                                                 sdo_ordinate_array (-180,90, 
                                                                     -180,0, 
                                                                     -180,-90))) = 'TRUE')
  LOOP
    east_geom  := split_geom (r.shape, 'EAST');
    west_geom  := split_geom (r.shape, 'WEST');
    multi_geom := split_geom (r.shape, 'MULTI');
  END LOOP;
END;


slight modification to prevent parts from sharing a coordinate:


create or replace
FUNCTION split_geom (in_geom SDO_GEOMETRY,
                                       east_west_or_multi VARCHAR2) 
  RETURN SDO_GEOMETRY DETERMINISTIC AS

  -- West and east boxes each cover 1/2 of the world
  --   Rotation must be counterclockwise when defining these polygons.
  west_box  SDO_GEOMETRY := sdo_geometry (2003,8307,null,
                              sdo_elem_info_array (1,1003,1),
                              sdo_ordinate_array (-179.999,89.999,
                                                  -179.999,-89.999,
                                                  -0.001,-89.999,
                                                  -0.001, 89.999,
                                                  -179.999,89.999));
  east_box  SDO_GEOMETRY := sdo_geometry (2003,8307,null,
                              sdo_elem_info_array (1,1003,1),
                              sdo_ordinate_array (0.001, 89.999,
                                                  0.001,-89.999,
                                                  179.999,-89.999,
                                                  179.999,89.999,
                                                  0.001, 89.999));
BEGIN
  IF upper (east_west_or_multi) = 'EAST'
  THEN
    return sdo_geom.sdo_intersection (in_geom, east_box, .05);
  ELSIF upper (east_west_or_multi) = 'WEST'
  THEN
    return sdo_geom.sdo_intersection (in_geom, west_box, .05);
  ELSE
    return sdo_util.append (sdo_geom.sdo_intersection (in_geom, east_box, .05),
                            sdo_geom.sdo_intersection (in_geom, west_box, .05));
  END IF;
END;

then call it with:

update FEATURE set SHAPE = 
    split_geom (shape, 'MULTI')
 
WHERE 
  sdo_anyinteract (
    shape, 
    sdo_geometry (
      2002,
      8307,
      NULL,
      sdo_elem_info_array (1,2,1), 
      sdo_ordinate_array (-180,90, -180,0, -180,-90)
      )
    ) = 'TRUE';

Combine Multiple Geometries into One

Create a single "unioned" geometry from the inputs. If the inputs are discontinuous, the output geometry will be a multipart polygon

 select sdo_aggr_union(sdoaggrtype(a.shape,0.05)) from NOS_BAG_TSQA a where SURVEY_ID = 'H11413';

Create an envelope which contains all the input geometries

 select sdo_aggr_mbr(a.shape) from NOS_BAG_TSQA a where  SURVEY_ID = 'H11413';

Another example which produces the multipart polygon output

 drop table POLY0;
 
 create table POLY0 (OBJECTID number(38), SHAPE MDSYS.SDO_GEOMETRY);
 
 insert into POLY0 values (
 1,
 MDSYS.SDO_GEOMETRY(
   2003, -- 2D polygon
   8307, -- SRID
   NULL,
   MDSYS.SDO_ELEM_INFO_ARRAY(1,1003,1), -- one polygon (exterior)
   MDSYS.SDO_ORDINATE_ARRAY(
     -10,-10,
     10,-10,
     10,10,
     -10,10,
     -10,-10
   )));
 
 insert into POLY0 values (
 2,
 MDSYS.SDO_GEOMETRY(
   2003, -- 2D polygon
   8307, -- SRID
   NULL,
   MDSYS.SDO_ELEM_INFO_ARRAY(1,1003,1), -- one polygon (exterior)
   MDSYS.SDO_ORDINATE_ARRAY(
     20,20,
     30,20,
     30,30,
     20,30,
     20,20
   )));
 
 insert into user_sdo_geom_metadata values (
  'POLY0',
  'SHAPE',
  MDSYS.SDO_DIM_ARRAY(
     MDSYS.SDO_DIM_ELEMENT('Longitude',-180,180,0.5),
     MDSYS.SDO_DIM_ELEMENT('Latitude',-90,90,0.5)
  ),8307);
 
 create index poly0_spidx on poly0(SHAPE) indextype is MDSYS.SPATIAL_INDEX;
 
 select objectid,sdo_geom.validate_geometry_with_context(p.shape,0.5) from poly0 p;
 
 select sdo_aggr_union(sdoaggrtype(a.shape,0.05)) from POLY0 a;
 
 select sdo_aggr_mbr(a.shape) from POLY0 a;    -- MDSYS.SDO_GEOMETRY(2003,8307,'null',
                                               -- MDSYS.SDO_ELEM_INFO_ARRAY(1,1003,3),
                                               -- MDSYS.SDO_ORDINATE_ARRAY(....

SDO_AGGR_UNION should not be used on lines, use SDO_AGGR_CONCAT_LINES, e.g.

 select
    sdo_util.simplify( sdo_aggr_concat_lines(a_shape), 250, 0.05)
  from (
    select a.shape a_shape from mb.mbinfo_file_tsql a
    where ngdc_id = 'NEW1492');

Note the above combines a simplify operation with the concatenation to produce a multipart line with reduced number of vertices

use SDO_JOIN to aggregate points into "spatial bins"

provided by Daniel Geringer

 DROP TABLE results;
 CREATE TABLE results  NOLOGGING AS
    SELECT b.objectid rp_objectid, c.objectid one_deg_object_id
    FROM TABLE (sdo_join('RANDOM_POINTS', 'SHAPE', 'ONE_DEG', 'SHAPE',
                         'mask=anyinteract')) a,
         random_points b,
         one_deg c
    WHERE a.rowid1 = b.rowid
      AND a.rowid2 = c.rowid;

 -- Query to find random points that interact with more than one cell
 SELECT rp_objectid, count(*)
    FROM results
    GROUP BY rp_objectid
    HAVING count(*) > 1;
Personal tools