DIXIS Blog

We are in the process of migrating an old GIS system. For our new systems we use POSTGIS. But this one still uses oracle. The data is spanning two countries: Belgium and the Netherlands. Our system does something awful: all data is stored in RD (the dutch coordinate system, using Oracle SRID 90112).

So how do we get data into the system: belgian data is entered as Lambert 72 (oracle srid 327680) and then transformed to 90112.

Our client uses a customised viewer that shows the data either in RD or Lambert72. Now we want to switch to a more generic solution, and show the data in WGS84. We are using oracle 11, so my initial process was the following

extract belgian data from tables, convert back to 327680 (SDO_CS.transform(geom, 327680))
set the SRID to 31370 (which is the correct/best srid for belgium --it has the correct transformation to wgs84) as follows: update be_geoms bg set bg.geom.sdo_srid = 31370 (so without transformation)
for dutch data I just set it to 28992
and then I transform both to WGS!

Easy! done! ready! However ... I was not ... The data was not positioned correctly. So I checked the definition in MDSYS.CS_SRS for both 28992 and 31370 and compared it to epsg.io and lo and behold: both where wrong. So now I had to update them.

Updating EPSG:31370

delete from mdsys.cs_srs where srid=31370;
Insert into MDSYS.CS_SRS (CS_NAME,SRID,AUTH_SRID,AUTH_NAME,WKTEXT,CS_BOUNDS,WKTEXT3D) values ('Belge 1972 / Belgian Lambert 72',31370,31370,'IGN Brussels www.ngi.be/html-files/french/0038.html','PROJCS["Belge 1972 / Belgian Lambert 72", GEOGCS ["Belge 1972", DATUM ["Reseau National Belge 1972 (EPSG ID 6313)", SPHEROID ["International 1924 (EPSG ID 7022)", 6378388.0, 297.0], -106.869,52.2978,-103.724,0.3366,-0.457,1.8422,-1.2747], PRIMEM ["Greenwich", 0.000000], UNIT ["Decimal Degree", 0.0174532925199433]], PROJECTION ["Lambert Conformal Conic"], PARAMETER ["Latitude_Of_Origin", 90.0], PARAMETER ["Central_Meridian", 4.3674866666666667], PARAMETER ["Standard_Parallel_1", 51.1666672333333333], PARAMETER ["Standard_Parallel_2", 49.8333339], PARAMETER ["False_Easting", 150000.013], PARAMETER ["False_Northing", 5400088.438], UNIT ["Meter", 1.0]]',null,'PROJCS[
  "Belge 1972 / Belgian Lambert 72",
  GEOGCS["Belge 1972",
    DATUM["Reseau National Belge 1972",
      SPHEROID[
        "International 1924",
        6378388.0,
        297.0,
        AUTHORITY["EPSG", "7022"]],
      TOWGS84[-106.869,52.2978,-103.724,0.3366,-0.457,1.8422,-1.2747],
      AUTHORITY["EPSG", "6313"]],
    PRIMEM["Greenwich", 0.000000, AUTHORITY["EPSG","8901"]],
    UNIT["degree (supplier to define representation)", 0.0174532925199433, AUTHORITY["EPSG", "9122"]],
    AXIS["Lat", NORTH],
    AXIS["Long", EAST],
    AUTHORITY["EPSG", "4313"]],
  PROJECTION ["Lambert Conformal Conic"],
  PARAMETER ["Latitude_Of_Origin", 90.0],
  PARAMETER ["Central_Meridian", 4.3674866666666667],
  PARAMETER ["Standard_Parallel_1", 51.1666672333333333],
  PARAMETER ["Standard_Parallel_2", 49.8333339],
  PARAMETER ["False_Easting", 150000.013],
  PARAMETER ["False_Northing", 5400088.438],
  UNIT["metre", 1.0, AUTHORITY["EPSG", "9001"]],
  AXIS["X", EAST],
  AXIS["Y", NORTH],
  AUTHORITY["EPSG", "31370"]]');

... and this worked and now my transformation for Lambert is correct!

Updating EPSG:28992

... proved to be a little trickier. I assumed I could just reuse the same method as for the belgian coordinate system (yes, I know, assume = ass-u-me).

I was unable to just delete or update 28992 because I got an error that a child record existed: ORA-02292 with reason COORD_OPERATION_FOREIGN_SOURCE. Googling this revealed nothing at all.

So I had to dig deeper. And deeper. Actually MDSYS.CS_SRS is actually a view which tries to update the underlying tables. And the TOWGS84 coordinates, which I had to change/update, are stored in SDO_DATUM. So after some searching, it actually proved to be quite easy. To updated the EPSG:28992, I just had to do:

update mdsys.sdo_datums set
  shift_x = 565.417,
  shift_y = 50.3319,
  shift_z = 465.552,
  rotate_x = -0.398957,
  rotate_y = 0.343988,
  rotate_z = -1.8774,
  scale_adjust = 4.0725
where datum_id = 6289;

EXECUTE SDO_CS.UPDATE_WKTS_FOR_EPSG_DATUM(6289);

My first initial (naive) assumption was that the SDO_CS.UPDATE_... functions would actually retrieve the latest EPSG definitions, unfortunately no such luck :) :)

Stuff like this makes me appreciate PostGIS even more.

When storing spatial data into Oracle, there are a few steps one needs to complete (as shown here too):

insert the data (obviously)
update the USER_SDO_GEOM_METADATA table. This table specifies for each column the bounding box and the SRID (coordinate system).
create a spatial index I have created a script to this automatically for me, once all tables are filled with their data (e.g. after import or after using FME to import your data). [sourcecode language="sql"] set serverout on DECLARE schema_orig varchar2(100) := upper('&user_orig'); CURSOR ctab is select TABLE_NAME from sys.dba_tables where owner = schema_orig; tn sys.dba_tab_columns.TABLE_NAME%TYPE ; CURSOR ctabcol is select column_name,DATA_TYPE from sys.dba_tab_columns where owner = schema_orig and table_name = tn ; collist varchar2(2000) ; query varchar2(2000) ; column_name varchar2(200); has_geometric_column boolean; col_count number; BEGIN -- dbms_output.enable(1000000); FOR ctabrec IN ctab LOOP tn := ctabrec.table_name ; has_geometric_column := false; col_count := 0; for ctabcolrec IN ctabcol LOOP if ctabcolrec.data_type = 'SDO_GEOMETRY' then has_geometric_column := true; col_count := col_count + 1; column_name := ctabcolrec.column_name; -- insert data into user_sdo_geom query := 'INSERT INTO USER_SDO_GEOM_METADATA (TABLE_NAME, COLUMN_NAME, DIMINFO, SRID) VALUES ('''||tn|| ''' , '''|| COLUMN_NAME ||''', MDSYS.SDO_DIM_ARRAY( MDSYS.SDO_DIM_ELEMENT(''X'', 0, 10000000, 0.000005),MDSYS.SDO_DIM_ELEMENT(''Y'', 0, 10000000, 0.000005), MDSYS.SDO_DIM_ELEMENT(''Z'', 0,10000000, 0.000005)), &srid)'; dbms_output.put_line(query); BEGIN execute immediate query; EXCEPTION WHEN OTHERS THEN dbms_output.put_line('**** failed to insert row into user_sdo_geom_metadata for '||tn||'('||column_name||')'); --rollback; ignore any errors! there will be existing columns END; -- create spatial index query := 'CREATE INDEX '||tn||'_'||to_char(col_count)||'_SX ON '||tn||'('|| COLUMN_NAME ||') INDEXTYPE IS MDSYS.SPATIAL_INDEX'; dbms_output.put_line(query); BEGIN execute immediate query; EXCEPTION WHEN OTHERS THEN dbms_output.put_line('**** failed to create spatial index for '||tn||'('||column_name||')'); --rollback; ignore any errors END; end if; END LOOP; END LOOP; commit; END; / [/sourcecode] The script has two parameters: the schema from where to scan all tables, and the srid (coordinate system) which needs to be filled in. Any thoughts on this? suggestions?