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

2024-05-22 Thread Rahkonen Jukka via gdal-dev
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 
Lähetetty: keskiviikko 22. toukokuuta 2024 18.24
Vastaanottaja: gdal-dev@lists.osgeo.org
Kopio: Rahkonen Jukka 
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 
mailto:jukka.rahko...@maanmittauslaitos.fi>>:
Hi,

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

-Jukka-

Lähettäjä: gdal-dev 
mailto:gdal-dev-boun...@lists.osgeo.org>> 
Puolesta Paul Meems via gdal-dev
Lähetetty: keskiviikko 22. toukokuuta 2024 17.31
Vastaanottaja: gdal-dev@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:


01 field data-EPSG32634.csv
wkbPoint




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.7014100071e+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:


01 field data-EPSG3857.csv
wkbPoint




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.7014100071e+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.7014100071e+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


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

2024-05-22 Thread Rahkonen Jukka via gdal-dev
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 
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 
mailto:jukka.rahko...@maanmittauslaitos.fi>>:
Hi,

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

-Jukka-

Lähettäjä: gdal-dev 
mailto:gdal-dev-boun...@lists.osgeo.org>> 
Puolesta Paul Meems via gdal-dev
Lähetetty: keskiviikko 22. toukokuuta 2024 17.31
Vastaanottaja: gdal-dev@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:


01 field data-EPSG32634.csv
wkbPoint




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.7014100071e+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:


01 field data-EPSG3857.csv
wkbPoint




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.7014100071e+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.7014100071e+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 

Re: [gdal-dev] [EXTERNAL] [BULK] Threading and Datasets

2024-05-22 Thread Meyer, Jesse R. (GSFC-618.0)[SCIENCE SYSTEMS AND APPLICATIONS INC] via gdal-dev
That’s how I’ve always done it and so far without problems.

From: gdal-dev  on behalf of Andrew Bell via 
gdal-dev 
Date: Wednesday, May 22, 2024 at 12:36 PM
To: gdal dev 
Subject: [EXTERNAL] [BULK] [gdal-dev] Threading and Datasets
CAUTION: This email originated from outside of NASA.  Please take care when 
clicking links or opening attachments.  Use the "Report Message" button to 
report suspicious messages to the NASA SOC.


Hi,

Are there any issues with opening the same input multiple times as distinct 
datasets so that simultaneous read access can be done by threads?

Thanks,

--
Andrew Bell
andrew.bell...@gmail.com
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev


[gdal-dev] Threading and Datasets

2024-05-22 Thread Andrew Bell via gdal-dev
Hi,

Are there any issues with opening the same input multiple times as distinct
datasets so that simultaneous read access can be done by threads?

Thanks,

-- 
Andrew Bell
andrew.bell...@gmail.com
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev


Re: [gdal-dev] building just the python bindings (cmake)

2024-05-22 Thread Even Rouault via gdal-dev


Le 22/05/2024 à 16:59, Greg Troxel a écrit :

I got 3.5 to work, so I'm focusing on the 3.9 upgrade.   It seems that
3.9, vs 3.5, has withdrawn the swig generated files from the distfile.

   With 3.9, you could also do a more manual approach (not sure the
   python_generated_files target is available in 3.5)
   cmake ..
   make GDAL  # build the library
   cd swig/python
   make python_generated_files  # generate SWIG generated files
   python setup.py build

That looks like it works.  Except that I omitted the 'make GDAL'  -- in
my case, gdal is already installed.   I don't see why one would need to
build, as long as the headers and libs are available in the system.
Does that seem reasonable, or am I missing something?
The  GDAL Python (or CSharp or Java) bindings, being in-tree, have a 
special status compared to other bindings in dedicated code repository, 
since they get built as one of the many components built by the build 
process, and do not rely on installed libs/headers in the default 
process. You've figured out it is still possible to build a wheel (or 
use the source pakage on pypi https://pypi.org/project/GDAL/ which are 
generated following the procedure documented in HOWTO-RELEASE) that will 
manage to find the installed libs/headers


It would be great if README in swig/python explained this.  The
"installation" instructions are more pointers to packaging systems, and
not so much guidance for people trying to build from source.  I realize
that almost all people who use gdal use it from packaging and do not
have the slightest idea how to build perhaps anything, but on the other
hand I would expect that a README in the sources is only read by people
who are trying to build it.

Feel free to contribute documentation :-)


Also in README, it talks about numpy being a depdendency but not
formally required.   I think two clarificaions would be useful:

   What is lost specifically, vs "many examples and utilities"?  If you
   use it with qgis, without numpy, are things missing there too?  I'm
   willing to believe it's a big deal, but it's hard to understand.

   I am guessing that the things that need numpy work fine if you install
   this built in an environment without numpy and then later also install
   numpy, but it is unclear if numpy is a run-time-only dependency or if
   it is needed at build time.


It is a built-time for the osgeo.gdal_array component that is for 
example used when using the osgeo.gdal.Dataset.ReadAsArray() method.


For generic purpose packages (unless you're very very tight on the 
binary footprint), I'd strongly suggest building the Python bindings 
with the numpy dependency, and make it a run-time one as well.



--
http://www.spatialys.com
My software is free, but my time generally not.

___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev


Re: [gdal-dev] building just the python bindings (cmake)

2024-05-22 Thread Greg Troxel via gdal-dev
I got 3.5 to work, so I'm focusing on the 3.9 upgrade.   It seems that
3.9, vs 3.5, has withdrawn the swig generated files from the distfile.

  With 3.9, you could also do a more manual approach (not sure the
  python_generated_files target is available in 3.5)
  cmake ..
  make GDAL  # build the library
  cd swig/python
  make python_generated_files  # generate SWIG generated files
  python setup.py build

That looks like it works.  Except that I omitted the 'make GDAL'  -- in
my case, gdal is already installed.   I don't see why one would need to
build, as long as the headers and libs are available in the system.
Does that seem reasonable, or am I missing something?

It would be great if README in swig/python explained this.  The
"installation" instructions are more pointers to packaging systems, and
not so much guidance for people trying to build from source.  I realize
that almost all people who use gdal use it from packaging and do not
have the slightest idea how to build perhaps anything, but on the other
hand I would expect that a README in the sources is only read by people
who are trying to build it.

Also in README, it talks about numpy being a depdendency but not
formally required.   I think two clarificaions would be useful:

  What is lost specifically, vs "many examples and utilities"?  If you
  use it with qgis, without numpy, are things missing there too?  I'm
  willing to believe it's a big deal, but it's hard to understand.

  I am guessing that the things that need numpy work fine if you install
  this built in an environment without numpy and then later also install
  numpy, but it is unclear if numpy is a run-time-only dependency or if
  it is needed at build time.
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev


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

2024-05-22 Thread Rahkonen Jukka via gdal-dev
Hi,

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

-Jukka-

Lähettäjä: gdal-dev  Puolesta Paul Meems via 
gdal-dev
Lähetetty: keskiviikko 22. toukokuuta 2024 17.31
Vastaanottaja: gdal-dev@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:


01 field data-EPSG32634.csv
wkbPoint




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.7014100071e+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:


01 field data-EPSG3857.csv
wkbPoint




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.7014100071e+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.7014100071e+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


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

2024-05-22 Thread Paul Meems via gdal-dev
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:


01 field data-EPSG32634.csv
wkbPoint




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.7014100071e+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:


01 field data-EPSG3857.csv
wkbPoint




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.7014100071e+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.7014100071e+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