Re: [gdal-dev] Kriging Interpolation
Hi Evan, I probably wasn't clear in my response to you. I did try optimizing IDW using max points and various radii. The results were mixed, but overall by using radius and max points, the process is slowed considerably. If it is helpful, I can give you some time differences when I get back to testing. Usually I can work around speed issues by reducing the area I am interpolating, but there does reach a point where the area is too small for interpolating and artifacts near the edges become an issue for visualization. As far as kriging goes, for testing and using the exact same dataset, I used the default settings in Arc. So: Kriging Method: Ordinary SemiVariogram: Spherical Search Radius: Variable Number of Points: 12 For test examples, I tired to keep things as simple as possible since the more complex solutions were providing similar visual result as the basic examples. I really appreciate you looking into the processes I provided. I plan to upgrade my Gdal in few days too. You suggested in a previous thread that doing so should solve some of the speed issues. I let you know how it turned out. Thanks again, Ian From: Even Rouault [even.roua...@spatialys.com] Sent: Tuesday, 26 June 2018 12:21 a.m. To: gdal-dev@lists.osgeo.org Cc: Ian Reese; Jeremy Palmer Subject: Re: [gdal-dev] Kriging Interpolation Ian, Looking at your script, for the IDW, I see you don't use the max_points, nor the radius1 and radius2 to limit the neighbourhoud into which the points are searched. For the linear search, I've profiled the run and found that an enormous amout of time was spent processing the points outside of the triangulation. Investigating, I found there was a missed optimization that I've added per https://github.com/OSGeo/gdal/commit/134b9c4e437e347aac48c459a18a933ed526be5b With that optimization, the runtime with 8 threads is 16 sec, vs 48 sec before (on a non-optimized debug build) And if you don't want the nearest neighbour search for points outside of the triangulation, you can also add :radius=0.0, and this will speed up things again. Down to 2 sec. Regarding Kriging, I guess you specified some neighbourood value in ArcGIS ? Otherwise as mentionned by Thomas, Knudsen, the computation time and memory requirements would be enormous (between O(n^2) and O(n^3) where n is the number of points). I gave a try to https://github.com/bsmurphy/PyKrige and used Ordinary Kriging. For the above reasons, it blows up with a MemoryError on the full layer. I also found that the extent of the output grid had also a strong influence on the memory requirements. Which is less expected, perhaps a side effect of that particular implementation. Anyway, if arbitrarily restricting to the first 1000 points of your layer (which is roughly a strip in the middle of the dataset) and to their extent (producing a 1321x45 raster with 1m resolution), it gives a result (2 minutes single-threaded, 1.8 GB RAM), which looks consistent with the one from ArcGIS. This is not exactly the same due to a non-moving neighbouroud having been used due to my simplified methodology, but for the 'inner' of the computed raster, this is quite similar. So for Kriging, a moving window approach limiting the number of input points would be absolutely needed for "large" point datasets, and the processing time would still be significantly slower than linear interpolation. Even > Hi Gdal, > > I'd like to know if it is possible to make Kriging Interpolation part of > gdal_grid. We work with a good deal of scattered point elevation data and > neither IDW or Linear interpolations(smooth or not smooth) produce the > desired visual results. Nor are they fast enough for the density of data > we commonly work with. Linear interpolation is a reliable and fast method, > however the triangulated results are difficult to work with visually. > > For some visual comparisons and reference, here is a link to a folder of > sample points(shp) for a test area, two examples of kriged data for that > test region, and a BASH script to reproduce IDW and Linear interpolations > of the same region. > > https://drive.google.com/open?id=1HKbK4wMmVN6JZSvhxBQDub1WXBLE710H > > Here is an example of the Kriging method used by ArcGIS: > http://desktop.arcgis.com/en/arcmap/latest/extensions/geostatistical-analyst > /kriging-in-geostatistical-analyst.htm > > Few graphs explaining best interpolation methods to use with given data. > Gridded vs. scattered. > https://www.neonscience.org/spatial-interpolation-basics > > If it is possible, we are curious to know what type of funding would be > needed. > > Cheers, > > Ian > > > > This message contains information, which may be in confidence and may be > subject to legal privilege. If you are not the intended recipient, you must > not peruse, use, disseminate, distribute or copy this message. If you have > received this message in error, please notify us
Re: [gdal-dev] Kriging Interpolation
Cheer Joakim! That I am not familiar with and will have a look From: Joaquim Manuel Freire Luís [jl...@ualg.pt] Sent: Tuesday, 26 June 2018 10:58 a.m. To: Ian Reese; Even Rouault; gdal-dev@lists.osgeo.org Cc: Jeremy Palmer Subject: RE: [gdal-dev] Kriging Interpolation I have plenty of methods for working around speed issues and have |>worked with any number of interpolations methods in GRASS, SAGA and |>Python. You have also the minimum curvature in GMT. Not lightning fast but quicker than krigging and normally better looking, specially on the extrapolation zones. Joaquim This message contains information, which may be in confidence and may be subject to legal privilege. If you are not the intended recipient, you must not peruse, use, disseminate, distribute or copy this message. If you have received this message in error, please notify us immediately (Phone 0800 665 463 or i...@linz.govt.nz) and destroy the original message. LINZ accepts no responsibility for changes to this email, or for any attachments, after its transmission from LINZ. Thank You. ___ gdal-dev mailing list gdal-dev@lists.osgeo.org https://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] Kriging Interpolation
I have plenty of methods for working around speed issues and have |>worked with any number of interpolations methods in GRASS, SAGA and |>Python. You have also the minimum curvature in GMT. Not lightning fast but quicker than krigging and normally better looking, specially on the extrapolation zones. Joaquim ___ gdal-dev mailing list gdal-dev@lists.osgeo.org https://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] Esri JSON Curves
> This would still be an issue with paging. Most of the actual datasets > of interest (to me) involve multiple HTTP requests to fetch many pages > of JSON data (usually 1000 - 1 features at a time depending on > server setup). The curves may easily not be on the first page and > re-running the whole request is expensive. Fair enough > > If my memories are right, you only need to set it for write support. So > > shouldn't be needed there. > > All I know is that without it set, the curved geometries were linearized > somewhere in OGR before they were passed to the PostGIS driver. You should double check and spot where this comes from. You can set a breakpoint on OGRGeometryFactory::curveToLineString which should be taken when linearization occurs. As an experience, I've disabled OLCCurveGeometries and ODsCCurveGeometries capability declaration in the CSV driver, and no linearization occurs when converting to PostGIS Even -- Spatialys - Geospatial professional services http://www.spatialys.com ___ gdal-dev mailing list gdal-dev@lists.osgeo.org https://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] Esri JSON Curves
On 06/22/2018 02:51 PM, Even Rouault wrote: > On vendredi 22 juin 2018 14:38:15 CEST James Klassen wrote: >> I finally got around to doing a rough implementation of parsing ESRIJSON >> with curved geometries. The code is very rough yet, but works enough to >> import a polygon layer from an ESRI FeatureService into PostGIS with the >> curves intact. Code is available at [1]. >> >> Main issues: >> >> I don't know how to detect if there will be curves in the source file >> without first parsing the whole file, so currently I just assume a curved >> geometry type. This might be a candidate for an open option. > The ESRI Json reader ingests everything in memory and builds the JSon tree > structure, so doing an initial pass by default shouldn't be a big deal. This would still be an issue with paging. Most of the actual datasets of interest (to me) involve multiple HTTP requests to fetch many pages of JSON data (usually 1000 - 1 features at a time depending on server setup). The curves may easily not be on the first page and re-running the whole request is expensive. >> My focus has been on polygon layers (because that's the dataset I need to >> read). I doubt it works outside of polygon layers at the moment. >> >> I had to enable the "Curve Geometries" capability in the top level GeoJSON >> driver (which may break some things for the GeoJSON/TopoJSON drivers). > If my memories are right, you only need to set it for write support. So > shouldn't be needed there. > All I know is that without it set, the curved geometries were linearized somewhere in OGR before they were passed to the PostGIS driver. ___ gdal-dev mailing list gdal-dev@lists.osgeo.org https://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] Motion: Promote GDAL 2.3.1 rc2 for release
+1 Tamas 2018-06-25 15:18 GMT+02:00 Even Rouault : > Hi, > > Motion: GDAL/OGR 2.3.1 rc2 is promoted to be the official 2.3.1 final > release. > > --- > > My vote: +1 > > Even > > -- > Spatialys - Geospatial professional services > http://www.spatialys.com > ___ > gdal-dev mailing list > gdal-dev@lists.osgeo.org > https://lists.osgeo.org/mailman/listinfo/gdal-dev ___ gdal-dev mailing list gdal-dev@lists.osgeo.org https://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] [EXTERNAL] Convert VRT+JPEG with gdalwarp
On lundi 25 juin 2018 13:10:16 CEST Kalev Julge wrote:
> Hi Doug,
>
> Thanks for the info. I know that JPEG is lossy but our datasets are very
> large. That’s why we use JPEG, which is good enough for our needs and the
> file sizes are smaller.
After your gdalwarp to VRT, do:
gdal_translate new.vrt new.jpg -of JPEG
You cannot directly wrap to JPEG, because of the driver limitation
See
https://trac.osgeo.org/gdal/wiki/FAQRaster#Whywontgdalwarporgdal_mergewritetomostformats
> It seems that GDAL doesn’t create new JPEG files without having pre-existing
> JPEG-s to modify.
Actually, that's right the opposite :-)
> If I try to modify old JPEG, I get a message “Output
> dataset exists, but cannot be opened in update mode.
That means you cannot modify an existing JPEG. You have to recreate a
completely new file.
> I also tried using “CreateCopy” but that produced syntax error. Used command
> on GDAL CMD:
“GDALDriver::CreateCopy ("new.jpg","old.jpg", FALSE, NULL,
> NULL, NULL)”
The second argument must be a GDALDataset* object, not a string, so do
something like
GDALDataset* srcDS = reinterpret_cast(
GDALOpen("old.jpg", GA_ReadOnly));
GDALDataset* outDS = jpegDriver->CreateCopy(
"new.jpg", srcDS, FALSE, NULL, NULL, NULL);
GDALClose(outDS);
GDALClose(srcDS);
Even
--
Spatialys - Geospatial professional services
http://www.spatialys.com
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev
[gdal-dev] Motion: Promote GDAL 2.3.1 rc2 for release
Hi, Motion: GDAL/OGR 2.3.1 rc2 is promoted to be the official 2.3.1 final release. --- My vote: +1 Even -- Spatialys - Geospatial professional services http://www.spatialys.com ___ gdal-dev mailing list gdal-dev@lists.osgeo.org https://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] [EXTERNAL] Convert VRT+JPEG with gdalwarp
Hi Doug,
Thanks for the info. I know that JPEG is lossy but our datasets are very large.
That’s why we use JPEG, which is good enough for our needs and the file sizes
are smaller.
It seems that GDAL doesn’t create new JPEG files without having pre-existing
JPEG-s to modify. If I try to modify old JPEG, I get a message “Output dataset
exists, but cannot be opened in update mode.
I also tried using “CreateCopy” but that produced syntax error. Used command on
GDAL CMD:
“GDALDriver::CreateCopy ("new.jpg","old.jpg", FALSE, NULL, NULL, NULL)”
Kalev
From: Doug Newcomb
Sent: Thursday, June 21, 2018 5:24 PM
To: Kalev Julge
Cc: gdal dev
Subject: Re: [EXTERNAL] [gdal-dev] Convert VRT+JPEG with gdalwarp
Kalev,
Just set the output file format to -f JPEG format ( .jpg) or -f JP2OpenJPEG
(.jp2)
The JPEG format will be lossey by default. Is that what you want?
See:
https://medium.com/planet-stories/a-gentle-introduction-to-gdal-part-1-a3253eb96082
https://medium.com/planet-stories/a-gentle-introduction-to-gdal-part-2-map-projections-gdalwarp-e05173bd710a
Doug
On Thu, Jun 21, 2018 at 10:00 AM kalevj
mailto:kalev.ju...@reach-u.com>> wrote:
Hello,
I used gdalwarp to transform a VRT+JPEG pair from one EPSG projection to
another. Everything worked well. This conversion changed the VRT file so
that reprojection takes place when the images are loaded. However, I would
need to reproject the JPEG so that there would not be a need to reproject
every time the images are loaded. Is it possible somehow?
I used the command:
"gdalwarp -s_srs EPSG:2006 -t_srs EPSG:3857 -of VRT -r near old.vrt new.vrt"
Thanks!
Kalev
--
Sent from: http://osgeo-org.1560.x6.nabble.com/GDAL-Dev-f3742093.html
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev
--
Doug Newcomb
USFWS
551F Pylon Dr
Raleigh, NC
919-856-4520 ext. 14 doug_newc...@fws.gov
-
NOTE: This email correspondence and any attachments to and from this sender is
subject to the Freedom of Information Act (FOIA) and may be disclosed to third
parties.
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] Kriging Interpolation
Ian, Looking at your script, for the IDW, I see you don't use the max_points, nor the radius1 and radius2 to limit the neighbourhoud into which the points are searched. For the linear search, I've profiled the run and found that an enormous amout of time was spent processing the points outside of the triangulation. Investigating, I found there was a missed optimization that I've added per https://github.com/OSGeo/gdal/commit/134b9c4e437e347aac48c459a18a933ed526be5b With that optimization, the runtime with 8 threads is 16 sec, vs 48 sec before (on a non-optimized debug build) And if you don't want the nearest neighbour search for points outside of the triangulation, you can also add :radius=0.0, and this will speed up things again. Down to 2 sec. Regarding Kriging, I guess you specified some neighbourood value in ArcGIS ? Otherwise as mentionned by Thomas, Knudsen, the computation time and memory requirements would be enormous (between O(n^2) and O(n^3) where n is the number of points). I gave a try to https://github.com/bsmurphy/PyKrige and used Ordinary Kriging. For the above reasons, it blows up with a MemoryError on the full layer. I also found that the extent of the output grid had also a strong influence on the memory requirements. Which is less expected, perhaps a side effect of that particular implementation. Anyway, if arbitrarily restricting to the first 1000 points of your layer (which is roughly a strip in the middle of the dataset) and to their extent (producing a 1321x45 raster with 1m resolution), it gives a result (2 minutes single-threaded, 1.8 GB RAM), which looks consistent with the one from ArcGIS. This is not exactly the same due to a non-moving neighbouroud having been used due to my simplified methodology, but for the 'inner' of the computed raster, this is quite similar. So for Kriging, a moving window approach limiting the number of input points would be absolutely needed for "large" point datasets, and the processing time would still be significantly slower than linear interpolation. Even > Hi Gdal, > > I'd like to know if it is possible to make Kriging Interpolation part of > gdal_grid. We work with a good deal of scattered point elevation data and > neither IDW or Linear interpolations(smooth or not smooth) produce the > desired visual results. Nor are they fast enough for the density of data > we commonly work with. Linear interpolation is a reliable and fast method, > however the triangulated results are difficult to work with visually. > > For some visual comparisons and reference, here is a link to a folder of > sample points(shp) for a test area, two examples of kriged data for that > test region, and a BASH script to reproduce IDW and Linear interpolations > of the same region. > > https://drive.google.com/open?id=1HKbK4wMmVN6JZSvhxBQDub1WXBLE710H > > Here is an example of the Kriging method used by ArcGIS: > http://desktop.arcgis.com/en/arcmap/latest/extensions/geostatistical-analyst > /kriging-in-geostatistical-analyst.htm > > Few graphs explaining best interpolation methods to use with given data. > Gridded vs. scattered. > https://www.neonscience.org/spatial-interpolation-basics > > If it is possible, we are curious to know what type of funding would be > needed. > > Cheers, > > Ian > > > > This message contains information, which may be in confidence and may be > subject to legal privilege. If you are not the intended recipient, you must > not peruse, use, disseminate, distribute or copy this message. If you have > received this message in error, please notify us immediately (Phone 0800 > 665 463 or i...@linz.govt.nz) and destroy the original message. LINZ > accepts no responsibility for changes to this email, or for any > attachments, after its transmission from LINZ. Thank You. -- Spatialys - Geospatial professional services http://www.spatialys.com ___ gdal-dev mailing list gdal-dev@lists.osgeo.org https://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] [gdal-announce] GDAL Support For PostGIS Out-DB Raster
Hi Even, Thanks a lot. This is great news! I will start investigating the fix for Out-DB rasters. I will keep you posted. Great work! On Mon, Jun 25, 2018 at 10:26 AM, Even Rouault wrote: > On lundi 25 juin 2018 10:18:22 CEST Osahon Oduware wrote: > > Hi Even, > > > > I want to find out if this GDAL fix will be release by the end of June. I > > use PostGIS version 2.3.2 at present and I will be building from source. > > Yes, this is in the release notes of the release candidate of GDAL 2.3.1 > that > was published last Friday: > https://trac.osgeo.org/gdal/wiki/Release/2.3.1-News > > -- > Spatialys - Geospatial professional services > http://www.spatialys.com > ___ gdal-dev mailing list gdal-dev@lists.osgeo.org https://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] [gdal-announce] GDAL Support For PostGIS Out-DB Raster
On lundi 25 juin 2018 10:18:22 CEST Osahon Oduware wrote: > Hi Even, > > I want to find out if this GDAL fix will be release by the end of June. I > use PostGIS version 2.3.2 at present and I will be building from source. Yes, this is in the release notes of the release candidate of GDAL 2.3.1 that was published last Friday: https://trac.osgeo.org/gdal/wiki/Release/2.3.1-News -- Spatialys - Geospatial professional services http://www.spatialys.com ___ gdal-dev mailing list gdal-dev@lists.osgeo.org https://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] [gdal-announce] GDAL Support For PostGIS Out-DB Raster
Hi Even, I want to find out if this GDAL fix will be release by the end of June. I use PostGIS version 2.3.2 at present and I will be building from source. On Thu, May 31, 2018 at 9:39 PM, Even Rouault wrote: > On jeudi 31 mai 2018 21:35:30 CEST Osahon Oduware wrote: > > Hi Even, > > > > Thanks for the prompt response. Please, I want to know what is the > minimum > > version of GDAL, > no released version yet. Will be in 2.3.1 (end of June probably) and 2.4.0 > (next year) > > PostgreSQL > irrelevant. Any PostgreSQL version supported by the PostGIS version you > use is > fine > > and Postgis > will be in Postgis 2.5.0. But not strictly needed. Existing versions will > work > too > > > that this outdb raster > > functionality is currently implemented. > > > > Best Regards > > > > On May 31, 2018 20:13, "Even Rouault" > wrote: > > > On jeudi 31 mai 2018 20:04:58 CEST Osahon Oduware wrote: > > > > Hi Even, > > > > > > > > Any update on this? I have been away for a while. > > > > > > Yep, support for out-db raster has been added in GDAL master, and > > > backported > > > to 2.3 branch as well > > > > > > See outdb_resolution option in http://gdal.org/frmt_postgisraster.html > > > > > > 2 new functions added to PostGIS 2.5dev as well to help reading the > out-db > > > raster directly from client-side if client and server see the same > raster. > > > > > > -- > > > Spatialys - Geospatial professional services > > > http://www.spatialys.com > > > -- > Spatialys - Geospatial professional services > http://www.spatialys.com > ___ gdal-dev mailing list gdal-dev@lists.osgeo.org https://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] Kriging Interpolation
Hi Ian, If you find IDW on N points slow, you will find Kriging on the same number of points unbearably slow, since for Kriging you will have to solve a system of N linear equations (essentially inverting an NxN matrix) to compute the Kriging weights. To set up the matrix you will even need to compute the same distance elements you need for IDW. So Kriging will never become faster than IDW. Hence, the only realistic way to use Kriging on large data sets is to base each estimate on local subsets - typically by working in a TIN model, just like in the linear interpolation case. And if your data set is dense, the difference between Kriging and linear interpolation will be miniscule. Kriging will, however, provide a variance estimate for the prediction, which may be realistic assuming you provide a realistic covariance function for the phenomenon under investigation. So to make any difference wrt the linear interpolation, you will need to include a larger number of data, e.g. by including neighbouring triangles of the TIN, or using another data structure to subset your observations. Some years ago I wrote an informal note about this subject. The title is “Data Assimilation for Updates of Digital Terrain Models” and it is available over at http://sdfe.dk/media/2916626/kms_technical_report_14.pdf /Thomas 2018-06-25 5:19 GMT+02:00 Ian Reese : > Hi Gdal, > > I'd like to know if it is possible to make Kriging Interpolation part of > gdal_grid. We work with a good deal of scattered point elevation data and > neither IDW or Linear interpolations(smooth or not smooth) produce the > desired visual results. Nor are they fast enough for the density of data > we commonly work with. Linear interpolation is a reliable and fast method, > however the triangulated results are difficult to work with visually. > > For some visual comparisons and reference, here is a link to a folder of > sample points(shp) for a test area, two examples of kriged data for that > test region, and a BASH script to reproduce IDW and Linear interpolations > of the same region. > > https://drive.google.com/open?id=1HKbK4wMmVN6JZSvhxBQDub1WXBLE710H > > Here is an example of the Kriging method used by ArcGIS: > http://desktop.arcgis.com/en/arcmap/latest/extensions/ > geostatistical-analyst/kriging-in-geostatistical-analyst.htm > > Few graphs explaining best interpolation methods to use with given data. > Gridded vs. scattered. > https://www.neonscience.org/spatial-interpolation-basics > > If it is possible, we are curious to know what type of funding would be > needed. > > Cheers, > > Ian > > -- > > This message contains information, which may be in confidence and may be > subject to legal privilege. If you are not the intended recipient, you must > not peruse, use, disseminate, distribute or copy this message. If you have > received this message in error, please notify us immediately (Phone 0800 > 665 463 or i...@linz.govt.nz) and destroy the original message. LINZ > accepts no responsibility for changes to this email, or for any > attachments, after its transmission from LINZ. Thank You. > > ___ > gdal-dev mailing list > gdal-dev@lists.osgeo.org > https://lists.osgeo.org/mailman/listinfo/gdal-dev > ___ gdal-dev mailing list gdal-dev@lists.osgeo.org https://lists.osgeo.org/mailman/listinfo/gdal-dev
