<html>
  <head>

    <meta http-equiv="content-type" content="text/html; charset=UTF-8">
  </head>
  <body>
    <p>Dear all Users and/or Participants of PostGIS,</p>
    <p>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.<br>
    </p>
    <p>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<br>
    </p>
    <p>I stored both in a PostgreSQL-DB with PostGIS extension and a
      tiling size of 40x40.<br>
    </p>
    <p>The target to create polygons and request raster data of both
      layers mentioned above covered by those polygons.<br>
      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. <br>
      Afterwards I created a query for the covered raster data:</p>
    <pre>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;

</pre>
    <p>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.<br>
      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.<br>
      <br>
      By the way, does someone have a hint for me, to further accelerate
      this query? <br>
    </p>
    <p><br>
      System setting:<br>
          DB: PostgreSQL version 12.1 <br>
          OS: Ubuntu 18.04<br>
          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</p>
    <p>Thanks in advance.</p>
    <p>Best Regards <br>
    </p>
    <p><br>
    </p>
    <p><br>
    </p>
  </body>
</html>