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

Wed May 22 10:09:15 PDT 2024

Hi,

I do not know what happens, but it seems that the algorithm finds rather few points within the search ellipse. The radius of the circle is 40 m and by looking at the data there are typically at least 30 points within that distance. However, for filling the whole raster I had to use "min_points=4".
I could repeat the behavior by opening the 3857 data into QGIS and running the invdist from there. And when I reprojected the point data into EPSG:32634 then QGIS created a good gridded raster. The difference in scale between EPSG:3857 and EPSG:32634 is not so big in Poland that it could be the reason IMHO.
For removing possible errors due to CSV format and VRT I materialized the datasets into OpenJUMP JML format. It could have been any other GIS format.

Unfortunately I cannot help you more. If I were a programmer I think I would try to catch how many points the algorithm finds within the search radius in EPSG:3857, compare how many points it finds from EPSG:32634 data, and finally consider if the results make sense.

-Jukka Rahkonen-

Lähettäjä: Paul Meems <bontepaarden at gmail.c
Lähetetty: keskiviikko 22. toukokuuta 2024 18.24
Vastaanottaja: gdal-dev at lists.osgeo.org
Kopio: Rahkonen Jukka <jukka.rahkonen at maanmittauslaitos.fi>
Aihe: Re: [gdal-dev] Gdal_grid produces a lot of empty cells

Sorry Jukka,

Forgot all about that ;)

I've attached the csv from the field in Poland with the original lon-lat values and the X and Y in EPSG:3857

Op wo 22 mei 2024 om 16:57 schreef Rahkonen Jukka <jukka.rahkonen at maanmittauslaitos.fi<mailto:jukka.rahkonen at maanmittauslaitos.fi>>:
Hi,

At least I cannot resolve your issue on paper. Please share some test data.

-Jukka-

Lähettäjä: gdal-dev <gdal-dev-bounces at lists.osgeo.org<mailto:gdal-dev-bounces at lists.osgeo.org>> Puolesta Paul Meems via gdal-dev
Lähetetty: keskiviikko 22. toukokuuta 2024 17.31
Vastaanottaja: gdal-dev at lists.osgeo.org<mailto:gdal-dev at lists.osgeo.org>
Aihe: [gdal-dev] Gdal_grid produces a lot of empty cells

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/9e1516de/attachment-0001.htm>