An example of “local reprojection” from WGS84 to Gauss-Boaga Rome40 with SQLite/SpatiaLite

A common problem for many Italian GIS users in archaeology is to reproject vector features from WGS84 to Gauss-Boaga Rome40 (national Spatial Reference System in Italy till 2012). For shapefiles this transformation can be performed by specific GIS tools, often based on PROJ library (http://trac.osgeo.org/proj/), by OGR utilities (i.e. ogr2ogr: http://www.gdal.org/ogr/) or others standalone software.

In SQLite/SpatiaLite the reprojection between 2 different Spatial Reference Systems (SRS) is provided by “Transform” function (for complete description and syntax see: http://www.gaia-gis.it/gaia-sins/spatialite-tutorial-2.3.1.html#t5): it’s possible to convert data in all the many Spatial Reference Systems stored in SpatiaLite’s “spatial_ref_sys” table and coded by SRID (Spatial Reference Identifier) and EPSG dataset. For example, to transform a vector from WGS84 to UTM ED50 it's enough to update geometry field in your table with the following SQL command:

UPDATE your_table SET geometry_field = Transform(geomWGS84, 23032);

where “geomWGS84” is a geometry field in WGS84 system (EPSG:4326) and “23032” is EPSG value for UTM ED50 zone 32 N.

The problem is that the reprojection between global systems often yields an approximate and imprecise result. For reaching a more accurate outcome a “local reprojection” is required: it is feasible using specific parameters of transformation (translation, rotation and scaling), different for each part of territory.

My example is about Veneto Region in Northern Italy. I needed to reproject some points, representing archaeological sites of this Region, from WGS84 to Gauss-Boaga Rome40 – West (EPSG:3003). I tested this transformation in SpatiaLite simply using the Spatial Reference Systems stored in spatial_ref_sys table; so I executed the SQL command

UPDATE my_table SET geometry_field = Transform(geom4326, 3003);

where “geom4326” is the geometry field recording the geometries of my points in WGS84 and “3003” is the EPSG value for Gauss-Boaga Rome40 W. The outcome was not good: the points converted were located till 80 meters far away from the correct position.

In order to reduce this displacement, I read this post of Flavio Rigolon about reprojection with ogr2ogr and I adapted that solution for SQLite/SpatiaLite (but I think it could work also with PostGIS in similar way).

In particular, I added a new Spatial Reference System (SRS) in spatial_ref_sys table:

INSERT INTO spatial_ref_sys VALUES(30033003,'epsg',30033003,'Veneto 3003','+proj=tmerc +ellps=intl +lat_0=0 +lon_0=9 +k=0.999600 +x_0=1500000 +y_0=0 +units=m +towgs84=-104.1,-49.1,-9.9,0.971,-2.917,0.714,-11.68','');

My new SRS is identified by SRID and EPSG value “30033003”, is called “Veneto 3003” and has the same geodetic attributes (projection, ellipsoid, units, etc.) of Gauss-Boaga Rome40, but with the addition of 7 trasformation parameters (translation + rotation + scaling) defined by attribute “towgs84”.

For testing precision of my new SRS I selected 5 points of Veneto Region identifiable in WGS84 maps and in Gauss-Boaga Rome40 maps (CTR = “Carta Tecnica Regionale” at scale 1:10.000). I transformed my points from WGS84 SRS (EPSG:4326) to Gauss-Boaga Rome40 SRS (EPSG:3003) and to my new SRS “Veneto 3003” (EPSG:30033003). In the following images you can estimate the different positions for points transformed respectively in global EPSG:3003 and in local EPSG:30033003.

These two simple plots visualize the displacement (in meters) between correct position and points transformed.

The min, max and mean value of x (longitude) error of transformation from EPSG:4326 to EPSG:3003 is: 18.57, 36.26, 23.23 meters;

the min, max and mean value of y (latitude) error of transformation from EPSG:4326 to EPSG:3003 is: 63.36, 71.95, 68.74 meters;

the min, max and mean value of x (longitude) error of transformation from EPSG:4326 to EPSG:30033003 is: 0.49, 15.14, 4.31 meters;

the min, max and mean value of y (latitude) error of transformation from EPSG:4326 to EPSG:30033003 is: 3.3, 11.6, 6.8 meters.

[Part of this error is due to my test data: the 5 points in WGS84 have been selected from Googleearth and not recorded with a GPS and the same points identified in CTR are affected by the resolution of paper maps (“errore di graficismo” in Italian). I hope to do a more accurate test in the coming months...].

A mean error between 4 and 7 meters is acceptable for my purposes and in general for many archaeological works: in fact this error is not far from the best accuracy of portable GPS device (often used in archaeological surveys) and certainly smaller than positioning inaccuracy of many archaeological sites found in 19^th or in the first half of 20^th century. More accurate parameters of transformation (the 7 roto-traslation and scaling parameters of towgs84) could reduce this error, in particular in the Western and Northern part of Region where the distance from correct position seems to be greater.

That's all. If you know other (and faster) methods or if you detect mistakes in my post, please let me know. Any suggestions are welcome.

Denis Francisci

P.S. To enlarge the first image, open it in a new tab or window!