[GRASS-dev] Re: [GRASS GIS] #1088: r.fillnulls: support other interpolation methods
#1088: r.fillnulls: support other interpolation methods
-+--
Reporter: kyngchaos| Owner: grass-dev@…
Type: enhancement | Status: new
Priority: normal | Milestone: 6.4.3
Component: Raster | Version: svn-develbranch6
Keywords: r.fillnulls |Platform: All
Cpu: All |
-+--
Comment(by mmetz):
Replying to [comment:26 glynn]:
Replying to [comment:24 mmetz]:
Ok, r.buffer now (r50830) uses squared distance for metric=squared and
regular distance in meters for metric=geodesic.
I think that r.grow.distance needs this, right?
{{{
--- raster/r.grow.distance/main.c (revision 50836)
+++ raster/r.grow.distance/main.c (working copy)
@@ -208,8 +208,11 @@
else
G_fatal_error(_(Unknown metric: '%s'), opt.met-answer);
-if (flag.m-answer)
+if (flag.m-answer) {
scale = G_database_units_to_meters_factor();
+ if (strcmp(opt.met-answer, squared) == 0)
+ scale *= scale;
+}
in_fd = Rast_open_old(in_name, );
}}}
Right. It worked in my test because map units were meters and scale thus
1. Patch applied in r50841.
Markus M
--
Ticket URL: https://trac.osgeo.org/grass/ticket/1088#comment:28
GRASS GIS http://grass.osgeo.org
___
grass-dev mailing list
grass-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/grass-dev
[GRASS-dev] Re: [GRASS GIS] #1088: r.fillnulls: support other interpolation methods
#1088: r.fillnulls: support other interpolation methods -+-- Reporter: kyngchaos| Owner: grass-dev@… Type: enhancement | Status: new Priority: normal | Milestone: 6.4.3 Component: Raster | Version: svn-develbranch6 Keywords: r.fillnulls |Platform: All Cpu: All | -+-- Comment(by mmetz): Replying to [comment:23 glynn]: Replying to [comment:14 mmetz]: PS: I fixed r.buffer in trunk, it could now be used (again) by r.fillnulls instead of r.grow. What was the problem? metric=geodesic does not return squared distances, and for other metrics the distances as returned by r.grow.distance were in map units, not meters, but r.buffer requires meter as unit. Because the fix is needlessly inefficient. The efficient way to compare Euclidean distance against a fixed value is distance!^2 limit!^2 rather than sqrt(distance!^2) limit. Calculating Euclidean distance invariably involves calculating distance-squared as an intermediate value; if you only need the value for a comparison, squaring one value is faster than taking the square root of the other. Ok, r.buffer now (r50830) uses squared distance for metric=squared and regular distance in meters for metric=geodesic. Markus M -- Ticket URL: https://trac.osgeo.org/grass/ticket/1088#comment:24 GRASS GIS http://grass.osgeo.org ___ grass-dev mailing list grass-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/grass-dev
[GRASS-dev] Re: [GRASS GIS] #1088: r.fillnulls: support other interpolation methods
#1088: r.fillnulls: support other interpolation methods -+-- Reporter: kyngchaos| Owner: grass-dev@… Type: enhancement | Status: new Priority: normal | Milestone: 6.4.3 Component: Raster | Version: svn-develbranch6 Keywords: r.fillnulls |Platform: All Cpu: All | -+-- Comment(by mmetz): Replying to [comment:22 hamish]: Replying to [comment:21 mmetz]: r.fillnulls should do just that, filling nulls, and not do resampling. I am not convinced (not even a little). No exceptions to the standard map creation conventions unless it is absolutely fundamental to what the module does. If this was a edit-in-place module like r.null maybe I could buy the argument, but it isn't. Either deal with all of the input map's native bounds, or use the current region, not half and half. Using the current region settings gives the user the choice of which one they want. In theory I would agree, but not everybody knows that for a 3 arcsec SRTM g.region -a res=00:00:03 is wrong and only g.region align=SRTM properly aligns the region to the raster map. If you want to respect the current region, you would need to remove the align=input option from g.region in r.fullnulls when you zoom to the holes. You could zoom to the raster buffer instead which has been created using the current region. Markus M -- Ticket URL: https://trac.osgeo.org/grass/ticket/1088#comment:25 GRASS GIS http://grass.osgeo.org ___ grass-dev mailing list grass-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/grass-dev
[GRASS-dev] Re: [GRASS GIS] #1088: r.fillnulls: support other interpolation methods
#1088: r.fillnulls: support other interpolation methods
-+--
Reporter: kyngchaos| Owner: grass-dev@…
Type: enhancement | Status: new
Priority: normal | Milestone: 6.4.3
Component: Raster | Version: svn-develbranch6
Keywords: r.fillnulls |Platform: All
Cpu: All |
-+--
Comment(by glynn):
Replying to [comment:24 mmetz]:
Ok, r.buffer now (r50830) uses squared distance for metric=squared and
regular distance in meters for metric=geodesic.
I think that r.grow.distance needs this, right?
{{{
--- raster/r.grow.distance/main.c (revision 50836)
+++ raster/r.grow.distance/main.c (working copy)
@@ -208,8 +208,11 @@
else
G_fatal_error(_(Unknown metric: '%s'), opt.met-answer);
-if (flag.m-answer)
+if (flag.m-answer) {
scale = G_database_units_to_meters_factor();
+ if (strcmp(opt.met-answer, squared) == 0)
+ scale *= scale;
+}
in_fd = Rast_open_old(in_name, );
}}}
--
Ticket URL: https://trac.osgeo.org/grass/ticket/1088#comment:26
GRASS GIS http://grass.osgeo.org
___
grass-dev mailing list
grass-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/grass-dev
[GRASS-dev] Re: [GRASS GIS] #1088: r.fillnulls: support other interpolation methods
#1088: r.fillnulls: support other interpolation methods -+-- Reporter: kyngchaos| Owner: grass-dev@… Type: enhancement | Status: new Priority: normal | Milestone: 6.4.3 Component: Raster | Version: svn-develbranch6 Keywords: r.fillnulls |Platform: All Cpu: All | -+-- Comment(by hamish): none of this gets away from the fact that r.buffer.py is consistently ~ 40x slower than r.buffer2 at all(?) region sizes. For my 8GB workstation the extra RAM needed by the classic r.buffer.c would mean I'd have to transition to the less ram-hungry r.grow.distance method at a region size of 9x9 (mem. req. here is 1 byte per region cell), which is, so far, still a corner case. r.buffer.py is an interesting exercise, but I'm not convinced that the r.grow.distance method needs to be anything more than a full example in the r.buffer help page, wrt alternate steps to take if it wants to use more RAM than you've got. regards, Hamish -- Ticket URL: https://trac.osgeo.org/grass/ticket/1088#comment:27 GRASS GIS http://grass.osgeo.org ___ grass-dev mailing list grass-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/grass-dev
[GRASS-dev] Re: [GRASS GIS] #1088: r.fillnulls: support other interpolation methods
#1088: r.fillnulls: support other interpolation methods -+-- Reporter: kyngchaos| Owner: grass-dev@… Type: enhancement | Status: new Priority: normal | Milestone: 6.4.3 Component: Raster | Version: svn-develbranch6 Keywords: r.fillnulls |Platform: All Cpu: All | -+-- Comment(by hamish): Replying to [comment:20 hamish]: hmmm, what if r.buffer.py checks rows*columns and decides which backend it should call based on the result. small regions get sent to one, large ones to another, best of both worlds. Replying to [comment:21 mmetz]: Small and large is relative to the amount of free memory. You would need to check for that first before you decide what is small and what is large. IMHO too complicated, I would just use the safer method. see attached patch for a rough idea. check free memory is not a road I want to go down, so set a default at 12000x12000 cells (I just picked a number), and if it's bigger than that use the r.grow.distance method, otherwise the fast but memory intensive method. Add an option to allow the user to override the default. 4.5 sec for .py vs. 0.140 sec for .c. I know I've got at least one script which runs r.buffer in a loop for thousands of starting points, I don't want to think how long it will take at 4.5sec each iteration. (actually 0.250 sec for r.buffer2 if you use the python wrapper) for a 12000x12000 region using r.buffer2: 17.4 sec, 156mb RAM. same using r.grow.distance: 10 minutes 7 second, 20mb RAM (takes forever reading and writing map, even on a fast RAID [it's CPU bound not IO bound]) 2x2 region: r.buffer2 takes 51 sec, uses 400mb RAM (1 byte per cell). I'm not going to bother to find out how long r.grow.distance would take. 4x4 region: r.buffer2 takes 4 min 7sec (~1 min to read, ~1 min to buffer, ~2 min to write), uses 1.6gb RAM. r.grow.distance, no idea, maybe days to complete. 95% of work will be smaller than 40k^2 cells. the r.grow.distance method seems like it is best left as a footnote + example in the r.buffer2 help page. r.fillnulls should do just that, filling nulls, and not do resampling. I am not convinced (not even a little). No exceptions to the standard map creation conventions unless it is absolutely fundamental to what the module does. If this was a edit-in-place module like r.null maybe I could buy the argument, but it isn't. Either deal with all of the input map's native bounds, or use the current region, not half and half. Using the current region settings gives the user the choice of which one they want. thanks, Hamish -- Ticket URL: https://trac.osgeo.org/grass/ticket/1088#comment:22 GRASS GIS http://grass.osgeo.org ___ grass-dev mailing list grass-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/grass-dev
[GRASS-dev] Re: [GRASS GIS] #1088: r.fillnulls: support other interpolation methods
#1088: r.fillnulls: support other interpolation methods -+-- Reporter: kyngchaos| Owner: grass-dev@… Type: enhancement | Status: new Priority: normal | Milestone: 6.4.3 Component: Raster | Version: svn-develbranch6 Keywords: r.fillnulls |Platform: All Cpu: All | -+-- Comment(by mmetz): Replying to [comment:17 hamish]: Replying to [comment:15 hamish]: any idea about the lambda_i setting? Replying to [comment:16 mmetz]: In my tests, 1 is usually way too high, 0.005 seems to do well in mountainous areas, 0.1 is better for nearly flat areas. in trunk and devbr6 I have changed it to 0.01 now. I would suggest to backport to 6.4svn ASAP before it harms any more data, but would like to test it with some meter-based LIDAR data first, as I've only tried with lat/lon srtm so far. Tested with LiDAR datasets, lambda=1 is too high, synced in 6.4 to trunk/devbr6 with default lambda=0.01. I think it is not easy to avoid segmentation artifacts with RST (very high npmin, very low segmax maybe). for r.fillnulls I'd expect the 3-cell buffer of most holes not to fall on segmentation boundaries. it wasn't that it had obvious segmentation lines in it, rather it just didn't follow the ridges and valleys in such a natural way. (qualitative eyeball analysis) I've been convinced for years that there must be a way to avoid the segmentation artifacts without massively overlapping the quadtree boxes.. still keeping an eye out. I had obvious segmentation lines with RST. From my testing with the bspline method I think that there is no way around massively overlapping quadtree boxes. Hmm. One of the two r.buffer's should probably be deactivated, after comparative testing with latlong and real projections. I get 4.3sec for a test in r.buffer.py, 0.140 sec for r.buffer.c, for a simple test starting at a single cell, with visually identical results. Actually I was pleasantly surprised that r.buffer.py did the right thing in lat/lon and made a pear-shaped buffer for the plate carree display, wider nearer to the poles. In trunk, r.buffer.c loads both input and output completely to memory, whereas r.buffer.py using r.grow.distance doesn't. Therefore r.buffer.c is faster for smaller maps but will freeze the system if it goes into swap space. r.buffer.py keeps the disk busy, but the system remains responsive. Markus M -- Ticket URL: https://trac.osgeo.org/grass/ticket/1088#comment:18 GRASS GIS http://grass.osgeo.org ___ grass-dev mailing list grass-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/grass-dev
[GRASS-dev] Re: [GRASS GIS] #1088: r.fillnulls: support other interpolation methods
#1088: r.fillnulls: support other interpolation methods -+-- Reporter: kyngchaos| Owner: grass-dev@… Type: enhancement | Status: new Priority: normal | Milestone: 6.4.3 Component: Raster | Version: svn-develbranch6 Keywords: r.fillnulls |Platform: All Cpu: All | -+-- Comment(by cmbarton): The trick on sin and sie (spline interval north and spline interval east) is to set them to the minimum NS (sin) and EW (sie) distance between points. A first run of v.surf.bspline -e will provide a decent estimate of this to use. If you set these too large, it will take forever to run. Michael -- Ticket URL: https://trac.osgeo.org/grass/ticket/1088#comment:19 GRASS GIS http://grass.osgeo.org ___ grass-dev mailing list grass-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/grass-dev
[GRASS-dev] Re: [GRASS GIS] #1088: r.fillnulls: support other interpolation methods
#1088: r.fillnulls: support other interpolation methods -+-- Reporter: kyngchaos| Owner: grass-dev@… Type: enhancement | Status: new Priority: normal | Milestone: 6.4.3 Component: Raster | Version: svn-develbranch6 Keywords: r.fillnulls |Platform: All Cpu: All | -+-- Comment(by hamish): Replying to mmetz: I have modified r.fillnulls in trunk to use r.resamp.bspline which is designed for raster interpolation instead v.surf.bspline. my test 2x3 deg srtm mosaic takes 40 seconds to complete with 'r.resamp.bspline method=bilinear' as the backend and 2 minutes using bicubic (slightly better results). v.surf.bspline was taking ~4-6 minutes for bicubic, and r.surf.rst takes about 4 minutes with the zoom-to-data trick. both rst(no zoom) and v.surf.bspline's bicubic took about 8 minutes in grass 6 (no OpenMP for both there [still inefficiently used, but helps a bit]). Replying to [comment:18 mmetz]: In trunk, r.buffer.c loads both input and output completely to memory, whereas r.buffer.py using r.grow.distance doesn't. Therefore r.buffer.c is faster for smaller maps but will freeze (ie degrade) the system if it goes into swap space. r.buffer.py keeps the disk busy, but the system remains responsive. hmmm, what if r.buffer.py checks rows*columns and decides which backend it should call based on the result. small regions get sent to one, large ones to another, best of both worlds. (as long as the two r.buffer2/r.grow.distance backends give identical results, which visually they do) btw, you'll need to move that del_temp_region back to where it was. with it still in place when you run r.patch you're cropping the patched result with the zoomed-in region. Oops. Changed in r50803. The grid geometry of the input raster is still maintained, though. yeah, but if the interp and the orig are patched together on the orig's native grid (I'm not sure if it makes a practical difference here, but no objection to that), it will still need a 'r.mapcalc resamp = patched step to follow the rule that output maps are created using the (real) current region's grid settings, not the input map's native alignment. the alternative is to decalre that r.fillnulls works on the input map's native region completely and ignores the current region settings completely. but that stops you from working on sub-windows of a larger dataset. what would it take to get 'g.region vect=' to work with a level 1 map? The routines are already in v.info (trunk). shall we move vector/v.info/level1.c into libvector as Vect_level1_describe(Plus_head *) or similar? seems quite reusable. Hamish -- Ticket URL: https://trac.osgeo.org/grass/ticket/1088#comment:20 GRASS GIS http://grass.osgeo.org ___ grass-dev mailing list grass-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/grass-dev
[GRASS-dev] Re: [GRASS GIS] #1088: r.fillnulls: support other interpolation methods
#1088: r.fillnulls: support other interpolation methods -+-- Reporter: kyngchaos| Owner: grass-dev@… Type: enhancement | Status: new Priority: normal | Milestone: 6.4.3 Component: Raster | Version: svn-develbranch6 Keywords: r.fillnulls |Platform: All Cpu: All | -+-- Comment(by mmetz): Replying to [comment:20 hamish]: Replying to [comment:18 mmetz]: In trunk, r.buffer.c loads both input and output completely to memory, whereas r.buffer.py using r.grow.distance doesn't. Therefore r.buffer.c is faster for smaller maps but will freeze (ie degrade) the system if it goes into swap space. r.buffer.py keeps the disk busy, but the system remains responsive. hmmm, what if r.buffer.py checks rows*columns and decides which backend it should call based on the result. small regions get sent to one, large ones to another, best of both worlds. Small and large is relative to the amount of free memory. You would need to check for that first before you decide what is small and what is large. IMHO too complicated, I would just use the safer method. btw, you'll need to move that del_temp_region back to where it was. with it still in place when you run r.patch you're cropping the patched result with the zoomed-in region. Oops. Changed in r50803. The grid geometry of the input raster is still maintained, though. yeah, but if the interp and the orig are patched together on the orig's native grid (I'm not sure if it makes a practical difference here, but no objection to that), it will still need a 'r.mapcalc resamp = patched step to follow the rule that output maps are created using the (real) current region's grid settings, not the input map's native alignment. the alternative is to decalre that r.fillnulls works on the input map's native region completely and ignores the current region settings completely. but that stops you from working on sub-windows of a larger dataset. Working with subwindows is supported, the subwindow is aligned to the input, though. Surface resampling is another issue and I would rather avoid implicit nearest-neighbour resampling. The standard resampling modules for surfaces are r.resamp.interp (anything but nearest) and r.resamp.stats. r.fillnulls should do just that, filling nulls, and not do resampling. Working with a subregion is fine of course. what would it take to get 'g.region vect=' to work with a level 1 map? The routines are already in v.info (trunk). shall we move vector/v.info/level1.c into libvector as Vect_level1_describe(Plus_head *) or similar? seems quite reusable. Could make sense, I would like to think a bit about it. Markus M -- Ticket URL: https://trac.osgeo.org/grass/ticket/1088#comment:21 GRASS GIS http://grass.osgeo.org ___ grass-dev mailing list grass-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/grass-dev
[GRASS-dev] Re: [GRASS GIS] #1088: r.fillnulls: support other interpolation methods
#1088: r.fillnulls: support other interpolation methods -+-- Reporter: kyngchaos| Owner: grass-dev@… Type: enhancement | Status: new Priority: normal | Milestone: 6.4.3 Component: Raster | Version: svn-develbranch6 Keywords: r.fillnulls |Platform: All Cpu: All | -+-- Comment(by mmetz): Replying to [comment:12 hamish]: * applied in devbr6 by MarkusN in r50114, with sie,sin = 3*res and the default lambda_i=1. * applied to trunk by luca in r50128 with sie,sin = 1*res and the default lambda_i=1. see above comments for MarkusM's recommendation for using si=2*res and lambda=0.005. The sie/sin settings are the result of guestimation and experimenting. I found that for evenly placed points, e.g. a vectorized raster without gaps, 1 * res gives nice results. With r.fillnulls, the points are not evenly spaced and sie/sin needs to be larger than 1*res in order to avoid artifacts. For very large gaps, support points are needed anyway by both RST and bspline. For moderately large gaps, 2*res or 3*res seems to work fine with bspline. The idea of adding bspline to r.fillnulls is to have a different interpolation method which should give different results, just like in r.resamp.interp. RST tends to provide more detail at the cost of overshoots, whereas bspline with method=bilinear does not produce overshoots and produces a smoother surface. Bspline with method=bicubic can produce results very similar to RST, but AFAIU the idea of bspline interpolation is to avoid overshoots and to produce a smoother surface. Therefore r.fillnulls could just as well have the method options rst,bspline with bspline using bilinear by default. The results would then be intentionally different. I have modified r.fillnulls in trunk to use r.resamp.bspline which is designed for raster interpolation instead v.surf.bspline. That avoids the vectorization step and essentially reduces the processing to r.resamp.bspline -n followed by r.patch. Markus M PS: I fixed r.buffer in trunk, it could now be used (again) by r.fillnulls instead of r.grow. -- Ticket URL: https://trac.osgeo.org/grass/ticket/1088#comment:14 GRASS GIS http://grass.osgeo.org ___ grass-dev mailing list grass-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/grass-dev
[GRASS-dev] Re: [GRASS GIS] #1088: r.fillnulls: support other interpolation methods
#1088: r.fillnulls: support other interpolation methods -+-- Reporter: kyngchaos| Owner: grass-dev@… Type: enhancement | Status: new Priority: normal | Milestone: 6.4.3 Component: Raster | Version: svn-develbranch6 Keywords: r.fillnulls |Platform: All Cpu: All | -+-- Comment(by hamish): Replying to [comment:12 hamish]: see above comments for MarkusM's recommendation for using si=2*res and lambda=0.005. Replying to [comment:14 mmetz]: The sie/sin settings are the result of guestimation and experimenting. any idea about the lambda_i setting? I suggest to change the default to 0.01, but I've only tested with lat/lon srtm data so far. (v.surf.bspline cross- validation output is now a bit easier to digest in 6.5svn, btw; will port to trunk soon) for that data res*3.0 seemed to do well for sie,sin. res*1.0 was not enough for the patched tiles I tried it with. RST tends to provide more detail at the cost of overshoots, fwiw in my tests yesterday v.surf.bspline bicubic with lambda=0.005 and si?=res*3 did the best job, v.surf.bspline with default lambda=1 wasn't totally horrible, but wasn't great, and rst was fairly obvious something had been patched in (but I didn't tune the tension at all). whereas bspline with method=bilinear does not produce overshoots and produces a smoother surface. but the downside of that is that it lowers the contrast/rounds down towards the mean more than a spline fit would.. which is an issue because at least for srtm the holes usually happen in the mountains. fwiw most of the overshoot errors from r.fillnulls are from the masked area afaict, so perhaps can be ignored. Bspline with method=bicubic can produce results very similar to RST, but AFAIU the idea of bspline interpolation is to avoid overshoots and to produce a smoother surface. Therefore r.fillnulls could just as well have the method options rst,bspline with bspline using bilinear by default. The results would then be intentionally different. I personally prefer bicubic to bilinear (will let a job run slower if it means a better end result), and I'd vote to keep method=rst,bicubic,bilinear instead of method= + method2= options. the user doesn't care about which middleware modules, are doing the work, just the beginning and the end. I have modified r.fillnulls in trunk to use r.resamp.bspline which is designed for raster interpolation instead v.surf.bspline. I'll try it out see how it does against the others. That avoids the vectorization step fwiw with -z that wasn't too costly, but happy to try out a new option and see if it does as well as v.surf.bspline's bicubic + lambda=0.005 (which _really_ did nicely). also, what would it take to get 'g.region vect=' to work with a level 1 map? right now you get a no-topology error if you try it with 'r.to.vect -b' ( so make the rst method prep work even smaller/faster/less ram). and essentially reduces the processing to r.resamp.bspline -n followed by r.patch. btw, you'll need to move that del_temp_region back to where it was. with it still in place when you run r.patch you're cropping the patched result with the zoomed-in region. the python lib create temp region fn includes a on-exit hook to automatically unset and delete itself when the module ends, so don't worry about leftovers. PS: I fixed r.buffer in trunk, it could now be used (again) by r.fillnulls instead of r.grow. you mean r.buffer.py :-) r.buffer2 (the classic C version) in trunk always worked bug free for many years, and is about 50x faster (race them with `time`). Perhaps the python version was an old experiment to see if some modules could be combined in a single engine? thanks, Hamish -- Ticket URL: https://trac.osgeo.org/grass/ticket/1088#comment:15 GRASS GIS http://grass.osgeo.org ___ grass-dev mailing list grass-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/grass-dev
[GRASS-dev] Re: [GRASS GIS] #1088: r.fillnulls: support other interpolation methods
#1088: r.fillnulls: support other interpolation methods -+-- Reporter: kyngchaos| Owner: grass-dev@… Type: enhancement | Status: new Priority: normal | Milestone: 6.4.3 Component: Raster | Version: svn-develbranch6 Keywords: r.fillnulls |Platform: All Cpu: All | -+-- Comment(by mmetz): Replying to [comment:15 hamish]: Replying to [comment:12 hamish]: see above comments for MarkusM's recommendation for using si=2*res and lambda=0.005. Replying to [comment:14 mmetz]: The sie/sin settings are the result of guestimation and experimenting. any idea about the lambda_i setting? In my tests, 1 is usually way too high, 0.005 seems to do well in mountainous areas, 0.1 is better for nearly flat areas. RST tends to provide more detail at the cost of overshoots, fwiw in my tests yesterday v.surf.bspline bicubic with lambda=0.005 and si?=res*3 did the best job, v.surf.bspline with default lambda=1 wasn't totally horrible, but wasn't great, and rst was fairly obvious something had been patched in (but I didn't tune the tension at all). I think it is not easy to avoid segmentation artifacts with RST (very high npmin, very low segmax maybe). whereas bspline with method=bilinear does not produce overshoots and produces a smoother surface. but the downside of that is that it lowers the contrast/rounds down towards the mean more than a spline fit would.. Sometimes this is a downside, sometimes not. For hydrological analysis, overshoots are a problem and must be dealt with somehow. Bspline with method=bicubic can produce results very similar to RST, but AFAIU the idea of bspline interpolation is to avoid overshoots and to produce a smoother surface. Therefore r.fillnulls could just as well have the method options rst,bspline with bspline using bilinear by default. The results would then be intentionally different. I personally prefer bicubic to bilinear (will let a job run slower if it means a better end result), and I'd vote to keep method=rst,bicubic,bilinear instead of method= + method2= options. the user doesn't care about which middleware modules, are doing the work, just the beginning and the end. Sounds good. Anyway, the manual should also state that r.fillnulls is an attempt for semi-automated NULL filling, and if in doubt this must be done manually to obtain the desired results. what would it take to get 'g.region vect=' to work with a level 1 map? The routines are already in v.info (trunk). and essentially reduces the processing to r.resamp.bspline -n followed by r.patch. btw, you'll need to move that del_temp_region back to where it was. with it still in place when you run r.patch you're cropping the patched result with the zoomed-in region. Oops. Changed in r50803. The grid geometry of the input raster is still maintained, though. PS: I fixed r.buffer in trunk, it could now be used (again) by r.fillnulls instead of r.grow. you mean r.buffer.py :-) r.buffer2 (the classic C version) in trunk always worked bug free for many years, and is about 50x faster (race them with `time`). Perhaps the python version was an old experiment to see if some modules could be combined in a single engine? Hmm. One of the two r.buffer's should probably be deactivated, after comparative testing with latlong and real projections. Markus M -- Ticket URL: https://trac.osgeo.org/grass/ticket/1088#comment:16 GRASS GIS http://grass.osgeo.org ___ grass-dev mailing list grass-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/grass-dev
[GRASS-dev] Re: [GRASS GIS] #1088: r.fillnulls: support other interpolation methods
#1088: r.fillnulls: support other interpolation methods -+-- Reporter: kyngchaos| Owner: grass-dev@… Type: enhancement | Status: new Priority: normal | Milestone: 6.4.3 Component: Raster | Version: svn-develbranch6 Keywords: r.fillnulls |Platform: All Cpu: All | -+-- Comment(by hamish): Replying to [comment:15 hamish]: any idea about the lambda_i setting? Replying to [comment:16 mmetz]: In my tests, 1 is usually way too high, 0.005 seems to do well in mountainous areas, 0.1 is better for nearly flat areas. in trunk and devbr6 I have changed it to 0.01 now. I would suggest to backport to 6.4svn ASAP before it harms any more data, but would like to test it with some meter-based LIDAR data first, as I've only tried with lat/lon srtm so far. I think it is not easy to avoid segmentation artifacts with RST (very high npmin, very low segmax maybe). for r.fillnulls I'd expect the 3-cell buffer of most holes not to fall on segmentation boundaries. it wasn't that it had obvious segmentation lines in it, rather it just didn't follow the ridges and valleys in such a natural way. (qualitative eyeball analysis) I've been convinced for years that there must be a way to avoid the segmentation artifacts without massively overlapping the quadtree boxes.. still keeping an eye out. Sounds good. Anyway, the manual should also state that r.fillnulls is an attempt for semi-automated NULL filling, and if in doubt this must be done manually to obtain the desired results. yeah, the man page needs quite a bit of work actually.. and I can't quite figure out what the WARNING is trying to get you to delta against with r.mapcalc to verify the interpolation. seems an impossible task, since if you knew the correct answer you there'd be no reason to interpolate. Hmm. One of the two r.buffer's should probably be deactivated, after comparative testing with latlong and real projections. I get 4.3sec for a test in r.buffer.py, 0.140 sec for r.buffer.c, for a simple test starting at a single cell, with visually identical results. Actually I was pleasantly surprised that r.buffer.py did the right thing in lat/lon and made a pear-shaped buffer for the plate carree display, wider nearer to the poles. Hamish -- Ticket URL: https://trac.osgeo.org/grass/ticket/1088#comment:17 GRASS GIS http://grass.osgeo.org ___ grass-dev mailing list grass-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/grass-dev
[GRASS-dev] Re: [GRASS GIS] #1088: r.fillnulls: support other interpolation methods
#1088: r.fillnulls: support other interpolation methods -+-- Reporter: kyngchaos| Owner: grass-dev@… Type: enhancement | Status: new Priority: normal | Milestone: 6.4.3 Component: Raster | Version: svn-develbranch6 Keywords: fillnulls|Platform: All Cpu: All | -+-- Changes (by hamish): * milestone: 6.5.0 = 6.4.3 Comment: * applied in devbr6 by MarkusN in r50114, with sie,sin = 3*res and the default lambda_i=1. * applied to trunk by luca in r50128 with sie,sin = 1*res and the default lambda_i=1. see above comments for MarkusM's recommendation for using si=2*res and lambda=0.005. shall we apply those coeff's? MM: are those numbers going to be dataset dependent, or are they generally good for the way r.fillnulls works by extracting the centers of the 3 cell buffer pixels as the vector starting points? * in 6.5svn I added a step which zooms into the minimum region containing all holes before running v.surf.rst. For my data this made the region about half the size and r.fillnulls ran twice as fast. beforehand v.surf.rst and v.surf.bspline bicubic methods took about the same amount of time (6.5svn versions..), so for my particular srtm data the v.surf.rst method was now twice as fast as the v.surf.bspline bicubic method. I tried running v.surf.bspline with the region zoom but it made no difference besides using 30 subregions (with a few empty) instead of the original 77 subregions (with a lot empty). actually for bspline with the region zoom it took about 10-15% longer, so I only applied it to the RST method, i.e. processing time was mostly independent from region settings. I suspect the original subregion segments were slightly more efficient than the zoomed version due to a bit of luck of where the points were what subregions were empty. both RST and bspline results were (very nearly) identical to without the zoom. nominating both the method= and the RST zoom patch for inclusion in 6.4.3. trialing bspline in grass7 now... rst time was about the same as in devbr6, a wee bit faster due to partial OpenMP support. Hamish -- Ticket URL: https://trac.osgeo.org/grass/ticket/1088#comment:12 GRASS GIS http://grass.osgeo.org ___ grass-dev mailing list grass-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/grass-dev
[GRASS-dev] Re: [GRASS GIS] #1088: r.fillnulls: support other interpolation methods
#1088: r.fillnulls: support other interpolation methods -+-- Reporter: kyngchaos| Owner: grass-dev@… Type: enhancement | Status: new Priority: normal | Milestone: 6.4.3 Component: Raster | Version: svn-develbranch6 Keywords: r.fillnulls |Platform: All Cpu: All | -+-- Changes (by hamish): * keywords: fillnulls = r.fillnulls Comment: trialing bspline in grass7 now... the zoom trick for bspline is the same as in grass65: almost the same but slightly slower. I changed the bspline coeffs to si=3*res and lambda_i = 0.005. Should we change the default option answer to that? I note in the cross-validation code there the maximum they try is 0.05, so maybe 1.0 is way out of useful range. with that, time to complete was not quite as fast as zoomed-RST, but about twice as fast as in grass65. (partly due to si= settings, partly due to using some OpenMP, and perhaps partly due to speedups in trunk that MM was talking about?) both RST and bspline results were (very nearly) identical to without the zoom. (there I meant r.univar results, but really that should be done on the cells-to-be patched in, not the final result where n_patched is quite a small percentage of the whole. an over-zoomed-in d.rast overlay didn't show discernible visual color change) Hamish -- Ticket URL: https://trac.osgeo.org/grass/ticket/1088#comment:13 GRASS GIS http://grass.osgeo.org ___ grass-dev mailing list grass-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/grass-dev
[GRASS-dev] Re: [GRASS GIS] #1088: r.fillnulls: support other interpolation methods
#1088: r.fillnulls: support other interpolation methods -+-- Reporter: kyngchaos| Owner: grass-dev@… Type: enhancement | Status: new Priority: normal | Milestone: 6.5.0 Component: Raster | Version: svn-develbranch6 Keywords: fillnulls|Platform: All Cpu: All | -+-- Comment(by neteler): This patch (I have updated it right now to current SVN) is yet interesting and should IMHO be submitted. -- Ticket URL: http://trac.osgeo.org/grass/ticket/1088#comment:11 GRASS GIS http://grass.osgeo.org ___ grass-dev mailing list grass-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/grass-dev
[GRASS-dev] Re: [GRASS GIS] #1088: r.fillnulls: support other interpolation methods
#1088: r.fillnulls: support other interpolation methods -+-- Reporter: kyngchaos| Owner: grass-...@… Type: enhancement | Status: new Priority: normal | Milestone: 6.5.0 Component: Raster | Version: svn-develbranch6 Keywords: fillnulls|Platform: All Cpu: All | -+-- Comment(by mmetz): Replying to [comment:9 mmetz]: Replying to [comment:8 kyngchaos]: A couple notes/questions: I didn't completely understand the description for the sparse option, but it didn't look like it was needed here. Now it looks to me like it's essentially a mask option, like the maskmap option in v.surf.rst. For now, the mask is implied when the results are patched into the original raster. I guess a mask option (if sparse is such) only makes sense if it reduces processing time by ignoring the masked areas. Right. Unfortunately, using sparse points as input for v.surf.bspline does not reduce processing time substantially, although it could be implemented similar to what I did for r.resamp.bspline. I just noticed that I did implement it, interpolation in v.surf.bspline is skipped if a sparse points vector is given but no sparse points fall into the current subregion. The advantage of the sparse points approach for filling NULL cells would be that more than the currently 3 cells around a gap are used for interpolation, providing (in theory) better results for larger gaps. Markus M -- Ticket URL: https://trac.osgeo.org/grass/ticket/1088#comment:10 GRASS GIS http://grass.osgeo.org ___ grass-dev mailing list grass-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/grass-dev
[GRASS-dev] Re: [GRASS GIS] #1088: r.fillnulls: support other interpolation methods
#1088: r.fillnulls: support other interpolation methods -+-- Reporter: kyngchaos| Owner: grass-...@… Type: enhancement | Status: new Priority: normal | Milestone: 6.5.0 Component: Raster | Version: svn-develbranch6 Keywords: fillnulls|Platform: All Cpu: All | -+-- Comment(by mmetz): In grass 7, it's already implemented in r.resamp.bspline. r.resamp.bspline interpolates a raster to the current resolution, filling NULL cells on the fly, or optionally only interpolates NULL cells. Available methods are bilinear and bicubic, tested with SRTM data, filling large gaps in the European Alps. For grass 6.x, I would use a different approach to r.fillnulls to follow the design of v.surf.bspline, because only interpolating NULL cells is already possible with v.surf.bspline: Create a new raster where NULL cells in the original surface raster are set to something else than NULL, others to NULL. Convert both raster maps to vector points. Use the vector points representing NULL cells as sparse points input for v.surf.bspline. Recommended settings for v.surf.bspline in this case are sie = 2 * ewres, sin = 2 * nsres, lambda = 0.005 (best in my tests, should be somewhere between 0.001 and 0.01). The output raster holds the interpolated NULL cells and can be patched with the original surface. Although all lidar tools work within the current region, none respects a MASK. I'm not sure how this could be implemented properly and efficiently, currently there will be IMHO harmless warnings about no points in the current subregion. Markus M -- Ticket URL: https://trac.osgeo.org/grass/ticket/1088#comment:7 GRASS GIS http://grass.osgeo.org ___ grass-dev mailing list grass-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/grass-dev
[GRASS-dev] Re: [GRASS GIS] #1088: r.fillnulls: support other interpolation methods
Thanks Markus, This is very close to the methods we used in the last month to rebuild early Holocene landscapes in eastern Spain. Your tips are very useful. I only wish we'd had them before we did a lot of experimentation. Cheers Michael __ C. Michael Barton Director, Center for Social Dynamics Complexity Professor of Anthropology, School of Human Evolution Social Change Arizona State University Tempe, AZ 85287-2402 USA voice: 480-965-6262 (SHESC), 480-727-9746 (CSDC) fax: 480-965-7671(SHESC), 480-727-0709 (CSDC) www:http://csdc.asu.edu, http://shesc.asu.edu http://www.public.asu.edu/~cmbarton On Jun 21, 2010, at 9:00 AM, grass-dev-requ...@lists.osgeo.org wrote: Message: 2 Date: Mon, 21 Jun 2010 08:29:31 - From: GRASS GIS t...@osgeo.org Subject: [GRASS-dev] Re: [GRASS GIS] #1088: r.fillnulls: support other interpolation methods To: undisclosed-recipients:; Message-ID: 052.ec39642cf05ae49be4ce809f1417f...@osgeo.org Content-Type: text/plain; charset=utf-8 #1088: r.fillnulls: support other interpolation methods -+-- Reporter: kyngchaos| Owner: grass-...@??? Type: enhancement | Status: new Priority: normal | Milestone: 6.5.0 Component: Raster | Version: svn-develbranch6 Keywords: fillnulls|Platform: All Cpu: All | -+-- Comment(by mmetz): In grass 7, it's already implemented in r.resamp.bspline. r.resamp.bspline interpolates a raster to the current resolution, filling NULL cells on the fly, or optionally only interpolates NULL cells. Available methods are bilinear and bicubic, tested with SRTM data, filling large gaps in the European Alps. For grass 6.x, I would use a different approach to r.fillnulls to follow the design of v.surf.bspline, because only interpolating NULL cells is already possible with v.surf.bspline: Create a new raster where NULL cells in the original surface raster are set to something else than NULL, others to NULL. Convert both raster maps to vector points. Use the vector points representing NULL cells as sparse points input for v.surf.bspline. Recommended settings for v.surf.bspline in this case are sie = 2 * ewres, sin = 2 * nsres, lambda = 0.005 (best in my tests, should be somewhere between 0.001 and 0.01). The output raster holds the interpolated NULL cells and can be patched with the original surface. Although all lidar tools work within the current region, none respects a MASK. I'm not sure how this could be implemented properly and efficiently, currently there will be IMHO harmless warnings about no points in the current subregion. Markus M -- Ticket URL: https://trac.osgeo.org/grass/ticket/1088#comment:7 GRASS GIS http://grass.osgeo.org ___ grass-dev mailing list grass-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/grass-dev
[GRASS-dev] Re: [GRASS GIS] #1088: r.fillnulls: support other interpolation methods
#1088: r.fillnulls: support other interpolation methods -+-- Reporter: kyngchaos| Owner: grass-...@… Type: enhancement | Status: new Priority: normal | Milestone: 6.5.0 Component: Raster | Version: svn-develbranch6 Keywords: fillnulls|Platform: All Cpu: All | -+-- Comment(by kyngchaos): Markus, see my patch. It's very similar to what you describe. A couple notes/questions: I didn't completely understand the description for the sparse option, but it didn't look like it was needed here. Now it looks to me like it's essentially a mask option, like the maskmap option in v.surf.rst. For now, the mask is implied when the results are patched into the original raster. I guess a mask option (if sparse is such) only makes sense if it reduces processing time by ignoring the masked areas. The description says 1 * resolution was good for regularly spaced points, which is what we get from a raster. Why 2* in your method? For the lambda, the docs are unclear what orders of numbers have effects, beyond the larger the number the more smoothing applied. The default is one. You say 0.001-0.01 work well. I figured 1 was something like 1*, or no smoothing. Can you clarify this option? -- Ticket URL: http://trac.osgeo.org/grass/ticket/1088#comment:8 GRASS GIS http://grass.osgeo.org ___ grass-dev mailing list grass-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/grass-dev
[GRASS-dev] Re: [GRASS GIS] #1088: r.fillnulls: support other interpolation methods
#1088: r.fillnulls: support other interpolation methods -+-- Reporter: kyngchaos| Owner: grass-...@… Type: enhancement | Status: new Priority: normal | Milestone: 6.5.0 Component: Raster | Version: svn-develbranch6 Keywords: fillnulls|Platform: Unspecified Cpu: Unspecified | -+-- Comment(by kyngchaos): ...I should make sure it works before attaching the patch... -- Ticket URL: http://trac.osgeo.org/grass/ticket/1088#comment:1 GRASS GIS http://grass.osgeo.org ___ grass-dev mailing list grass-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/grass-dev
[GRASS-dev] Re: [GRASS GIS] #1088: r.fillnulls: support other interpolation methods
#1088: r.fillnulls: support other interpolation methods -+-- Reporter: kyngchaos| Owner: grass-...@… Type: enhancement | Status: new Priority: normal | Milestone: 6.5.0 Component: Raster | Version: svn-develbranch6 Keywords: fillnulls|Platform: Unspecified Cpu: Unspecified | -+-- Comment(by kyngchaos): Now the patch works. I forgot to mention, the idea is that the rst method may not work well on large areas (blocks from segmentation don't align well). bilinear/bicubic seems to process smoother. -- Ticket URL: http://trac.osgeo.org/grass/ticket/1088#comment:2 GRASS GIS http://grass.osgeo.org ___ grass-dev mailing list grass-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/grass-dev
[GRASS-dev] Re: [GRASS GIS] #1088: r.fillnulls: support other interpolation methods
#1088: r.fillnulls: support other interpolation methods -+-- Reporter: kyngchaos| Owner: grass-...@… Type: enhancement | Status: new Priority: normal | Milestone: 6.5.0 Component: Raster | Version: svn-develbranch6 Keywords: fillnulls|Platform: All Cpu: All | -+-- Changes (by hamish): * platform: Unspecified = All * cpu: Unspecified = All Comment: thanks, this has been on the todo list for a long time (I 1/2 thought there was already a ticket for it, but guess not). here's a question (maybe for Markus Metz): what method should be the default in grass 7? porting this patch to python should be trivial. cheers, Hamish -- Ticket URL: http://trac.osgeo.org/grass/ticket/1088#comment:3 GRASS GIS http://grass.osgeo.org ___ grass-dev mailing list grass-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/grass-dev
[GRASS-dev] Re: [GRASS GIS] #1088: r.fillnulls: support other interpolation methods
#1088: r.fillnulls: support other interpolation methods -+-- Reporter: kyngchaos| Owner: grass-...@… Type: enhancement | Status: new Priority: normal | Milestone: 6.5.0 Component: Raster | Version: svn-develbranch6 Keywords: fillnulls|Platform: All Cpu: All | -+-- Comment(by kyngchaos): One thing that would be nice in v.surf.bspline is for it to use the MASK, like v.surf.rst does, otherwise there can be a lot of warnings about no data in a subregion when I want to mask off areas that I want to stay null. (there's a note in the patch about the MASK not applying, but I set it anyways for the future possibility) I was going to add a feature request for this, but I'm burnt out today from testing the processing of the patch. -- Ticket URL: http://trac.osgeo.org/grass/ticket/1088#comment:4 GRASS GIS http://grass.osgeo.org ___ grass-dev mailing list grass-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/grass-dev
[GRASS-dev] Re: [GRASS GIS] #1088: r.fillnulls: support other interpolation methods
#1088: r.fillnulls: support other interpolation methods
-+--
Reporter: kyngchaos| Owner: grass-...@…
Type: enhancement | Status: new
Priority: normal | Milestone: 6.5.0
Component: Raster | Version: svn-develbranch6
Keywords: fillnulls|Platform: All
Cpu: All |
-+--
Comment(by hamish):
MASKs used in GRASS are only used when _reading_ existing raster maps.
They do not mask writing. If v.surf.rst respects MASKs at all it would be
somewhat surprising news to me and should be well documented.
As the output raster map is only written for the current region, ignoring
points outside of this box (and maybe a little buffer around it) can save
some needless work, and so v.surf.rst does do that.
Hamish
ps- note to self:
{{{
- g.rename /dev/null
+ g.rename --quiet
}}}
--
Ticket URL: http://trac.osgeo.org/grass/ticket/1088#comment:5
GRASS GIS http://grass.osgeo.org
___
grass-dev mailing list
grass-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/grass-dev
[GRASS-dev] Re: [GRASS GIS] #1088: r.fillnulls: support other interpolation methods
#1088: r.fillnulls: support other interpolation methods -+-- Reporter: kyngchaos| Owner: grass-...@… Type: enhancement | Status: new Priority: normal | Milestone: 6.5.0 Component: Raster | Version: svn-develbranch6 Keywords: fillnulls|Platform: All Cpu: All | -+-- Comment(by kyngchaos): Sorry, I was being vague, v.surf.rst has an option to set a mask for writing (though I didn't realize that MASKs only apply to reading, thanks, though it makes sense). -- Ticket URL: http://trac.osgeo.org/grass/ticket/1088#comment:6 GRASS GIS http://grass.osgeo.org ___ grass-dev mailing list grass-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/grass-dev
