[gdal-dev] Gdal_grid produces a lot of empty cells

[Paul Meems]

Hello List,

I have a CSV-file with sensor data and GPS coordinates in WGS84.
I convert these GPS coordinates to a projected projection. Mostly a UTM
zone, depending on the location.
These X and Y values are added to the CSV-file and saved with a different
name.
Then I create a VRT-file to get the proper Z-column:
<OGRVRTDataSource>
    <OGRVRTLayer name="01 field data-EPSG32634">
        <SrcDataSource>01 field data-EPSG32634.csv</SrcDataSource>
        <GeometryType>wkbPoint</GeometryType>
        <GeometryField encoding="PointFromColumns" x="X" y="Y"
z="Countrate"/>
    </OGRVRTLayer>
</OGRVRTDataSource>

This VRT file is sent to Gdal_grid to create a tiff-file using these
arguments:
-a invdist:power=2:
   smoothing=20:
   min_points=16:
   max_points=56:
   radius1=40:
   radius2=40:
   nodata=1.701410000000000071e+38
-txe 273360 274940
-tye 5559440 5561020
-tr 12.5 12.5
-a_srs EPSG:32634
"01 field data-EPSG32634-TC.vrt" "tmpmrmq4d.tiff"

This tiff-file looks like this: https://pasteboard.co/tKP1ZJ1cAY2r.png

This process works fine and we've been using this for several years now.

Recently we added a process step to export one of the sensor data to a
PostGIS database and with the help of Geoserver and OpenLayers we've
created a very simple webpage.
In this step, things are getting weird. The WebMap is in EPSG:3857
(Pseudo/Web-Mercator).
I start with the original CSV-file with the sensor data and add the X and Y
columns again, now in EPSG:3857. Then I create a similar vrt-file:
<OGRVRTDataSource>
    <OGRVRTLayer name="01 field data-EPSG3857">
        <SrcDataSource>01 field data-EPSG3857.csv</SrcDataSource>
        <GeometryType>wkbPoint</GeometryType>
        <GeometryField encoding="PointFromColumns" x="X" y="Y"
z="Countrate"/>
    </OGRVRTLayer>
</OGRVRTDataSource>

And call gdal_grid with similar arguments:
-a invdist:power=2:
   smoothing=20:
   min_points=16:
   max_points=56:
   radius1=40:
   radius2=40:
   nodata=1.701410000000000071e+38
-txe 1984580 1987000
-tye 6471320 6473740
-tr 12.5 12.5
-a_srs EPSG:3857
"01 field data-EPSG3857-TC.vrt" "tmpun25fo.tiff"

But now the result is: https://pasteboard.co/bcMp7rP5eyff.png

A huge amount of empty cells. I don't understand this.
Both EPSG:32634 (WGS 84 / UTM zone 34N) and EPSG:3857 are in meters so I
thought the same settings should work.
But in fact only with very different settings I can get a proper grid:
 -a invdist:power=2:
   smoothing=60:
   min_points=4:
   max_points=56:
   radius1=80:
   radius2=80:
   nodata=1.701410000000000071e+38
-txe 1984580 1987000
-tye 6471320 6473740
-tr 20 20
-a_srs EPSG:3857
"01 field data-EPSG3857-TC.vrt" "tmpsvlhgy.tiff"

I got these values through trial and error.
https://pasteboard.co/7pypGsynS09K.png

But now we loose a lot of variation because the search radius and grid cell
size are much higher.

To make it even more weird. Other fields nearby don't have this problem and
I also have this problem with other projections like:

Belgium (EPSG:3857): https://pasteboard.co/ZXLTsEuvpZ3s.png
Belgium (EPSG:31370): https://pasteboard.co/yNpF6yujF3tV.png

The Netherlands (EPSG:3857): https://pasteboard.co/8PpDJrDi1Wxq.png
The Netherlands (EPSG:28992): https://pasteboard.co/38SlnHVCog4m.png

Norway (EPSG:3857): https://pasteboard.co/q69jQHRK40c6.png
Norway (EPSG: 25833): https://pasteboard.co/gs9HNAmnaxSJ.png

We do have enough data points: https://pasteboard.co/4CmBs7kgQGVO.png
https://pasteboard.co/J8ZbMHKeErZu.png

Thanks if you made it all the way here ;)

My questions are:

   - Am I using Gdal_Grid in the proper way or do I need to change one or
   more arguments?
   - What might be an explanation for the empty cells when using EPSG:3857?


Thanks,

Paul Meems
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/gdal-dev/attachments/20240522/57c738fb/attachment.htm>