[GRASS-dev] Re: [GRASS-user] mask in PERMANENT affect on other mapsets during reprojection

GRASS GRASS developer
1

[CC to grass-dev]

Jarekj wrote:

The situation is as follows:
In location GRS mapset PERMANET I have MASK (with projection of 2180
based on GRS elipsoid)
in second mapset let say jarekjaI have map which is to be reprojected to
other location (WGS). In that mapset no MASK is present.
In location WGS mapset PERMANENT I reporject map from location GRS and
I have only fragment allowed by MASK from GRS PERMANENT It is intended
or some bug?

The MASK file has to exist in the current mapset in the current
location.

However: r.proj keeps switching the current location between the two
locations. It switches to the source location when opening and reading
the input map, and switches back to the target location when opening
and writing the output map.

Where it looks for the mask is determined by which is the current one
at the point that it looks. This point isn't well defined; it normally
happens as a side effect of some other function.

AFAICT, it uses the source location when looking for the mask.
However, it doesn't change the mapset; i.e. it uses the current mapset
(not the mapset= value), which may not even exist in the source
location.

This is quite clearly a bug; it should switch mapsets when it switches
locations. At present the code only uses the mapset= value when
opening the input map, overlooking the fact that it is also used
internally, e.g. when looking for the mask.

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

2

Jarekj wrote:

The situation is as follows:
In location GRS mapset PERMANET I have MASK (with projection of 2180
based on GRS elipsoid)
in second mapset let say jarekjaI have map which is to be reprojected to
other location (WGS). In that mapset no MASK is present.
In location WGS mapset PERMANENT I reporject map from location GRS and
I have only fragment allowed by MASK from GRS PERMANENT It is intended
or some bug?

side-suggestion: before running r.proj run v.in.region in the source location then v.proj from the taget location. then
"g.region vect=region_box -p" and check that the new resolution is ok and that rows x cols is no less than in the source region, slightly more if the reprojection rotates the box significantly. and/or try the "-n" flag with r.proj to see if that makes any help.

Hamish

3

I always do as you suggest (in general other methods (with cs2cs) in grass are to complicated), but it not concern the problem I mentioned
Well, I found this during testing prototype of script to reproject any raster in one step (without changing locations) so I must to know If it is intended (for reasons unknown for me) and Script must additionally check if no mask is present in PERMANENT or it is a bug and will be removed in the feature
sorry I didn't add grass version: latest grass7.0 svn

Jarek

Hamish pisze:

Jarekj wrote:
  

The situation is as follows:
In location GRS mapset PERMANET I have MASK (with projection of 2180 based on GRS elipsoid)
in second mapset let say jarekjaI have map which is to be reprojected to other location (WGS). In that mapset no MASK is present.
In location WGS mapset PERMANENT I reporject map from location GRS and I have only fragment allowed by MASK from GRS PERMANENT It is intended or some bug?
    
side-suggestion: before running r.proj run v.in.region in the source location then v.proj from the taget location. then
"g.region vect=region_box -p" and check that the new resolution is ok and that rows x cols is no less than in the source region, slightly more if the reprojection rotates the box significantly. and/or try the "-n" flag with r.proj to see if that makes any help.

Hamish

4

So as I understand the problem will heppen when I have mask in mapset in source which is named as my target mapset in target location: ussually it will concern PERMANENT?

Thanks for information
Jarek

Glynn Clements pisze:

[CC to grass-dev]

Jarekj wrote:

The situation is as follows:
In location GRS mapset PERMANET I have MASK (with projection of 2180 based on GRS elipsoid)
in second mapset let say jarekjaI have map which is to be reprojected to other location (WGS). In that mapset no MASK is present.
In location WGS mapset PERMANENT I reporject map from location GRS and I have only fragment allowed by MASK from GRS PERMANENT It is intended or some bug?
    
The MASK file has to exist in the current mapset in the current
location.

However: r.proj keeps switching the current location between the two
locations. It switches to the source location when opening and reading
the input map, and switches back to the target location when opening
and writing the output map.

Where it looks for the mask is determined by which is the current one
at the point that it looks. This point isn't well defined; it normally
happens as a side effect of some other function.

AFAICT, it uses the source location when looking for the mask. However, it doesn't change the mapset; i.e. it uses the current mapset
(not the mapset= value), which may not even exist in the source
location.

This is quite clearly a bug; it should switch mapsets when it switches
locations. At present the code only uses the mapset= value when
opening the input map, overlooking the fact that it is also used
internally, e.g. when looking for the mask.

5

Hamish:

> side-suggestion: before running r.proj run v.in.region
> in the source location then v.proj from the taget location.

Jarekj wrote:

I always do as you suggest (in general other methods (with
cs2cs) in grass are to complicated),

sorry, I don't understand. is that to say that GRASS is easier or harder
than other methods? (ie, if harder, what areas need improving and what is
a clearer (proj4/gdal) method to use?)

Script must additionally check if no mask is present in PERMANENT or
it is a bug and will be removed in the feature

....

So as I understand the problem will heppen when I have mask in mapset in
source which is named as my target mapset in target location: ussually
it will concern PERMANENT?

* Maps are searched for in the current mapset search path. (g.mapsets)

* PERMANENT is, by default, in the map search path. (hence the name)

* If you like you can use g.mapsets to remove PERMANENT from the (source
  mapset's) search path, then later re-add it, but that may leave you
  with a damaged mapset if the script exits before the path is restored.

* You can switch between mapsets/locations within your script using the
  g.mapset (without an "s" at the end) module.

* When raster maps are read, the mapset search path is searched for a map
  called MASK. If a map by that name is found in the search path then it
  is applied as the map of interest is read from the disk.

* As Glynn noted, the module internally switches the current mapset, and
  so the for the time when the source mapset becomes the current mapset
  if there is a MASK map in the (source) search path, it will be used.
  (AFAIU)

work arounds include -

o temporarily switching into the source mapset with g.mapset (no "s"),
  removing everything but the current mapset from the source mapset's
  search path with g.mapsets (with an "s"), then switching back to the
  target mapset before running r.proj, then after r.proj is done going
  back and adding PERMANENT+any others removed back into the search path.
  As hinted above I think this solution to be rather poor.

o temporarily switching into the source mapset with g.mapset (no "s"),
  and then run "g.findfile element=cell file=MASK". The result (and exit
  code) will tell you if there is a map called MASK in the search path.
  What you do then is up to you. Maybe switch back to the original
  mapset/location and exit with an error or proceed with a warning.

o change the method used by the r.proj code. (as suggested by Glynn)

hope that makes the "why" a bit clearer.

Hamish

6

H:

work arounds include -

...

o change the method used by the r.proj code. (as suggested
by Glynn)

sorry, after re-reading Glynn's post that isn't what he said at all.

It is a bit muddy to me if the source map should observe any MASK. I think
it only makes sense to leave the question open for the source mapset MASK,
as a MASK in the target mapset is involved with writing not reading.

As for source mapset MASKs, in my personal opinion it is better to treat
the source map as a stand-alone external data source, and ignore any MASKs.
Sort of like the old r.out.gdal.sh script.

Hamish

7

Hamish pisze:

Hamish:
  

side-suggestion: before running r.proj run v.in.region
in the source location then v.proj from the taget location.
      
Jarekj wrote:
  
I always do as you suggest (in general other methods (with
cs2cs) in grass are to complicated),
    
sorry, I don't understand. is that to say that GRASS is easier or harder
than other methods? (ie, if harder, what areas need improving and what is
a clearer (proj4/gdal) method to use?)

In general in grass help is also suggrested method with manual recalculating extension of grass region with cs2cs (and there is also wrapper for that program), the method you suggested is also in help, both method are "grass" method

Many Thanks for information you added below, It will help me to finish my work faster
Jarek

  

Script must additionally check if no mask is present in PERMANENT or
it is a bug and will be removed in the feature
    

....
  

So as I understand the problem will heppen when I have mask in mapset in
source which is named as my target mapset in target location: ussually
it will concern PERMANENT?
    
* Maps are searched for in the current mapset search path. (g.mapsets)

* PERMANENT is, by default, in the map search path. (hence the name)

* If you like you can use g.mapsets to remove PERMANENT from the (source
  mapset's) search path, then later re-add it, but that may leave you
  with a damaged mapset if the script exits before the path is restored.

* You can switch between mapsets/locations within your script using the
  g.mapset (without an "s" at the end) module.

* When raster maps are read, the mapset search path is searched for a map
  called MASK. If a map by that name is found in the search path then it
  is applied as the map of interest is read from the disk.

* As Glynn noted, the module internally switches the current mapset, and
  so the for the time when the source mapset becomes the current mapset
  if there is a MASK map in the (source) search path, it will be used.
  (AFAIU)

work arounds include -

o temporarily switching into the source mapset with g.mapset (no "s"),
  removing everything but the current mapset from the source mapset's
  search path with g.mapsets (with an "s"), then switching back to the
  target mapset before running r.proj, then after r.proj is done going
  back and adding PERMANENT+any others removed back into the search path.
  As hinted above I think this solution to be rather poor.

o temporarily switching into the source mapset with g.mapset (no "s"),
  and then run "g.findfile element=cell file=MASK". The result (and exit
  code) will tell you if there is a map called MASK in the search path.
  What you do then is up to you. Maybe switch back to the original
  mapset/location and exit with an error or proceed with a warning.

o change the method used by the r.proj code. (as suggested by Glynn)

hope that makes the "why" a bit clearer.

Hamish

8

Hamish wrote:

> Script must additionally check if no mask is present in PERMANENT or
> it is a bug and will be removed in the feature
....
> So as I understand the problem will heppen when I have mask in mapset in
> source which is named as my target mapset in target location: ussually
> it will concern PERMANENT?

* Maps are searched for in the current mapset search path. (g.mapsets)

That reminds me of another issue: the mapset search path is determined
by the contents of the SEARCH_PATH file in the current directory, not
the environment, so it isn't affected by G__create_alt_env() and
G__switch_env(). Also, the contents are cached, so switching mapset or
location won't change any search path.

That isn't relevant to the MASK, though, as that is always obtained
from the current mapset (G_mapset()), not the mapset search path.

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

9

Jarekj wrote:

So as I understand the problem will heppen when I have mask in mapset in
source which is named as my target mapset in target location:

Correct.

If you try to reproject "srcmap in srcmapset in srclocation" to
"dstmap in dstmapset in dstlocation", the source map will be masked
with "MASK in dstmapset in srclocation", which is a bug.

In many cases, "dstmapset in srclocation" won't exist, but if
dstmapset is PERMANENT, it will.

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

10

On Mon, Jun 30, 2008 at 9:19 AM, Jarosław Jasiewicz <jarekj@amu.edu.pl> wrote:

In general in grass help is also suggrested method with manual recalculating
extension of grass region with cs2cs (and there is also wrapper for that
program), the method you suggested is also in help, both method are "grass"
method

My personal method is (if I understand correctly what you mean) if I dom't
know the corners in the target location:

1. go into source location, set region as I need
2. run there v.in.region to generate region box vector map

3. go into target location
4. run v.proj to reproject region box vector map from source location
5. run g.region vect=region_box res=value -pa
6. run r.proj to reproject desired raster map.

Like this I don't need cs2cs or m.proj, through the "box" trick I get the
corners easily reprojected.
(see also r.proj manual page for this hint)

Markus