[postgis-users] Extract specific area by ROI polygon

[klos]

Dear all Users and/or Participants of PostGIS,

I am working with a digital elevation point cloud (mean point desnity 
4-10 pts/m², positional precision +/- 30 cm, height precision +/- 15 cm, 
area: 1000x1000m). I interpolated a raster file with an resolution of 1m 
x 1m. Further denoted as DEM.

Furthermore, I have a 4 channel, high-precision, orthographic picture 
(mean positional precision is +/- 20-30 cm, horizontal resolution is 
20cm x 20cm) for the same area. Further denoted as DOP

I stored both in a PostgreSQL-DB with PostGIS extension and a tiling 
size of 40x40.

The target to create polygons and request raster data of both layers 
mentioned above covered by those polygons.
I created a polygon with an extend of e.\,g. 32 x 32m. This polygon 
become shifted over the area of interest, so that the whole area is 
covered by many polygons.
Afterwards I created a query for the covered raster data:

WITH
     ref AS(
             SELECT
             ST_Transform(rast,25832) AS align_to
             FROM geodata.dsm LIMIT 1
     ),
     aligned AS(
         SELECT
             ST_Union(
                 ST_Transform(
                         ST_Clip(
                                 dsm.rast,
                                     p.geom
                             ),
                             align_to
                     )
         ) as rast1,
             ST_Union(
                 ST_AsRaster(
                         p.geom,
                             align_to
                     )
             ) rast2
             FROM
             geodata.dsm AS dsm,
             public.classification_result_overlap_35cm_ss8 AS p
         CROSS JOIN ref
         WHERE
         p.gid=1 AND ST_Intersects(dsm.rast, p.geom)
     )
     SELECT
         ST_MapAlgebra(rast1, rast2,'[rast1]') as raster FROM aligned;

After a visual inspection of the results, it became obvious, that the 
mapping between the queried DEM and DOP is not 100% exact. It became 
more evident, when I made the polygon size smaller, lets say 8 x 8m.
I queried with ST_GeoReference the upper left coordinates of the queried 
DEM and DOP data and determined with pyproj.Geod(...).inv the distance 
between the upper left coordinates. I figured out, that the shift 
between the upper left coordinates of the DEM and DOP data ranges 
between 0 and 1.2. Might it be, that I did something wrong within my 
query? I thought maybe the max difference between DEM and DOP could 
result through the different horizontal resolution of the datasets. In 
this case, the maximum difference between the data could be \sqrt{2} m.

By the way, does someone have a hint for me, to further accelerate this 
query?


System setting:
     DB: PostgreSQL version 12.1
     OS: Ubuntu 18.04
     Extension: POSTGIS="3.0.0 r17983" [EXTENSION] PGSQL="120" 
GEOS="3.8.0-CAPI-1.13.1 " PROJ="6.3.0" GDAL="GDAL 3.1.0dev-437c92a187, 
released 2020/01/05" LIBXML="2.9.9" LIBJSON="0.12.1" RASTER

Thanks in advance.

Best Regards



-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/postgis-users/attachments/20200106/dd082257/attachment.html>