Example 6-1 Setting Up the GCP Georeferencing Model

For example, if you have a Landsat image in a plain area and want to georeference it, you might choose the Quadratic Polynomial geometric model. For that purpose, assuming you have collected 9 GCPs (at least 6 GCPs in this case) and 3 check points, you can set up the GCPs and store them in the GeoRaster's metadata using the code in Example 6-1.

DECLARE
   gr1         sdo_georaster;
   georefModel SDO_GEOR_GCPGEOREFTYPE;
   GCPs        SDO_GEOR_GCP_COLLECTION;
BEGIN
   SELECT georaster INTO gr1 from georaster_table WHERE georid=1 FOR UPDATE;
   GCPs := SDO_GEOR_GCP_COLLECTION( 
                 SDO_GEOR_GCP('1', '', 1, 
                      2, sdo_number_array(25, 73),
                      2, sdo_number_array(237036.9, 897987.2),
                      NULL, NULL),
                 SDO_GEOR_GCP('2', '', 1, 
                      2, sdo_number_array(100, 459),
                      2, sdo_number_array(237229.6, 897949.7),
                      NULL, NULL),
                 SDO_GEOR_GCP('3', '', 1, 
                      2, sdo_number_array(362, 77),
                      2, sdo_number_array(237038.9, 897818.8),
                      NULL, NULL),
                 SDO_GEOR_GCP('4', '', 1, 
                       2, sdo_number_array(478, 402),
                       2, sdo_number_array(237201.06, 897760.56),
                       NULL, NULL),
                 SDO_GEOR_GCP('5', '', 1, 
                       2, sdo_number_array(167, 64),
                       2, sdo_number_array(237032.02, 897916.26),
                       NULL, NULL),
                 SDO_GEOR_GCP('6', '', 1, 
                       2, sdo_number_array(101, 257),
                       2, sdo_number_array(237128.9, 897949.3),
                       NULL, NULL),
                 SDO_GEOR_GCP('7', '', 1, 
                       2, sdo_number_array(235, 501),
                       2, sdo_number_array(237250.9, 897882.2),
                       NULL, NULL),
                 SDO_GEOR_GCP('8', '', 1, 
                       2, sdo_number_array(423, 214),
                       2, sdo_number_array(237107.3, 897788.0),
                       NULL, NULL),
                 SDO_GEOR_GCP('9', '', 1, 
                       2, sdo_number_array(127, 178),
                       2, sdo_number_array(237089.0, 897936.5),
                       NULL, NULL),
                 SDO_GEOR_GCP('10', '', 2, 
                       2, sdo_number_array(131, 425),
                       2, sdo_number_array(237212.8, 897934.2),
                       NULL, NULL),
                 SDO_GEOR_GCP('11', '', 2, 
                       2, sdo_number_array(299, 111),
                       2, sdo_number_array(237055.7, 897850.4),
                       NULL, NULL),
                 SDO_GEOR_GCP('12', '', 2, 
                       2, sdo_number_array(329, 253),
                       2, sdo_number_array(237126.9,897835.4), 
                       NULL, NULL) );
   georefModel := SDO_GEOR_GCPGEOREFTYPE('QuadraticPolynomial', GCPs.count, GCPs, NULL);
   -- Set and store the GCP georeference model into the GeoRaster object's metadata
   sdo_geor.setGCPGeorefModel(gr1, georefModel);
   UPDATE georaster_table SET georaster=gr1 WHERE georid=1;
   COMMIT;
END;
/

Example 6-2 Generating the Functional Fitting Model Using GCPs

After using the code in Example 6-1, you can generate the functional fitting model coefficients by using the code in Example 6-2.

DECLARE
   gr1  sdo_georaster;
   rms  sdo_number_array;
BEGIN
   SELECT georaster INTO gr1 from georaster_table WHERE georid=1 FOR UPDATE;
   -- georeference the image using the GCPs stored in the image's metadata
   rms := sdo_geor.georeference(gr1, null, 26986, 0, 'TRUE'); 
   UPDATE georaster_table SET georaster=gr1 WHERE georid=1;
   COMMIT;
END;
/

The steps in Example 6-1 and Example 6-2 can be combined without the need to pre-set the GCPs into the GeoRaster object's metadata (see the example for SDO_GEOR.georeference in SDO_GEOR Package Reference). The returned value array of SDO_GEOR.georeference in Example 6-2 contains RMS values and residuals for each GCP. Using these, you can examine the solution accuracy and identify erratic GCPs. If the accuracy is not satisfactory, recheck all GCPs to make sure they are accurate and add more GCPs as necessary, and then run the script or scripts again.

The GCP support in GeoRaster enables you to spatially reference any non-geospatial images and rasters also.

After geospatial images are georeferenced, you can process those images, such as applying rectification, reprojection, and mosaicking, and spatially querying and subsetting the rasters using geometry polygons in different coordinate systems.

Example 6-3 Image Reprojection

Example 6-3 reprojects a raster image that had been loaded into a GeoRaster object with SRID 4326, but needs to be reprojected to have the same SRID as other images previously stored with SRID 23619.

DECLARE
  gr1 sdo_georaster;
  gr2 sdo_georaster;
BEGIN
  select raster into gr1 from georaster_load_table where georid = 10;
  delete from georaster_table where georid = 54;
  insert into georaster_table 
         values(54,'reprojected', sdo_geor.init()) 
         returning georaster into gr2;
  sdo_geor.reproject(inGeoRaster   => gr1,
                     pyramidLevel  => 0,
                     cropArea      => null,
                     layerNumbers  => null,
                     resampleParam => 'resampling=BILINEAR',
                     storageParam  => null,
                     outSRID       => 32619,
                     outGeoraster  => gr2);
  update georaster_table set georaster = gr2 where georid = 54;
  commit;
END;

The same operation can be accomplished by the SDO_GEOR.rectify procedure, producing similar results. The SDO_GEOR.rectify procedure offers more capabilities and flexibility than SDO_GEOR.reproject; for example, the input and output SRID can be the same and users can specify the precise resolution of the output (see Image Rectification).

If a GeoRaster object does not have an associated SRS, the process for georeferencing and rectifying it is explained in Georeferencing GeoRaster Objects and Image Rectification.

Parallel reprojection is supported with the SDO_GEOR_AGGR.mosaicSubset procedure.

Example 6-4 Image Rectification

Example 6-3 rectifies an aerial image that had been loaded into GeoRaster and later georeferenced with GCPs (see Advanced Georeferencing). The image is rectified so that the output GeoRaster object has the same SRS and resolution of an existing GeoRaster object. The image is to be restricted to the area of existing GeoRaster object, and the pixels should be perfectly aligned with the existing GeoRaster object.

DECLARE
  gr_src sdo_georaster;
  gr_ref sdo_georaster;
  gr_out sdo_georaster;
BEGIN
  select raster into gr_src from georaster_load_table where georid = 15;
  select raster into gr_ref from georaster_table where georid = 1;
  delete from georaster_table where georid = 2;
  insert into georaster_table 
         values(2, 'rectified', sdo_geor.init()) 
         returning georaster into gr_out;
  sdo_geor.rectify(inGeoRaster      => gr_src,
                   pyramidLevel     => null,
                   elevationParam   => null,
                   dem              => null,
                   outSRID          => sdo_geor.getModelSRID(gr_ref),
                   outModelCoordLoc => null,
                   cropArea         => sdo_geor.generateSpatialExtent(gr_ref),
                   polygonClip      => null,
                   layerNumbers     => null,
                   outResolutions   => sdo_geor.getSpatialResolutions(gr_ref),
                   resolutionUnit   => 'unit=meters',
                   referencePoint   => sdo_geor.getModelCoordinate(gr_ref,
                                       0, sdo_number_array(-0.5,-0.5)),
                   resampleParam    => null,
                   storageParam     => null,
                   outGeoraster     => gr_out);
  update georaster_table set georaster = gr_out where georid = 2;
  commit;
END;

Rectification output can be significantly improved if information about elevation is passed to the SDO_GEOR.rectify procedure. (See Image Orthorectification for more information about elevation.)

Parallel rectification is supported with the SDO_GEOR_AGGR.mosaicSubset procedure.

Example 6-8 Image Warping

The following example uses the SDO_GEOR_SRS information from one GeoRaster image (gr1) as a reference to transform an existing GeoRaster object (gr2) into a new (warped) GeoRaster object (gr3). Thus, the third GeoRaster object is a “copy” (actually, a transformation) of the second GeoRaster object, but reflects the same georeferencing as the first GeoRaster object.

DECLARE
  srs sdo_geor_srs;
  gr1 sdo_georaster;
  gr2 sdo_georaster;
  gr3 sdo_georaster;
BEGIN
  select georaster into gr1 from georaster_table where georid = 1;
  select georaster into gr2 from georaster_table where georid = 2;

  srs := sdo_geor.getSRS(gr1); -- get the SRS from image 1.

  insert into georaster_table values(3, 'Warped Object',
         sdo_geor.init('imagery_rdt')) 
         returning georaster  into gr3;

  sdo_geor.warp( inGeoRaster      => gr2,
                 pyramidLevel     => null,
                 outSRS           => srs, -- apply SRS to warp transformation
                 cropArea         => null,
                 dimensionSize    => null,
                 layerNumbers     => null,
                 elevationParam   => null,
                 resampleParam    => ‘resampling=AVERAGE4’,
                 storageParam     => ‘pyramid=true’,
                 outGeoRaster     => gr3,
                 bgValues         => sdo_number_array(0,0,0),
                 parallelParam    => ‘parellel=4’ );

  update georaster_table set georaster = gr3 where georid = 3;
  commit;
END;

Example 6-9 Image Scaling Using SDO_GEOR.scaleCopy

This example performs rescaling by using SDO_GEOR.scaleCopy and specifying the scaleFactor value as 2. The input image will have 2 times more rows and 2 times more columns than the original, and the values will be resampled by the average16 algorithm. Note that the image will be 4 times larger than the original.

DECLARE
  gr_src  sdo_georaster;
  gr_out  sdo_georaster;
BEGIN
  select georaster into gr_src from georaster_table where georid = 7;
  -- Rescale
  delete from georaster_table where georid = 9;
  insert into georaster_table values(9, 're-scaled by scaleCopy',
              sdo_geor.init('rdt_4',9)) 
              returning georaster into gr_out;
  sdo_geor.scaleCopy(inGeoRaster      => null,
                     scaleFactor      => 'scaleFactor=2',
                     resampleParam    => 'resampling=AVERAGE16',
                     storageParam     => null,
                     outGeoraster     => gr_out);
  update georaster_table set georaster = gr_out where georid = 9;
  commit;
END;
/

Example 6-10 Image Scaling Using SDO_GEOR.rectify

This example performs rescaling by using SDO_GEOR.rectify and specifying the outResolutions parameter. The input image is already rectified, and the output will have the same SRID as the input.

DECLARE
  gr_src     sdo_georaster;
  gr_out     sdo_georaster;
BEGIN
  select georaster into gr_src from georaster_table where georid = 7;
  -- Rescale
  delete from georaster_table where georid = 10;
  insert into georaster_table values(10, 're-scaled by rectify',
              sdo_geor.init('rdt_4',10)) 
              returning georaster into gr_out;
  sdo_geor.rectify(inGeoRaster      => null,
                   pyramidLevel     => null,
                   elevationParam   => null,
                   dem              => null,
                   outSRID          => null,
                   outModelCoordLoc => null,
                   cropArea         => null,
                   polygonClip      => null,
                   layerNumbers     => null,
                   outResolutions   => sdo_number_array(1.2,1.2),
                   resolutionUnit   => null,
                   referencePoint   => null,
                   resampleParam    => 'resampling=CUBIC',
                   storageParam     => null,
                   outGeoraster     => gr_out,
                   parallelParam    => 'parallel=4');
  update georaster_table set georaster = gr_out where georid = 10;
  commit;
END;
/

Example 6-11 Rescaling Using SDO_GEOR.affineTransform

This example performs rescaling by using the SDO_GEOR.affineTransform procedure and specifying he scales parameter as sdo_number_array(2, 2), indicating that the image will be enlarged 2 times on the rows dimension and 2 times on the columns dimension.

DECLARE
  gr1 sdo_georaster;
  gr2 sdo_georaster;  
BEGIN
  select georaster into gr1 from georaster_table where georid = 1;

  insert into georaster_table values(2, 'Rotated 90 left',
         sdo_geor.init('rdt0',2)) returning georaster into gr2;

  sdo_geor.affineTransform(inGeoRaster   => gr1,
                           translation   => null,
                           scales        => sdo_number_array(2,2),
                           rotatePt      => null,
                           rotateAngle   => null,
                           shear         => null,
                           reflection    => null,
                           storageParam  => null,
                           outGeoraster  => gr2,
                           parallelParam => 'parallel=4');

  update georaster_table set georaster = gr2 where georid = 2;
  commit;
END;

Example 6-12 Parallel Generation of Pyramids

For faster pyramid generation, you can specify the parallelParam parameter with the SDO_GEOR.generatePyramid procedure. In Example 6-12, the degree of parallelism is set to 4. (The actual performance improvement for pyramid generation depends on the number of CPUs available to Oracle Database.)

DECLARE
  gr sdo_georaster;
BEGIN
 
  SELECT georaster INTO gr 
    FROM georaster_table WHERE georid = 6 FOR UPDATE;
 
  -- Generate pyramids.
  sdo_geor.generatePyramid(gr, 'rLevel=5, resampling=NN', null, 'parallel=4');
 
  -- Update the original GeoRaster object.
  UPDATE georaster_table SET georaster = gr WHERE georid = 6;
 
  COMMIT;
END;
/

To enable parallel processing of the pyramid generation, SDO_GEOR.generatePyramid performs an implicit commit operation. If an error during the call, the GeoRaster object may in an invalid state. If this occurs, use SDO_GEOR.deletePyramid to remove the newly generated and upper pyramid levels of the GeoRaster object.

Example 6-13 Partial Updating of Pyramids

You can partially update pyramids by using the SDO_GEOR.updateRaster procedure. In Example 6-13, the target GeoRaster object at a specified area (targetArea is specified as area) is updated by another GeoRaster object. The updateUpperPyramids parameter is set to true, so the upper pyramids of the target GeoRaster object are only partially updated at the specified area. In other words, the upper pyramid levels are not regenerated in full, but only the cells in that targetArea are regenerated, and thus performance is improved.

DECLARE
  gr1 sdo_georaster;
  gr2 sdo_georaster;
  area sdo_number_array := sdo_number_array(-200,-50,201,162);
BEGIN
  SELECT georaster INTO gr2 FROM georaster_table WHERE georid=0 FOR UPDATE;
  SELECT georaster INTO gr1 FROM georaster_table WHERE georid=1;
  SDO_GEOR.updateRaster(gr2, 0, null, area, gr1, 0, null, 'true');
  UPDATE GEORASTER_TABLE SET georaster=gr2 WHERE georid=0;
  COMMIT;
END;
/

Other techniques to speed up and automate the pyramiding process include batch processing and concurrent processing can be used. To batch pyramid many images in a certain area, see the example in Querying and Searching GeoRaster Objects. To process many batches concurrently, you can start different database sessions

Example 6-14 Vegetation Index Computation

Example 6-14takes a Landsat 7 ETM+ image and computes the NDVI with parallelism. The result is stored as another raster of floating number data type. Note that in the GeoRaster algebra language, band numbering starts with 0, so the formula translates into the expression: ({3}-{2})/({3}+{2}).

DECLARE 
  geor1    SDO_GEORASTER; 
  geor2    SDO_GEORASTER; 
EBGIN
  -- Source ETM+ image
  select georaster into geor1 from georaster_table where georid = 2; 
  -- Store NDVI 
  select georaster into geor2 from georaster_table where georid = 3 for update; 
  sdo_geor_ra.rasterMathOp(geor1,
       SDO_STRING2_ARRAY('({3}-{2})/({3}+{2})'),
       'celldepth=32bit_real',geor2, null, null, ‘parallel=4’);
  update georaster_table set georaster = geor2 where georid = 3;
  commit;
end;
/

In addition to NDVI, there are many other vegetation indexes in the area of remote sensing. Many of these can be similarly computed using the GeoRaster raster algebra.

Example 6-15 Tasseled Cap Transformation

Example 6-15 takes the DN data of a Landsat 5 TM image as input, executes the TCT using the GeoRaster raster algebra with parallelism, and creates a new image holding the results.

DECLARE
    gr1 sdo_georaster;
    gr2 sdo_georaster;
    ret varchar2(32);
BEGIN
    select georaster into gr1 from georaster_table where georid = 2;
    select georaster into gr2 from georaster_table where georid = 4 for update;
    sdo_geor_ra.rasterMathOp(
      gr1,
      SDO_STRING2_ARRAY(
         '0.3561*{0}+0.3972*{1}+0.3904*{2}+0.6966*{3}+0.2286*{4}+0.1596*{6}',
         '(-0.3344)*{0}-0.3544*{1}-0.4556*{2}+0.6966*{3}-0.0242*{4}-0.2630*{6}',
         '0.2626*{0}+0.2141*{1}+0.0926*{2}+0.0656*{3}-0.7629*{4}-0.5388*{6}',
         '0.0805*{0}-0.0498*{1}+0.1950*{2}-0.1327*{3}+0.5752*{4}-0.7775*{6}',
         '(-0.7252)*{0}-0.0202*{1}+0.6683*{2}+0.0631*{3}-0.1494*{4}-0.0274*{6}',
         '0.4000*{0}-0.8172*{1}+0.3832*{2}+0.0602*{3}-0.1095*{4}+0.0985*{6}' ),
      'celldepth=32BIT_REAL',
      gr2, null, null, ‘parallel=4’); 
    update georaster_table set georaster = gr2 where georid = 4;
    commit;
END;
/

You can also use the same raster algebra language to add code in Example 6-15 to convert the 32-bit floating number image into an 8-bit integer image and to apply image stretching (described in Image Stretching) on the resulting TCT image to generate a new GeoRaster object for visualization and analysis.

In addition to using the optimized implementation of raster algebra algorithms and the embedded parallel processing, you can further take advantage of the Oracle grid computing infrastructure to quickly compute NDVI or apply TCT on thousands of images stored in the GeoRaster database.

Example 6-16 Band Merging

Example 6-16 includes two examples. The first example assumes there are eight GeoRaster objects, each of which contains only one band loaded from a single-band Landsat ETM+ image file in GeoTIFF format. The number of the band in each GeoRaster object is the same as the GEORID column value for the GeoRaster object. The example merges all bands into a single GeoRaster object to create a complete ETM+ scene.

DECLARE
    gr1 sdo_georaster;
BEGIN
    select georaster into gr1 from georaster_table where georid = 1 for update;
    for rec in (select georaster from georaster_table 
                     where georid >= 2 and georid <= 8 
                     order by georid)
    loop
       sdo_geor.mergelayers(gr1, rec.georaster);
    end loop;
    update georaster_table set georaster = gr1 where georid = 1;
    commit;
END;
/

The second example assumes there are eight GeoRaster objects, each of which contains three bands. The example picks up one band from each GeoRaster object and merges them into a single 8-band GeoRaster object in parallel.

DECLARE
  geor       SDO_GEORASTER;
  geo_array  SDO_GEORASTER_ARRAY;
BEGIN
  SELECT georaster INTO geor FROM georaster_table WHERE georid = 0 for update;
  geo_array:=SDO_GEORASTER_ARRAY();
  for rec in (select georaster from georaster_table 
                     where georid >= 1 and georid <= 8 
                     order by georid)
  loop
     geo_array.extend(1);
     geo_array(geo_array.last):=rec.georaster;
  end loop;
  sdo_geor_ra.rasterMathOp(geo_array,SDO_STRING2_ARRAY('{0,0}','{1,1}','{2,2}','{3,0}','{4,1}','{5,2}','{6,0}','{7,1}',),null,geor,'false',null,'parallel=4');
  UPDATE georaster_table SET georaster = geor WHERE georid = 0;
  COMMIT;
END;
/

Example 6-17 Appending One Image to Another Image

DECLARE
    gr1 sdo_georaster;
    gr1 sdo_georaster;
BEGIN
     select georaster into gr1 from georaster_table where georid = 1 for update;
     select georaster into gr2 from georaster_table where georid = 2;
     sdo_geor_aggr.append(gr1, gr2, 0, 'nodata=true');
     update georaster_table set georaster = gr1 where georid= 1;
     commit;
END;
/

Example 6-18 SDO_GEOR.mosaic (Table and Column Name)

Example 6-18 shows the SDO_GEOR.mosaic procedure.

DECLARE
  gr sdo_georaster;
BEGIN
  INSERT INTO georaster_table (georid, georaster) 
      VALUES (12, sdo_geor.init('rdt_1'))
      RETURNING georaster INTO gr;
  sdo_geor.mosaic('grtab', 'grobj', gr, 'blocking=optimalpadding blocksize=(512,512,1)');
  UPDATE georaster_table SET georaster=gr WHERE id=12;
END;
/

In the real world, however, the source images are often collected under different circumstances so as to have different resolutions or large areas of overlap, or using a different georeference system. In such cases, you can use the SDO_GEOR_AGGR.mosaicSubset procedure to mosaic these source images into one uniform mosaicked image. Compared to SDO_GEOR.mosaic, the SDO_GEOR_AGGR.mosaicSubset procedure provides more features and options:

• The source images do not have to be in the same coordinate system (SRID) and have the same georeferencing information or resolutions.

• The source images can be mosaicked on a user-specified pyramid level.

• The source images can be mosaicked on user-specified bands.

• The output images can have a different coordinate system and resolution than the input images (outSRID and outResolutions parameters).

• You have more control on the output of the overlapping area through the mosaicParam parameter: commonPointRule can specify which cell value to use for the output at the overlapping area, and NODATA can indicate whether to consider the NODATA value at the overlapping area.

• The output mosaicked image can be aligned at a specified point (the reference point). The source image can be resampled in order to align with the reference point if the source image is out of alignment more than the resampleTolerance value specified in mosaicParam.

• If there is small gap between the source images that is less than 2 pixels wide, it can be filled using the neighboring pixel values when fillGap is true in mosaicParam.

• Limited color balancing (linear stretching and normalization) is supported.

• Parallel processing is supported to speed up the mosaicking process.

Example 6-19 SDO_GEOR_AGGR.mosaicSubset

Example 6-19 uses SDO_GEOR_AGGR.mosaicSubset to mosaic all the source images from two GeoRaster tables (GRTAB1 and GRTAB2) into a large mosaicked image in SRID 4326 with a resolution of 30 meters on the x and y dimensions.

DECLARE
  resolutions sdo_number_array;
  gr sdo_georaster;
BEGIN
    insert into georaster_table (georid, georaster)
        values (10, sdo_geor.init('RDT_1',10))
         returning georaster into gr;
 
    resolutions := sdo_number_array(30, 30);
    sdo_geor_aggr.mosaicSubset('grtab1, grtab2', 'grobj, grobj',
                               0, 4326, null, null, null,
                               null, null, null, resolutions, 'unit=meter',
                               'commonPointRule = end, nodata=true, resampleTolerance=0.2, resampling=bilinear, fillGap=true',
                               'blocking=optimalpadding blocksize=(512, 512, 3)', gr, null, 'parallel=4');
 
      update georaster_table set georaster = gr where georid=10;
      commit;                                                              
END;
/

In Example 6-19:

• Any source image that is not rectified is rectified; any source image that is not in SRID 4326 is reprojected to SRID 4326.

• Any source image that has a resolution other than 30 meters is scaled to a resolution of 30 meters.

• The nodata keyword in the mosaicParam parameter is specified as true, which means the NODATA values in the overlapping area are not considered.

• The resampleTolerance keyword in the mosaicParam parameter is specified as 0.2, which means that if the source image is offset from the target by more than 0.2 pixel, the source image is resampled.

• The resampling method is specified as bilinear in the mosaicParam parameter.

• The degree of parallelism is specified as 4 in the parallelParam parameter.

You can call SDO_GEOR_AGGR.validateForMosaicSubset before calling SDO_GEOR_AGGR.getMosaicSubset to make sure that the source images can be mosaicked.