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

Wed May 22 11:34:41 PDT 2024

Hi,

I forgot to say that I think it is not a good idea to use gdal_grid with EPSG:3857 at all. The scale of Web Mercator is correct only at the equator. In your data from Poland, if the real ellipsoidal distance between two points on the ground is 10 meters, the cartesian distance in EPSG:3857 claims that the distance is 15 meters. I think it would be better to compute the grid in a CRS that has approximately correct scale, and re-project the raster into Web Mercator if needed. That said, perhaps there is a bug in the gdal_grid utility.

-Jukka Rahkonen-

Lähettäjä: Paul Meems <bontepaarden at gmail.com>
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/5407685e/attachment-0001.htm>