SAP HANA faces Coastline Paradox

[Vitaliy-R]

It is finally time to get back to last month's Get to know your neighbors with SAP HANA and to the comment from db48526ccb71414183a64aa46e31e784:

“Yeah, these measurements are all significantly shorter than the official border lengths, but in the right order…”

Let's compare results from SAP HANA to information from Wikipedia for Germany:



So, maybe the problem is with the scale of open data I loaded into SAP HANA, and so the area of Germany will be as well proportionally smaller comparing to the official number of 357,022 km²? Let's check:

select round(a.shape.ST_SRID(4326).ST_Area('kilometer')) as "area"

from "TESTSGEO"."cntry00"  a

where a.CNTRY_NAME like 'Germany';



--Result in square kilometers is 357,221

Hmmm, for the area the difference between results is only 0,0557%. So it is not a scale problem.



To better understand what's going on, let's get coastlines zip file from Generalized Coastlines  datasets available from OpenStreetMap website:



When uncompressed we get sets of shapefiles at 8 different zoom levels.

hxeadm@hxehost:/usr/sap/HXE/HDB90/work> mkdir osm_coastlines

hxeadm@hxehost:/usr/sap/HXE/HDB90/work> cd osm_coastlines/

hxeadm@hxehost:/usr/sap/HXE/HDB90/work/osm_coastlines> wget http://data.openstreetmapdata.com/coastlines-generalized-3857.zip

hxeadm@hxehost:/usr/sap/HXE/HDB90/work/osm_coastlines> unzip coastlines-generalized-3857.zip 

hxeadm@hxehost:/usr/sap/HXE/HDB90/work/osm_coastlines> cd coastlines-generalized-3857/

hxeadm@hxehost:/usr/sap/HXE/HDB90/work/osm_coastlines/coastlines-generalized-3857> ls -l *.shp

-rw-r--r-- 1 hxeadm sapsys   172884 Mar 26 11:02 coastlines_z1.shp

-rw-r--r-- 1 hxeadm sapsys   329740 Mar 26 11:02 coastlines_z2.shp

-rw-r--r-- 1 hxeadm sapsys   868932 Mar 26 11:02 coastlines_z3.shp

-rw-r--r-- 1 hxeadm sapsys  2051324 Mar 26 11:02 coastlines_z4.shp

-rw-r--r-- 1 hxeadm sapsys  4842884 Mar 26 11:02 coastlines_z5.shp

-rw-r--r-- 1 hxeadm sapsys 10507132 Mar 26 11:02 coastlines_z6.shp

-rw-r--r-- 1 hxeadm sapsys 19768956 Mar 26 11:02 coastlines_z7.shp

-rw-r--r-- 1 hxeadm sapsys 46651260 Mar 26 11:02 coastlines_z8.shp

Already by looking at size of files we can see that z8 zoom should contains much more details comparing to previous files.

Let's load these files. Please note shapefiles are created using Web Mercator projection, i.e. spatial reference system called EPSG:3857. This reference system is not available out-of-the-box in SAP HANA Express. Create this SRS first, if missing in "PUBLIC"."ST_SPATIAL_REFERENCE_SYSTEMS".

--Create EPSG:3857 SRS

CREATE SPATIAL REFERENCE SYSTEM "WGS 84 / Pseudo-Mercator" 

IDENTIFIED BY 3857 

DEFINITION 'PROJCS["WGS 84 / Pseudo-Mercator",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Mercator_1SP"],PARAMETER["central_meridian",0],PARAMETER["scale_factor",1],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["X",EAST],AXIS["Y",NORTH],EXTENSION["PROJ4","+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext  +no_defs"],AUTHORITY["EPSG","3857"]]' 

ORGANIZATION "EPSG" IDENTIFIED BY 3857 

TRANSFORM DEFINITION '+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext  +no_defs'   

TYPE PLANAR 

COORDINATE X BETWEEN -20037508.342789248 AND 20037508.342789248 

COORDINATE Y BETWEEN -20048966.104014635 AND 20048966.104014624  

TOLERANCE DEFAULT 

SNAP TO GRID DEFAULT 

POLYGON FORMAT 'EvenOdd' STORAGE FORMAT 'Internal';



--Now load selected files

IMPORT "TESTSGEO"."coastlines_z1"

AS SHAPEFILE 

FROM '/usr/sap/HXE/HDB90/work/osm_coastlines/coastlines-generalized-3857/coastlines_z1'

WITH REPLACE SRID 3857 THREADS 4;



IMPORT "TESTSGEO"."coastlines_z4"

AS SHAPEFILE 

FROM '/usr/sap/HXE/HDB90/work/osm_coastlines/coastlines-generalized-3857/coastlines_z4'

WITH REPLACE SRID 3857 THREADS 4;



IMPORT "TESTSGEO"."coastlines_z8"

AS SHAPEFILE 

FROM '/usr/sap/HXE/HDB90/work/osm_coastlines/coastlines-generalized-3857/coastlines_z8'

WITH REPLACE SRID 3857 THREADS 4;

Now let's look at one specific case: Iceland. Wikipedia says that Iceland "lies between latitudes 63 and 68°N, and longitudes 25 and 13°W".

Select from tables only geographies within the bounding box POLYGON ((-27 63,-11 63,-11 67,-27 67,-27 63)). Once again, remember about transformations between SRS 3857 used to store geographies in coastlines tables, and SRS 4326 (GPS coordinates) used to define the bounding box.

select ST_UnionAggr("SHAPE").ST_Transform(4326).ST_asSVG()

 from "TESTSGEO"."coastlines_z1"

 where 

 "SHAPE".ST_CoveredBy(new ST_Polygon('POLYGON ((-27 63,-11 63,-11 67,-27 67,-27 63))',4326).ST_Transform(3857))=1;

 

 select ST_UnionAggr("SHAPE").ST_Transform(4326).ST_asSVG()

 from "TESTSGEO"."coastlines_z4"

 where 

 "SHAPE".ST_CoveredBy(new ST_Polygon('POLYGON ((-27 63,-11 63,-11 67,-27 67,-27 63))',4326).ST_Transform(3857))=1;

 

 select ST_UnionAggr("SHAPE").ST_Transform(4326).ST_asSVG()

 from "TESTSGEO"."coastlines_z8"

 where 

 "SHAPE".ST_CoveredBy(new ST_Polygon('POLYGON ((-27 63,-11 63,-11 67,-27 67,-27 63))',4326).ST_Transform(3857))=1;

Below are all three visualized and color coded with black for zoom details z1, red for z4, and green for z8.



Or for selected fragment:



Calculating the coastline length for these three selected zoom levels of details.

select round(ST_UnionAggr("SHAPE").ST_Transform(4326).ST_length('kilometer'))

 from "TESTSGEO"."coastlines_z1"

 where 

 "SHAPE".ST_CoveredBy(new ST_Polygon('POLYGON ((-27 63,-11 63,-11 67,-27 67,-27 63))',4326).ST_Transform(3857))=1;

--The result is 1329

 

select round(ST_UnionAggr("SHAPE").ST_Transform(4326).ST_length('kilometer'))

 from "TESTSGEO"."coastlines_z4"

 where 

 "SHAPE".ST_CoveredBy(new ST_Polygon('POLYGON ((-27 63,-11 63,-11 67,-27 67,-27 63))',4326).ST_Transform(3857))=1;

--The result is 2712

 

select round(ST_UnionAggr("SHAPE").ST_Transform(4326).ST_length('kilometer'))

 from "TESTSGEO"."coastlines_z8"

 where 

 "SHAPE".ST_CoveredBy(new ST_Polygon('POLYGON ((-27 63,-11 63,-11 67,-27 67,-27 63))',4326).ST_Transform(3857))=1;

--The result is 5001

For the same country the coastline length is between 1329 km and 5001 km depending on the granularity of details stored in the geometry! Why then to bother about shapes with lover level of details? Well, for most calculations and visualizations lower level of details may be sufficient while offering reduced memory footprint and much faster calculation. So, it is a balancing act between precision and performance.

Even measuring the real physical length of coastlines (and borderlines too) is a challenge called the Coastline Paradox, which is nicely explained in this video.



Enjoy!

---

‘Till next Tuesday, then!
-Vitaliy, aka @Sygyzmundovych