ASCII art map of the World: the Web Mercator projection with SAP HANA

[Vitaliy-R]

This week I decided to go back to ASCII art maps calculated with SAP HANA. It is another #GeospatialTuesday entry, and - although published on Wednesday - I still can prove it was done on Tuesday: https://twitter.com/Sygyzmundovych/status/981287806252314624 😉

So, to the point. Last year I published two posts on generating ASCII map with an equirectangular projection and with a polar azimuthal equidistant projection. Looking more closely on the first one, you could see it looks different from the maps you see regularly on the web, like OpenStreetMap, Google Maps, Bing Map etc. That's because all of them are using different map projection, called Mercator, or - to be more precise - Web Mercator projection with Spatial Reference System id 3857.



The good news: in my previous post I showed how to create SRS 3857, how to download files from OpenStreetMap-derived data and how to import them to SAP HANA using this SRS. But today we will use land polygons files (not coastlines files), so that we can build ASCII map.

Spatial Reference Systems 3857 and 1000004326 have some differences, as you can see.

select srs_id, srs_name, axis_order, min_x, max_x, min_y, max_y, linear_unit_of_measure, round_earth

from "PUBLIC"."ST_SPATIAL_REFERENCE_SYSTEMS" 

where srs_id in (3857, 1000004326);

Web Mercator's unit of measure is not a planar degree, but a meter, and therefore coordinates are different. Additionally Web Mercator is limited to bounds - 85°S to 85°N. I would probably need a separate post to discuss all these in more details, but for the moment let's take this as a given.

It means that previously used SQLScript procedure must be updated to replace some hard-coded values with parameters now. Here it is.

CREATE TYPE "NATURAL_EARTH"."TT_SHAPES" AS TABLE ( "SHAPE" ST_GEOMETRY );



CREATE PROCEDURE "NATURAL_EARTH"."ASCII_MAP_PLANAR_FLEX" 

(IN input_shape "NATURAL_EARTH"."TT_SHAPES", IN resolution INT, 

IN srs INT DEFAULT 1000004326, IN primeMeridian FLOAT DEFAULT 0.0, IN northUp INT DEFAULT 1)

 LANGUAGE SQLSCRIPT AS

 loopLon, loopLat, calcLon, calcLat FLOAT;

 minX, maxX, minY, maxY, totalX FLOAT;

 meridianShift INTEGER;

 stringSurface, pointSurface STRING;

 BEGIN

 	select "MIN_X" into minX from "PUBLIC"."ST_SPATIAL_REFERENCE_SYSTEMS" where srs_id = :srs;

 	select "MAX_X" into maxX from "PUBLIC"."ST_SPATIAL_REFERENCE_SYSTEMS" where srs_id = :srs;

 	select "MIN_Y" into minY from "PUBLIC"."ST_SPATIAL_REFERENCE_SYSTEMS" where srs_id = :srs;

 	select "MAX_Y" into maxY from "PUBLIC"."ST_SPATIAL_REFERENCE_SYSTEMS" where srs_id = :srs;

	totalX = :maxX - :minX;

	

    IF northUp <> 1 THEN northUp := -1; END IF;

    IF resolution > 1000 THEN resolution := 1000; END IF;

    IF resolution < 20 THEN resolution := 20; END IF;

    IF primeMeridian < minX OR primeMeridian > maxX THEN primeMeridian := 0.0; END IF;

    loopLat := :maxY;

    WHILE loopLat>= :minY DO

        loopLon := :minX;

        calcLat := loopLat * northUp;

        stringSurface := '';

        WHILE loopLon<= :maxX DO

            calcLon := loopLon * northUp;

            select new st_point('POINT ('||:calcLon||' '||:calcLat||')', :srs).ST_CoveredBy("SHAPE") into pointSurface 

            FROM :input_shape where "SHAPE".ST_SRID() = :srs;

            stringSurface := concat (stringSurface, pointSurface);

            loopLon := loopLon+totalX/resolution;

        END WHILE;

        meridianShift := mod(primeMeridian*northUp +totalX, totalX) * resolution/totalX;

        stringSurface := concat (right(stringSurface, resolution-meridianShift), left(stringSurface, meridianShift+1));

        INSERT INTO "NATURAL_EARTH"."DRAW_EARTH" VALUES (:calcLat, :stringSurface);

        loopLat := loopLat-totalX/resolution;

    END WHILE;

 END;

Let's create table with a single geometry of all land polygons.

CREATE COLUMN TABLE "TESTSGEO"."EARTH3857" ("SHAPE" ST_GEOMETRY(3857) ) ;

INSERT INTO "TESTSGEO"."EARTH3857" SELECT ST_UnionAggr("SHAPE") FROM "TESTSGEO"."land_polygons_z1";



SELECT "SHAPE".st_asSVG() FROM "TESTSGEO"."EARTH3857";

SVG shape looks properly.



Let's try running our procedure then!

truncate table "NATURAL_EARTH"."DRAW_EARTH";



CALL "NATURAL_EARTH"."ASCII_MAP_PLANAR_FLEX"(

	INPUT_SHAPE => "TESTSGEO"."EARTH3857",

	RESOLUTION => 120,

	SRS => 3857

);



SELECT replace(replace("SSTRING", '0', ' '), '1', '*') as "ASCIIMAP" FROM "NATURAL_EARTH"."DRAW_EARTH";

Just what we needed! Ok, I did cut a few lines at the bottom, so Antarctica is somewhat smaller on the screenshot.



But does it still work for equirectangular projection, e.g. for Australian case?

TRUNCATE TABLE "NATURAL_EARTH"."DRAW_EARTH";



CALL "NATURAL_EARTH"."ASCII_MAP_PLANAR_FLEX"(

	INPUT_SHAPE => "NATURAL_EARTH"."NE_EARTH",

	RESOLUTION => 180,

	SRS => 1000004326,

	primeMeridian => 158, 

	northUp => 0

);



SELECT replace(replace("SSTRING", '0', ' '), '1', '*') as "ASCIIMAP" FROM "NATURAL_EARTH"."DRAW_EARTH";

---

‘Till next time,
-Vitaliy, aka @Sygyzmundovych