[GRASS-dev] latlon support

GRASS GRASS developer
1

When I was updating v.in.gshhs (see user list) I came across a strange phenomenon in GRASS latlon support when crossing the datum border. The map extends reported by v.info are wrong because East is larger than 180 (about 191) and converted to (East - 360) which is only 10 degrees east of West. The shorelines are spanning 370 degrees of longitude and not 10 degrees. I know GRASS doesn't support crossing the datum border, i.e. going beyond 180E or 180W, but I don't understand why not. The GSHHS shoreline vector for example crosses the datum border and clipping or converting coordinates of vertices would create artifacts. Geographic features can cross the datum border, so why not supporting a range -360 to 360 for longitude? I am aware that 181 is equal to -179 (in latlon), but going from 179 to 181 is different from going from 179 to -179. Mathematically the difference in the first case is 2 degrees longitude, in the latter case 358 degrees longitude. The distance between two points depends on whether you go east or west (or north or south). So if the first point is on 179 lon and the second is on 181 lon, this would indicate to go east, while when the second point is on -179 lon this would indicate to go west from the first point. Does that make sense?
How much havoc would it cause to support longitudes from -360 to 360 in latlon?

Markus M

2

Markus Metz wrote:

When I was updating v.in.gshhs (see user list) I came across a strange
phenomenon in GRASS latlon support when crossing the datum border. The
map extends reported by v.info are wrong because East is larger than 180
(about 191) and converted to (East - 360) which is only 10 degrees east
of West. The shorelines are spanning 370 degrees of longitude and not 10
degrees. I know GRASS doesn't support crossing the datum border, i.e.
going beyond 180E or 180W, but I don't understand why not. The GSHHS
shoreline vector for example crosses the datum border and clipping or
converting coordinates of vertices would create artifacts. Geographic
features can cross the datum border, so why not supporting a range -360
to 360 for longitude? I am aware that 181 is equal to -179 (in latlon),
but going from 179 to 181 is different from going from 179 to -179.
Mathematically the difference in the first case is 2 degrees longitude,
in the latter case 358 degrees longitude. The distance between two
points depends on whether you go east or west (or north or south). So
if the first point is on 179 lon and the second is on 181 lon, this
would indicate to go east, while when the second point is on -179 lon
this would indicate to go west from the first point. Does that make sense?
How much havoc would it cause to support longitudes from -360 to 360 in
latlon?

The biggest problem is likely to be that, once you open this can of
worms, you'll just run away in horror. I know I did.

IMHO, the optimal solution for lat/lon is to remove the region
constraints, and require the data to be "Euclidified" either on input
or by the low-level retrieval functions (e.g. G_get_raster_row etc).

Higher-level code would then treat all coordinates as conforming to
the rules of Euclidean geometry.

E.g. you could freely set the region to west=-360, east=+360, and get
two copies of the globe. The points 0,-360 and 0,+360 would be treated
as being separated by twice the equatorial circumference. And so on.

Expecting every piece of code which deals with coordinates to
explicitly handle both Euclidean and non-Euclidean cases is asking too
much, IMHO. Actually, not IMHO; you just have to look at the code to
see how often it doesn't happen.

E.g. although the new display library is based around the same path
construction primitives as PostScript, PDF, and cairo, it can't just
feed the coordinates directly to the driver.

Instead, it has to store the entire path locally, so that it can draw
multiple copies (once with the RHS of the shape appearing on the LHS
of the display, and once with the LHS of the shape appearing on the
RHS of the display).

--
Glynn Clements <glynn@gclements.plus.com>

3

Glynn Clements wrote:

Markus Metz wrote:

When I was updating v.in.gshhs (see user list) I came across a strange phenomenon in GRASS latlon support when crossing the datum border. The map extends reported by v.info are wrong because East is larger than 180 (about 191) and converted to (East - 360) which is only 10 degrees east of West. The shorelines are spanning 370 degrees of longitude and not 10 degrees. I know GRASS doesn't support crossing the datum border, i.e. going beyond 180E or 180W, but I don't understand why not. The GSHHS shoreline vector for example crosses the datum border and clipping or converting coordinates of vertices would create artifacts. Geographic features can cross the datum border, so why not supporting a range -360 to 360 for longitude? I am aware that 181 is equal to -179 (in latlon), but going from 179 to 181 is different from going from 179 to -179. Mathematically the difference in the first case is 2 degrees longitude, in the latter case 358 degrees longitude. The distance between two points depends on whether you go east or west (or north or south). So if the first point is on 179 lon and the second is on 181 lon, this would indicate to go east, while when the second point is on -179 lon this would indicate to go west from the first point. Does that make sense?
How much havoc would it cause to support longitudes from -360 to 360 in latlon?
    
The biggest problem is likely to be that, once you open this can of
worms, you'll just run away in horror. I know I did.
  

OK, I won't start hacking around here :slight_smile: I thought about it a bit more and it really is a big mess.

IMHO, the optimal solution for lat/lon is to remove the region
constraints, and require the data to be "Euclidified" either on input
or by the low-level retrieval functions (e.g. G_get_raster_row etc).
  

As long as the input is ok, I don't see a problem with GRASS modules, in theory the output should conform as well to Euclidean geometry as the input. But here the problem starts. I get all sorts of dirty data, in particular shapefiles, that I am supposed to work with, and cleaning them up is sometimes really a pain... After cleaning I have topologically correct vectors, but that doesn't guarantee that they are "Euclidified". If I understand you right I would have to check for every single line and every single boundary if every single vertex is correct. In practice I can only do an educated guess if, to stick to my previous example, it is correct that one point has lon 179 and the next point has lon 181, i.e. crossing the datum border is right, or if it is correct to go instead from lon 179 to lon -179, i.e. go over Greenwich meridian and not over the datum border. I guess that for most geographic features it will be correct to cross the datum border, as for the Siberian shorelines, and GRASS seems to do so (I haven't tested everything, just v.info, the display, and "Measure distance" in the display manager).
I have changed v.in.gshhs accordingly and both v.info and the display (wxgui) work as expected (correct).

Higher-level code would then treat all coordinates as conforming to
the rules of Euclidean geometry.
  

I think this is the only way, higher-level code must assume that the coordinates given by the input are correct. The alternative is IMHO out of discussion (trying to guess if coordinates should be converted or not).

E.g. you could freely set the region to west=-360, east=+360, and get
two copies of the globe. The points 0,-360 and 0,+360 would be treated
as being separated by twice the equatorial circumference. And so on.
  

That can be intended but I am afraid not all people generating spatial data think about this effect. It may be safer to convert coordinates to the 180 deg lon and 90 deg lat limits and take the shortest way between two points.

Expecting every piece of code which deals with coordinates to
explicitly handle both Euclidean and non-Euclidean cases is asking too
much, IMHO. Actually, not IMHO; you just have to look at the code to
see how often it doesn't happen.

E.g. although the new display library is based around the same path
construction primitives as PostScript, PDF, and cairo, it can't just
feed the coordinates directly to the driver.

Instead, it has to store the entire path locally, so that it can draw
multiple copies (once with the RHS of the shape appearing on the LHS
of the display, and once with the LHS of the shape appearing on the
RHS of the display).

Not afraid of worms but of horror,

Markus M

4

Markus Metz wrote:

> IMHO, the optimal solution for lat/lon is to remove the region
> constraints, and require the data to be "Euclidified" either on input
> or by the low-level retrieval functions (e.g. G_get_raster_row etc).

As long as the input is ok, I don't see a problem with GRASS modules, in
theory the output should conform as well to Euclidean geometry as the
input. But here the problem starts. I get all sorts of dirty data, in
particular shapefiles, that I am supposed to work with, and cleaning
them up is sometimes really a pain... After cleaning I have
topologically correct vectors, but that doesn't guarantee that they are
"Euclidified". If I understand you right I would have to check for every
single line and every single boundary if every single vertex is correct.
In practice I can only do an educated guess if, to stick to my previous
example, it is correct that one point has lon 179 and the next point has
lon 181, i.e. crossing the datum border is right, or if it is correct to
go instead from lon 179 to lon -179, i.e. go over Greenwich meridian and
not over the datum border.

By "Euclidified", I mean that the data is in a form such that code
which process it can treat -179 and +181 as being distinct points 360
units apart.

E.g. for a global map 360 degrees wide and centred on Greenwich,
Euclidification would involve either copying or moving the Chukchi
Peninsula over to the left-hand side of the map, so that it doesn't
get chopped off by the right-hand edge.

Similarly, to use r.resamp.interp on a global DEM without getting a
gap at the 180th meridian, the longitude bounds would need to be set
to 360 degrees plus 4 cells, so that the cells around the 180th
meridan appear on both sides of the map, as the code will assume that
the columns either side of column n are simply n+1 and n-1.

[r.resamp.interp can perform the enlargement itself, but libgis has to
permit it, and not insist that -182..+182 means +178..+182.]

> E.g. you could freely set the region to west=-360, east=+360, and get
> two copies of the globe. The points 0,-360 and 0,+360 would be treated
> as being separated by twice the equatorial circumference. And so on.

That can be intended but I am afraid not all people generating spatial
data think about this effect. It may be safer to convert coordinates to
the 180 deg lon and 90 deg lat limits and take the shortest way between
two points.

My point is that this needs to be done in the lowest-level libraries,
not in higher-level libraries or the modules.

High-level code must be free to assume that eastings increase
monotonically from left to right, that if a point has an easting of x,
then the point which is d units to its east has an easting of x+d, and
so on.

Requiring code to explicitly use specialised (spherical-aware)
functions in place of +,-,<,>,== etc is not a realistic strategy.

--
Glynn Clements <glynn@gclements.plus.com>

5

Glynn Clements wrote:

By "Euclidified", I mean that the data is in a form such that code
which process it can treat -179 and +181 as being distinct points 360
units apart.
  

It seems that grass treats -179 and 179 as being 2 units apart, I do not have to change it to 179 and 181. Latlon coords are apparently not always but sometimes treated as not being "Euclidified", when getting map extends they are forced into the 180 degree limits, when writing out a vector feature they are left as they are.

  

... It may be safer to convert coordinates to the 180 deg lon and 90 deg lat limits and take the shortest way between two points.
    
My point is that this needs to be done in the lowest-level libraries,
not in higher-level libraries or the modules.
  

I had to do that in the module, the vector libraries do not convert latlon coords to these limits although these limits are required by other library code. You can write a vector in latlon with whatever coordinates you like, but such a vector may not be processed properly. This does not apply to raster maps because these are (almost) always restricted to the current region settings.

High-level code must be free to assume that eastings increase
monotonically from left to right, that if a point has an easting of x,
then the point which is d units to its east has an easting of x+d, and
so on.
  

This is currently not the case with library code, at least with regard to map extends. If the Western edge of a map is at -179 and the Eastern edge is at +181, the map extend is not calculated as (East - West) but as (East - 360 - West).

Requiring code to explicitly use specialised (spherical-aware)
functions in place of +,-,<,>,== etc is not a realistic strategy.
  

It seems the display does so when zooming out or panning. The Chukchi Peninsula and other features are copied around and displayed on both sides, I haven't done that in the vector, there all features are present only on one side. This is not a complaint about the display, I like it that way.

I think my questioning started because the region settings including map extends are restricted to the 180 degree lon and 90 degree lat limits whereas vector operations can exceed these limits causing problems later on. To make grass handling of latlon coordinates consistent all around, either the limits on region settings must be relaxed (taking care of all the associated risks) or latlon vector coords must be forced into the latlon limits by some lower-level handling.

The currently allowed mix is confusing.

6

Markus Metz wrote:

> Requiring code to explicitly use specialised (spherical-aware)
> functions in place of +,-,<,>,== etc is not a realistic strategy.

It seems the display does so when zooming out or panning. The Chukchi
Peninsula and other features are copied around and displayed on both
sides,

This is currently done in the display library (lib/display/draw2.c).

I haven't done that in the vector, there all features are present
only on one side. This is not a complaint about the display, I like it
that way.

I think my questioning started because the region settings including map
extends are restricted to the 180 degree lon and 90 degree lat limits
whereas vector operations can exceed these limits causing problems later
on. To make grass handling of latlon coordinates consistent all around,
either the limits on region settings must be relaxed (taking care of all
the associated risks) or latlon vector coords must be forced into the
latlon limits by some lower-level handling.

Alas, it isn't just the coordinates, but the topology. E.g. for
Antartica, you have to add south, east and west edges to get a closed
region in Euclidean space.

The currently allowed mix is confusing.

Indeed.

--
Glynn Clements <glynn@gclements.plus.com>

7

Glynn Clements wrote:

Markus Metz wrote:

I think my questioning started because the region settings including map extends are restricted to the 180 degree lon and 90 degree lat limits whereas vector operations can exceed these limits causing problems later on. To make grass handling of latlon coordinates consistent all around, either the limits on region settings must be relaxed (taking care of all the associated risks) or latlon vector coords must be forced into the latlon limits by some lower-level handling.
    
Alas, it isn't just the coordinates, but the topology. E.g. for
Antartica, you have to add south, east and west edges to get a closed
region in Euclidean space.
  

I think the poles and the datum border are two different problems, correct me if I'm wrong.
Using the GSHHS shorelines, topology seems to be ok for Eurasiafrica when I create a proper area with centroid. At least building topology does not report errors and the display is ok. Not so for Antarctica, both topology and display are wrong. For Antarctica a different projection is needed.
A rather crude description of what we see when we use latlon is a plane created by putting an axis through the poles and then folding the globe open along the datum border. Therefore I guess it is less of a problem to traverse the datum border because the axis has been set first, then the folding-up longitude. Maybe the location of the axis is crucial and determines the degrees of freedom for moving around in this pseudo-projection (only the latter one of two determinants allows changing, the one set last, longitude).

  

The currently allowed mix is confusing.
    
Indeed.
  
How about reversing the current policy? Restrict both raster and vector maps to the 180 degree lon and 90 degree lat limit, that avoids duplicate features within a map. Allow the region to be set anywhere between -360 and 360 lon but restrict the width of the computational region to 360 degrees lon to avoid duplicates. Wrap map contents to the current computational region if any map contents exist that fall into the computational region. This would allow traversing the datum border without gaps or breaks, e.g. for raster resampling. The display is not affected, apart from "zoom to the current computational region".

I'm afraid you have thought of that possibility before and found flaws in it. Apart from rewriting a lot of code...

8

Markus Metz wrote:

How about reversing the current policy? Restrict both raster and vector
maps to the 180 degree lon and 90 degree lat limit, that avoids
duplicate features within a map. Allow the region to be set anywhere
between -360 and 360 lon but restrict the width of the computational
region to 360 degrees lon to avoid duplicates.

That is problematic for computations which use a window (e.g.
r.neighbors, r.proj, r.mfilter etc) on a whole-earth map. For that,
you need 360 degrees plus the width of the window.

[The same issue applies at the poles, but people don't normally care
about the poles if they're working in lat/lon.]

Wrap map contents to the
current computational region if any map contents exist that fall into
the computational region. This would allow traversing the datum border
without gaps or breaks, e.g. for raster resampling. The display is not
affected, apart from "zoom to the current computational region".

I'm afraid you have thought of that possibility before and found flaws
in it. Apart from rewriting a lot of code...

AFAICT, the most important case to solve is when the region is set to
the Bering strait, and a vector map has Alaska at ~180W and Sibera at
~180E. Some of the data will need a 360-degree shift to obtain the
correct relative position.

--
Glynn Clements <glynn@gclements.plus.com>

9

Glynn Clements wrote:

AFAICT, the most important case to solve is when the region is set to
the Bering strait, and a vector map has Alaska at ~180W and Sibera at
~180E. Some of the data will need a 360-degree shift to obtain the
correct relative position.

the same problem applies to data which is not land based!
fwiw, my longitude >170 and in my oceanography work I wrap past 180 on
a fairly regular basis. (~90% is raster for me, certainly others work
with point data, ship tracks, 200nm buffer zones, and sonar swaths)

I fully understand that it is a total pain in the neck to program the
longitude wrap-around successfully, but to the end users who have to
deal with it it's a real godsend and a real draw-card for GRASS. Not
much other GIS software can do it nearly as well as GRASS 5/6's raster
map libs, Xmonitors, d.zoom, etc. And so I would argue it is well worth
the (albeit large) effort to keep the libraries wrapping. The current
raster map code works remarkably well for most tasks.

see also the earlier bug reports for lat/lon wrap support for the vector
libs and r.resamp.(?? whatever one Glynn was fighting with some months
ago which exposed the neighbors' complication).

Markus Metz wrote:

Restrict both raster and vector maps to the 180 degree lon

note that the planetary science people prefer to keep lon as 0-360,
not -180 thru 180. (search ML archives for "Mars")

Also, AFAIR, for vector maps d.what.vect report correct land areas for
lat/lon polygons. So the vector lib isn't completely lat/lon braindead.

Hamish

10

Glynn Clements wrote:

Markus Metz wrote:

How about reversing the current policy? Restrict both raster and vector maps to the 180 degree lon and 90 degree lat limit, that avoids duplicate features within a map. Allow the region to be set anywhere between -360 and 360 lon but restrict the width of the computational region to 360 degrees lon to avoid duplicates.
    
That is problematic for computations which use a window (e.g. r.neighbors, r.proj, r.mfilter etc) on a whole-earth map. For that,
you need 360 degrees plus the width of the window.
  

Do you really mean "360 degrees plus the width of the window"? That would be max 720 degrees? The status quo is to restrict the width of the current window to 360 degrees, maximum possible extends are 180W to 180E. So you say the current status quo is problematic? This limit applies to the computational region, a module is free to override that for internal calculations, but as soon as a module calls G_put_raster_row(), the current computational region is respected, right?

  

Wrap map contents to the current computational region if any map contents exist that fall into the computational region. This would allow traversing the datum border without gaps or breaks, e.g. for raster resampling. The display is not affected, apart from "zoom to the current computational region".

I'm afraid you have thought of that possibility before and found flaws in it. Apart from rewriting a lot of code...
    
AFAICT, the most important case to solve is when the region is set to
the Bering strait, and a vector map has Alaska at ~180W and Sibera at
~180E. Some of the data will need a 360-degree shift to obtain the
correct relative position.

That is apparently done not only for display but also for vector operations.

Anyway, I just realized (embarrassing that I did so only now) that I can set the window to East = 170W and West = 170E which gives me a window width of 20 degrees going over the datum border. There is no error that East must be east of West, this is a valid window setting. A vector is properly displayed. I can convert that vector to a raster, that raster is also properly displayed. The raster map extends are East = 170W and West = 170E, seamlessly going over the datum border, great!
East = 190E and West = 170E is not a valid window setting though, g.region exits with an error.
So with the current latlon handling of GRASS everything is just fine as long as the map coordinates are within limits. My only complaint would now be that somewhere deep down it should be enforced for vectors to stick to these limits. Further on, if there are any modules that require East to be east of West they do so wrongly for latlon, this is correctly handled by g.region.

Replying to Hamish:

Markus Metz wrote:

> > Restrict both raster and vector maps to the 180 degree lon

> note that the planetary science people prefer to keep lon as 0-360,
> not -180 thru 180. (search ML archives for "Mars")

The status quo is that both raster and vector maps must stick to the 180 degree lon limit, this is required by GRASS, it seems. At least everything I tried works fine when keeping vector lon values to the 180 limit.

> Also, AFAIR, for vector maps d.what.vect report correct land areas for
> lat/lon polygons. So the vector lib isn't completely lat/lon braindead.

Yes, but playing around here with world shorelines gives me the impression that vector operations work better when the coords stick to the limits. The libraries really do the wrapping very well!

There is however a problem with the new vdigit, lines and borders are not displayed properly when the window is set to East = 170W and West = 170E. The old tcltk v.digit displays lines and borders properly.

11

Markus Metz wrote:

>> How about reversing the current policy? Restrict both raster and vector
>> maps to the 180 degree lon and 90 degree lat limit, that avoids
>> duplicate features within a map. Allow the region to be set anywhere
>> between -360 and 360 lon but restrict the width of the computational
>> region to 360 degrees lon to avoid duplicates.
>>
>
> That is problematic for computations which use a window (e.g.
> r.neighbors, r.proj, r.mfilter etc) on a whole-earth map. For that,
> you need 360 degrees plus the width of the window.
>
Do you really mean "360 degrees plus the width of the window"? That
would be max 720 degrees?

Wrong window :wink:

I'm talking about the buffer/kernel/etc "window", i.e. the size=
parameter for r.neighbors, 4 cells for r.proj method=cubic, etc.

E.g. for r.proj method=cubic, you need 360 degrees plus 4 cells.

The status quo is to restrict the width of the
current window to 360 degrees, maximum possible extends are 180W to
180E. So you say the current status quo is problematic?

It means that modules which operate upon a moving window cannot create
output which fully covers 360 degrees. Such modules need a full window
of non-null cells to generate a non-null result, so they leave a
half-window border of null cells on each side of the region where the
window bounds stray outside of the region bounds.

This limit
applies to the computational region, a module is free to override that
for internal calculations, but as soon as a module calls
G_put_raster_row(), the current computational region is respected, right?

The region used by the libraries defaults to the contents of the WIND
file (or WIND_OVERRIDE or GRASS_REGION), but can be changed by the
code. Code can even change it back and forth while operating, e.g.
r.resamp.{interp,stats} set the region to the input map's native
resolution for reading and the current (WIND) resolution for writing.

It's debatable whether you would want to do this in the case where
you're enlarging the window for the above reasons. It would be simpler
and more efficient to just generate a map which is slightly larger
than the current region, given that the excess will consist entirely
of null cells.

--
Glynn Clements <glynn@gclements.plus.com>

12

Glynn Clements wrote:

Markus Metz wrote:

Do you really mean "360 degrees plus the width of the window"? That would be max 720 degrees?
    
Wrong window :wink:

I'm talking about the buffer/kernel/etc "window", i.e. the size=
parameter for r.neighbors, 4 cells for r.proj method=cubic, etc.

E.g. for r.proj method=cubic, you need 360 degrees plus 4 cells.

Oh, that window :slight_smile:

The status quo is to restrict the width of the current window to 360 degrees, maximum possible extends are 180W to 180E. So you say the current status quo is problematic?
    
It means that modules which operate upon a moving window cannot create
output which fully covers 360 degrees. Such modules need a full window
of non-null cells to generate a non-null result, so they leave a
half-window border of null cells on each side of the region where the
window bounds stray outside of the region bounds.
  

Two stupid questions:
1) If the resampling method only needs the values of the moving window, can't you use a temporary array where you feed all the cell values in and assume NULL when no cell value is available (e.g. outside current region window)? Such a temporary array would not be bothered with the region window, resampling can be done in that array and the resulting value for the target cell put back to the array/structure used for writing the output map.
2) If the resampling method needs the full map in memory (also applies to segmented mode), I think RST needs that, I have the same question as for question 1) only that now the array needs to have (region window + resampling window) size in both dimensions (rows and columns).
Sorry, I haven't followed your work on resampling, so I don't know if this can work at all.

The region used by the libraries defaults to the contents of the WIND
file (or WIND_OVERRIDE or GRASS_REGION), but can be changed by the
code. Code can even change it back and forth while operating, e.g. r.resamp.{interp,stats} set the region to the input map's native
resolution for reading and the current (WIND) resolution for writing.
  

This is still what I mean with computational window, be it set by the user or by the module. I haven't looked it up, but assumed that when a module calls G_set_window(), also when using WIND_OVERRIDE, G_set_window() checks if the requested window is legal, in case of latlon if it is within the 180 lon and 90 lat limits.

I know you want to get r.resamp.{interp,stats} right, but I'm really more interested in the vector handling. Because 1) I think you will succeed and I'm far from knowing enough to help you in any way with that and 2) the vectors are a bit neglected since Radim Blazek left (I'm not stepping on anybody's feet?) so I start nagging a bit :wink: Here too I don't know enough to make useful (precise) suggestions for the libraries, but I'm looking for opportunities to improve them, and this latlon 180 degree limit is one such opportunity. In my experience vector operations really work much better when the coords stick to these limits. My initial complaint about the datum border is obsolete, I learned that the datum border is not a problem, for both raster and vector operations.

Markus M

13

On Thu, Jan 29, 2009 at 10:37 PM, Markus Metz
<markus.metz.giswork@googlemail.com> wrote:

I know you want to get r.resamp.{interp,stats} right, but I'm really more
interested in the vector handling. Because 1) I think you will succeed and
I'm far from knowing enough to help you in any way with that and 2) the
vectors are a bit neglected since Radim Blazek left (I'm not stepping on
anybody's feet?)

You are most welcome! The vector part needs definitely a new boost.

thanks for your efforts,
Markus

14

Markus Metz wrote:

> It means that modules which operate upon a moving window cannot create
> output which fully covers 360 degrees. Such modules need a full window
> of non-null cells to generate a non-null result, so they leave a
> half-window border of null cells on each side of the region where the
> window bounds stray outside of the region bounds.
>
Two stupid questions:

1) If the resampling method only needs the values of the moving window,
can't you use a temporary array where you feed all the cell values in
and assume NULL when no cell value is available (e.g. outside current
region window)? Such a temporary array would not be bothered with the
region window, resampling can be done in that array and the resulting
value for the target cell put back to the array/structure used for
writing the output map.

For a 360-degree lat/lon map, you don't want to pad with nulls, you
want to use the wrapped cells.

E.g. for r.neighbors size=5 with ewres=1, you want 364 columns, so
that you get 360 valid columns.

--
Glynn Clements <glynn@gclements.plus.com>