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
_______________________________________________
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev
