spatialite_dem
Not logged in

Back to main SpatiaLite-Tools Wiki page


spatialite_dem - Introduction

Goals of Tool spatialite_dem

Create a Dem-Database from 1 or more xyz-files

  1. this can be a single xyz-file

  2. a directory of xyz-files

  3. a list-file, each line containing one xyz-file-name,

    whereby list file must be in the directory that contains the listed files

Notes: Before importing, the first line of each file will read

, checking if 3 double values can be properly converted
  • if yes: a valid xyz-file is assumed and the file-name and point x/y are stored in a memory table

During importing, this memory table is read, sorting by point_y ASC, point_x ASC

Goal is to import all points sorted as:
  • y='South to North' and x='West to East'

The created table will contain the x/y/z values as doubles and a POINTZ geometry with the given srid.

Then a RecoverGeometryColumn and CreateSpatialIndex is done. For my 1.116 billion Berlin-Database
  • importing took about 8 hours (the same as the sqlite3 csv import of the same data) to compleate and
  • creating the SpatialIndex took a further 2 days and 8 hours.

Sort command for a xyz-file:

In theory;
  • you should recieve your .xyz files properly sorted
In praxis;
  • this is not always the case (in my case, out of 279 files, at least 1 was not correctly sorted)

For Linux, a correct usage would be:
  • sort -n -k2 -k1 392_5810.xyz -o 392_5810.sort.xyz
    • -n, --numeric-sort : Compare according to string numerical value.
    • -k2 : Start at the second column (y-Value).
    • -k1 : End at the first column (x-Value).
    • input_file.xyz : the input file to read.
    • -o : the output file to create. must not be the same as the input file

For gdal, it is crucial that the .xyz- file is sorted in the Y-Direction
    gdal_translate -ot Float32 -a_srs EPSG:25833 392_5810.xyz 2007.392_5810.dhhn92.25833.tif
      ERROR 1: Ungridded dataset: At line 348, change of Y direction

where as the correctly sorted input.xyz, worked correctly
    gdal_translate -ot Float32 -a_srs EPSG:25833 392_5810.sort.xyz 2007.392_5810.dhhn92.25833.tif

For spatialite_dem, it is only cosmetic
    no attempt will be made by spatialite_dem to sort the data inside input file



Query of Dem-Database (-fetch_yx)

This is very swift, mostly around 3 milli-seconds
spatialite_dem -v  -fetchz_xy  24700.55278283251 20674.74537357586
Dem: srid 25833
Dem: extent min x/y(370000.0000000,5798000.0000000)
	    max x/y(415999.0000000,5837999.0000000)
Dem: extent width(45999.0000000)
	   height(39999.0000000)
Dem: rows_count(utm_point) 1116000000
Dem: resolution(utm_point) 1.6486685
Dem: geometry_type(1001) has_z1 has_m0
Dem: spatial_index_enabled1
Dem '/home/mj10777/000_links/build_berlin_data/build.berlin_admin/geometry_admin/source_db/2007.berlin.dhhn92.db'
 TABLEberlin_dhhn92_2007 with GEOMETRY-Columnutm_point
-W-> -rdem was set. Using: resolution(0.5000000), overriding the calculated value: 1.6486685
FetchZ modus: with default_srid[3068]  x24700.5527828 y20674.7453736 has_m0
FetchZ modus: with     dem_srid[25833] x391427.4994515 y5819297.5054299 z33.5600000
>> time needed: 07 milli-secs

Note:

without '-v'

  • only the z-value is returned and can be used in a bash script.


Configuration of default Dem-Database to use and default srid for queries

can be saved with -save_conf to file default 'spatialite-dem.conf'
With:


spatialite_dem will read that file first, setting the default values

The Database I created was called '2007.berlin.dhhn92.db' and placed in the '/long/path/to/file/' directory.

Using the -sniff parameter, the settings can be checked:
spatialite_dem  -sniff -d /long/path/to/file/berlin_admin_geometries.db -t berlin_street_segments -g soldner_segment -ddem /long/path/to/file/2007.berlin.dhhn92.db -tdem berlin_dhhn92_2007 -gdem utm_point

resulting also with

Source: srid 3068 ... Source Database: has passed all checks. ... Dem: srid 25833 ... Dem: resolution(utm_point) 1.6486685 ... Dem Database: has passed all checks. Sniffing modus: All pre-conditions have been fulfilled. to start update, use the '-updatez' parameter. to save dem-conf, use the '-save_conf' parameter.

With SPATIALITE_DEM having been set, using -save_conf, would save the config to the file pointed to in SPATIALITE_DEM.

Since I will use this Database for 99% of my needs, I have added the export command to my local .bashrc.

So when calling (with out -v):
spatialite_dem -fetchz_xy  24700.55278283251 20674.74537357586
33.5600000
with out having to give:

over and over again.

Important is setting/editing the dem_resolution/-rdem parameter in the conf file

Since I forgot to use the '-rdem' parameter when using '-save_conf':
-rdem 0.50

this would then have to be changed in the file:

# -- -- ---------------------------------- --
# Area around given point to Query Dem-Geometry in units of Dem-Srid
# -> Rule: a Dem with 1m resolution: min=0.50
dem_resolution=0.50
# -- -- ---------------------------------- --

In this case the Database contains only the Open-Data points of Berlin,
    for the area of Brandenburg (which has no OpenData policy) there are no points, but are, never the less, inside the Extent.

dem_berlin_area.png
    So although the resolution is 1.0 meters, where dem points exist, the calculated resolution for the extent is: 1.6486685.
    For this reason the -rdem parameter, overriding the calculated resolution, should always be used.

Notes:

-v
  • will give out more information about what the tool is doing
    • and should not be used in a bash script
-sniff
  • will simulate what the tool will do, where possible
    • useful to check if the given parameters are correct
-save_conf
  • will save the used configuration to the given file
    • useful to the usage of the tool, using an alternative defaut parameter values
-rdem
  • will override the calculated resolution
    • which, in most cases, will be needed
  • Area around a given point to search for: 1.0/2=0.5


Updating a Database SpatialTables with Z-Values from Dem-Database (-updatez)

Preconditions:

That a copy from an origianl TABLE has been made that contains a Geometry of the Z or ZM Geometry type (using CastToXYZ or CastToXYZM).
  • The the Z-Value of the point(s) must be 0.
    • all Geometry-Types are supported: POINT, LINESTRING, POLYGON, GEOMETRYCOLLECTION
      including all MULTI-Types of POINT, LINESTRING and POLYGON
  • The srid of the Geometry can be different that the srid of the Dem-Points
    • the Z-Values will be taken over without any 'height' transformation
The update Process will not proceed if any preconditiions are not fulfilled.

-sniff
  • can be used to check that these pre-conditions are fulfiled

A Sql-Script to convert one TABLE from an XY-based Database to XYZ-based Database could look like this:
With the Database berlin_admin_geometries.xyz.db
Third we will remove the -v command, so that -updatez without any information messages
 spatialite_dem -updatez --db-path berlin_admin_geometries.xyz.db -t berlin_polygons -g soldner_ring
 spatialite_dem -updatez --db-path berlin_admin_geometries.xyz.db -t berlin_polygons -g soldner_linestring
 spatialite_dem -updatez --db-path berlin_admin_geometries.xyz.db -t berlin_polygons -g soldner_center

(As one can see, one sees nothing)

Process:

If the srid of the Geometry is different that the srid of the Dem-Points
  • The Geometry will be transformed, as a copy, to the srid of the Dem-Points.
    • the found Z-Point will be replaced in the original Geometry without any 'height' transformation
For any Geometry where the Z-Value is <> 0
  • the Geometry will be ignored, assuming the value has already been set
    • allowing for a swift update to added Geometries
For any Geometry where the Z-Value has not changed
  • the Geometry will not be updated
For any Geometry that contain a M-Value
  • the M-Value will remain unchanged

Notes:

Return code of tool:
  • 0 correct
  • 1 error
    • as expected on a unix system (bash)
The (internal) functions return:
  • 1 correct
  • 0 error
    • as expected in a c based program
Depending on what is needed. different running settings can be used
  • Logical checks with Messages, but without running the command
  • Messages while running
  • Running only with a return code and possibly (as for -fetchz) with a result
    • for use with bash or other calling applications



Sample Photos

Sample Points
The Database used here used srid 25833 (ETRS89 / UTM zone 33N)
but displays the map with srid 3068 (DHDN / Soldner Berlin)

    Sample Points

  • SRID=3068;POINT(24700.55278283251 20674.74537357586)
  • SRID=25833;POINT(391427.4990121186 5819297.504889704)

id_demx_utmy_utmZx_soldnery_soldnerdistance_meters
630745428391427.0000005819298.00000033.56000024700.04396020675.2302040.702887
630745429391428.0000005819298.00000033.60000024701.04403320675.2492300.703667
630709428391427.0000005819297.00000033.62000024700.06298620674.230126 0.710572
630709429391428.0000005819297.00000033.61000024701.06306020674.2491530.711344

Sql used to create the sample points

Sql used to create the above list
    SELECT
      id_dem, point_x AS x_utm, point_y AS y_utm,
      point_z AS Z,
      ST_X(ST_Transform(utm_point,3068)) AS x_soldner,
      ST_Y(ST_Transform(utm_point,3068)) AS x_soldner,
      ST_Distance(utm_point,ST_Transform(MakePoint(24700.55278283251,20674.74537357586,3068),25833)) AS distance_meters
    FROM "berlin_dhhn92_2007"
    WHERE
    (
     -- First condition
     ( -- Out of the 1.116.000.000 records, select only those that are in range 4 will be returned
      -- for a SpatialView, the Primary-Key must be used, otherwise ROWID can be used
      id_dem IN
      (
       -- Use the Spatialite internal logic to simplify the SQL-Query, avoiding rounding errors that could miss a valid geometry
       SELECT ROWID FROM SpatialIndex
       WHERE
       (
        -- Use the created index for the given TABLE
        -- To query an ATTACHED Database, replace 'main' with the used schema-name
        -- > 'DB=main.' for a non-ATTACHED Database is not mandatory
        (f_table_name = 'DB=main.berlin_dhhn92_2007') AND
        -- Use the given Geometry-Column mandatory only where there is more than 1 Geometry-Column
        (f_geometry_column = 'utm_point') AND
        -- search an area 0.5 meters around the given point (= 1 meter width/height)
        (search_frame = ST_Buffer(ST_Transform(MakePoint(24700.55278283251,20674.74537357586,3068),25833),0.5))
       )
      )
     )
    )
    ORDER BY distance;
    

    -- spatialite_gui: 4 rows fetched in 00:00:00:079 from 1.116 billion records



Screenshot showing the test points (Scale 1:15)

street_scale_15.png
Agua Point : 24700.55278283251 20674.74537357586

Since the SPATIALITE_DEM is set to a config file where default_srid=3068 and -rdem=0.5 is set, a simplified form of the command can be used:

spatialite_dem -fetchz_xy 24700.55278283251 20674.74537357586
33.5600000
(Upper-Left green circle)

spatialite_dem -fetchz_xy 24700.554 20674.745
33.6000000
(Upper-Right green circle)

spatialite_dem -fetchz_xy 24700.55 20674.70
33.6200000
(Lower-Left green circle)

spatialite_dem -fetchz_xy 24700.60 20674.70
33.6100000
(Lower-Right green circle)

This street intersection is not a cobble stoned street, so where does the 6 cm difference come from?


Screenshot showing the test points (Scale 1:100)

street_scale_100-25833.png

Determin if DEM-Points reflect the heights of buildings or ground level


spatialite_dem --save_conf

    # -- -- ---------------------------------- --
    # For use with spatialite_dem
    # -- -- ---------------------------------- --
    # export SPATIALITE_DEM/long/path/to/file/2007.berlin.dhhn92.conf
    # -- -- ---------------------------------- --
    # Full path to Spatialite-Database containing a Dem-POINTZ (or POINTZM) Geometry
    dem_path/long/path/to/file/2007.berlin.dhhn92.db
    # Table-Name containing a Dem-POINTZ (or POINTZM) Geometry
    dem_table=berlin_dhhn92_2007
    # Geometry-Column containing a Dem-POINTZ (or POINTZM) Geometry
    dem_geometry=utm_point
    # Srid of  Dem-Geometry
    dem_srid=25833
    # -- -- ---------------------------------- --
    # Area around given point to Query Dem-Geometry in units of Dem-Srid
    # -> Rule: a Dem with 1m resolution: min=0.50
    dem_resolution=0.50
    # -- -- ---------------------------------- --
    # Default Srid of Queries against Dem-Geometry -fetchz_xy, -updatez
    default_srid=3068
    # -- -- ---------------------------------- --
    # Count of rows in Dem-Geometry
    dem_rows_count=1116000000
    # Min X of Dem-Geometry
    dem_extent_minx=370000.0000000
    # Max X of Dem-Geometry
    dem_extent_maxx=415999.0000000
    # Min Y of Dem-Geometry
    dem_extent_miny=5798000.0000000
    # Max Y of Dem-Geometry
    dem_extent_maxy=5837999.0000000
    # Width of Dem-Area in Srid-Units
    dem_extent_width=45999.0000000
    # Height of Dem-Area in Srid-Units
    dem_extent_height=39999.0000000
    # -- -- ---------------------------------- --
    

Process:

The main goal of the conf file is to simplify the usage of the tool
  • The assumption is that main Dem-Database will often be used
    • therefore the Dem-Input parameters, that are needed for all command forms, should be automated
  • Also the default srid to use and the Radius (area around a given point) should be automated
Creation of conf file:
  • First create a valid command with all needed parameters:

    • spatialite_dem -ddem /home/mj10777/000_links/build_berlin_data/build.berlin_admin/geometry_admin/source_db/2007.berlin.dhhn92.db \
       -tdem berlin_dhhn92_2007 \
       -gdem utm_point \
       -rdem 0.50 \
       -default_srid 3068 \
       -sniff
      

      Dem: srid 25833 Dem: default_srid3068 Dem: extent min x/y(370000.0000000,5798000.0000000) max x/y(415999.0000000,5837999.0000000) Dem: extent width(45999.0000000) height(39999.0000000) Dem: rows_count(utm_point) 1116000000 Dem: resolution(utm_point) 1.6486685 Dem: geometry_type(1001) has_z1 has_m0 Dem: spatial_index_enabled1 Dem '/home/mj10777/000_links/build_berlin_data/build.berlin_admin/geometry_admin/source_db/2007.berlin.dhhn92.db' TABLEberlin_dhhn92_2007 with GEOMETRY-Columnutm_point -W-> -rdem was set. Using: resolution(0.5000000), overriding the calculated value: 1.6486685


  • If the result is correct:
    • Set the Environment variable SPATIALITE_DEM to a valid path/filename of your choice:
      • such as /long/path/to/file/database-name.conf
    • Add a -save_conf to the command and run:
    • Add the export SPATIALITE_DEM= ... to your local bashrc if desired.
    Inside a bash script
  • You can switch from one Dem-Database to another by setting SPATIALITE_DEM to another, existing, conf file.


spatialite_dem --help

    usage: spatialite_dem ARGLIST
    ==============================================================
    -h or --help                    print this help message
    ========================== Parameters ========================
      -- -- ---------------- Dem-Data Database ---------------- --
    -ddem or --dem-path  pathname to the SpatiaLite Dem DB
    -tdem or --table-dem table_name SpatialTable or SpatialView
    -gdem or --geometry-dem-column col_name the Geometry column
    	 must be a POINT Z or a POINT ZM type
    -rdem or --dem-resolution of the dem points while searching
    	 the automatic resolution calculation is based on the row_count
    	 within the extent, which may not be correct!
    	 Use '-rdem' to set a realistic value
    
    -- -- -------------- Source-Update-Database ----------------- -- -d or --db-path pathname to the SpatiaLite DB -t or --table table_name, must be a SpatialTable -g or --geometry-column the Geometry column to update must be a Z or a ZM Dimension type use CastToXYZ(geom) or CastToXYZM(geom) to convert -- -- --------------- General Parameters ---------------- -- -mdem or --copy-m 0=no, 1= yes [default if exists] -default_srid or --srid for use with -fetchz -fetchz_xy x- and y-value for use with -fetchz -v or --verbose messages during -updatez and -fetchz -save_conf based on active -ddem , -tdem, -gdem and -srid when valid
    -- -- -------------------- Notes: ---------------------- -- -I-> the Z value will be copied from the nearest point found -I-> the Srid of the source Geometry and the Dem-POINT can be different -I-> when -fetchz_xy is used in a bash script, -v should not be used the z-value will then be retured as the result
    -- -- -------------------- Conf file: ------------------- -- -I-> if 'SPATIALITE_DEM' is set with the path to a file -I--> 'export SPATIALITE_DEM=/long/path/to/file/berlin_dhh92.conf' -I-> then '-save_conf' save the config to that file -I-> this file will be read on each application start, setting those values -I--> the parameters for : which Dem-Database and Geometry and the default_srid to use for queries -> would then not be needed
    -- -- ---------------- Importing .xyz files: ------------------- -- -I-> a single xyz.file or a directory containing .xyz files can be given for directories: only files with the extension .xyz will be searched for -I-> a single list.file inside a directory containing .xyz files can be given each line containing the file-name that must exist in that directory -I-> validty checks are done before importing xyz-files the first line may contain only 3 double values (point_x/y/z) if valid, the file-name and the point_x/y points are stored when importing, the list will be read based of the y/x points
    -- -- ---------------- Sorting .xyz files ---------------------- -- -I-> xyz.files should be sorted: y='South to North' and x='West to East': sort -n -k2 -k1 input_file.xyz -o output_file.sort.xyz =========================== Commands =========================== -sniff default analyse settings without UPDATE of z-values -updatez Perform UPDATE of z-values -fetchz Perform Query of z-values using -fetchz_x_y and default_srid will be assumed when using -fetchz_x_y -create_dem create Dem-Database using -ddem,-tdem, -gdem and -srid for the Database -d as a dem.xyz file -import_xyz import another .xyz file into a Dem-Database created with -create_dem these points will not be sorted, but added to the end =========================== Sample =========================== --> with 'SPATIALITE_DEM' set: spatialite_dem -fetchz_xy 24700.55278283251 20674.74537357586 33.5600000 ==============================================================



Back to main SpatiaLite-Tools Wiki page