[GRASSLIST:4844] Re: Problems importing Data into lat/lon with negtiv south value

[Glynn Clements]

Paul Kelly wrote:

> >> This doesn't really explain what looks like a bug, but I think it is a
> >> side-effect of GDAL trying to set the current region to match the imported
> >> data, which IMHO I don't think it should do. E.g. what if you were running
> >> an analysis in the background for which the region set was very important,
> >> and then imported a new file and r.in.gdal suddenly changed the current
> >> region?
> >
> > r.in.gdal doesn't modify the WIND file. If you use the -e flag, it
> > modifies the DEFAULT_WIND file in the PERMANENT mapset[1], otherwise
> > it only modifies the process' internal region.
> 
> I don't disagree that this is the actual behaviour, but on reading though 
> the source code it seems to me this wasn't the intended behaviour and it 
> is trying to write out the bounds of the imported map as the current 
> region (i.e WIND). But G_set_window() doesn't change anything because of 
> some caching it does or something (I looked into it a long time ago and 
> don't remember all the details now).
> 
> Unless I'm misunderstanding what G_set_window() is supposed to do but I 
> think it's for setting the current region.

G_set_window() sets the process' internal region (i.e. G__.window),
not the WIND file. It doesn't modify the WIND file, and isn't meant
to.

If you want to set the mapset's current region (i.e. the WIND file),
you need to use G_put_window(). If you want to modify some other file
(i.e. PERMANENT/DEFAULT_WIND or a cellhd file), you need to use
G__put_window().

Quite a few programs use G_set_window() to set a window other than
that defined in the WIND file.

E.g. programs which operate on the raw raster data without clipping,
padding or resampling (e.g. r.proj, r.bilinear) copy the raster's
cellhd to the current region:

	G_get_cellhd(mapname, mapset, &window);
	G_set_window(&window);

Similarly, import (r.in.*) programs set up a region with the same
number of rows and columns as the raster being imported. If
georeferencing information is available, they will set the boundaries
from that, otherwise they typically place the south-west corner at
(0,0) and the north-east corner at (W,H) (i.e. nsres == ewres == 1).

However, only those programs which wish to make persistent changes
(e.g. g.region, d.zoom) use G_put_window().

I'm fairly certain that r.in.gdal's behaviour (both with and without
the -e flag) is intentional.

> >> I think it should only do this if the -e flag to extend the region is
> >> given, or perhaps there should be a new flag specifically to set the
> >> region to match the imported file.
> >
> > I'm not sure that there would be much point, as you can just use
> > "g.region rast=..." if you want that behaviour.
> 
> Yes I agree, except from what I've noticed it seems r.in.gdal was designed 
> to actually do this. But maybe these lines should really be deleted then:
> 
> >> If you comment out the lines in raster/r.in.gdal/main.c:
> >>      if(G_set_window (&cellhd) < 0)
> >>          exit(3);

No. Any raster import program has to at least set the window such that
the rows/cols values match the actual number of rows and columns. 
Either that, or it will have to manually resample the image.

G_put_raster_row() etc don't do any resampling. They expect the buffer
to contain <cols> values, and they expect to be called exactly <rows>
times.

Most raster programs take their region (and hence grid dimensions)
from the WIND file, but import programs take the grid dimensions from
the file being imported, and so have to set the region to match.

> >> and re-compile does it work?
> >> Although seeing that doesn't fix the underlying bug, it might run into
> >> problems later I suppose.
> >>
> >> Anybody else have an opinion on the region setting and extending behaviour
> >> of r.in.gdal and v.in.ogr?
> >
> > If r.in.gdal is importing a georeferenced raster, it should attempt to
> > set the correct bounds and resolution (i.e. cellhd) for that raster,
> > rather than using (0,0,W,H). It shouldn't attempt to set the current
> > region or the default region; that should be left to g.region (which
> > should have an option to set the default region).
> >
> > I don't think that this is an issue for v.in.ogr; AFAIK, vector maps
> > don't have a region (cellhd).
> 
> If you are using v.in.ogr to create a new location based on the projection
> of the vector map being imported, it has to create a DEFAULT_WIND, and it is
> good to populate it with meaningful region extents, which can be got through
> some OGR function (get extents or something like that).

Right; if you're creating a new location, it makes sense.

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