[GRASS-user] Re-organizing Project

GRASS GRASS user
1

On Dec 22, 2009, at 1:57 PM, grass-user-request@lists.osgeo.org wrote:

Date: Tue, 22 Dec 2009 11:45:56 -0800 (PST)
From: Rich Shepard <rshepard@appl-ecosys.com>
Subject: [GRASS-user] Re-organizing Project
To: grass-users@lists.osgeo.org
Message-ID:
       <alpine.LNX.2.00.0912221132470.14561@salmo.appl-ecosys.com>
Content-Type: TEXT/PLAIN; format=flowed; charset=US-ASCII

  Since I'm having all sorts of issues getting usable data in one place so I
can start modeling, I should step back and reorganize what I'm doing. I've
just picked up that mapsets can be used for different issues or themes as
well as different users. This leads to several questions for which I don't
find answers in the Book or Web site; they may be in there but I've missed
them.

  For the current project, all data is in Oregon. All source data files have
the same projection information, but each has different bounds. I'm thinking
that I should put all source data in a single location, called 'Oregon.' The
location's PERMANENT mapset, DEFAULT_WIND and PROJ_INFO will come from the
state boundary outline (an .e00 file). That will be placed in the 'statebnd'
mapset and should have the same WIND as the DEFAULT_WIND in PERMANENT.
County outlines will be placed in the 'cntybnd' mapset which should have the
same WIND, too.

  Then I use v.proj and r.proj to copy each vector map and the raster DEM
map to their own mapsets within the Oregon location. The projection
information should all the same, but each will have different WIND
boundaries.

  If I'm correct on the above, then it should be simpler and correct to
create a project mapset in the same location since it, too, has the same
projection, but bounds the smallest area.

  With the above scheme, do I need to run g.region to set the bounds for
each mapset/map theme? I ask because some of the maps (e.g., stream/rivers,
hydrologic units, soil mapping units) may extend past the state or county
boundaries.

  Re-learning after 5-6 years of not using GRASS at all will get me back to
speed so I can complete this project.

Thanks,

Rich

Rich,

Maybe I can help with suggesting a work flow. To start with, I want to explain a couple features about regions and projecting. Some of this has been explained before, but maybe this can put it all in one place.

Regions:
Regions are a set of 4 extents (6 for volumes) and resolution within a given geographic projection (I include latlon and xy here). The idea of regions is to define and limit where many GIS computational operations take place. Most importantly, RASTER data are READ only from within the current region (this is the default situation, but sometimes can be overridden for some operations). VECTOR data may or may not be read from only within the current region (this is NOT the default, but can be set for some operations).

There are various ways to set regions, but the most direct and reliable is to use g.region. You can set each extent and the resolution independently, or set the region to match one or more maps. To set a region to match the extents and resolution of a raster map named "map1", you run g.region rast=map1. All subsequent raster operations are constrained to only read data from within the extents (i.e., the NSEW boundaries) of map1. Data will be read at the resolution of map1 too.

You can set the region extents to match the maximum of multiple maps by giving g.region a list: g.region rast=map1,map2,map3 will set the region to the maximum extents of the 3 maps combined. The same thing goes for vector maps: g.region vect=vmap1,vmap2,vmap3. However, I'm not sure what happens with g.region rast=map1,map2,map3 vect=vmap1,vmap2,vmap3. I think that rast= will override vect= or vice versa.

Projecting (r.proj and v.proj):
If RMapSrc is a raster map in the source location, r.proj will read the the cells of RMapSrc **that lie within the extents of the current region of the target location** and create a new map, RMapTgt, **at the resolution of the current region** whose cell values are calculated from a transformation function that maps information from a map in one projection into another. No data are read from the cells of RMapSrc that lie outside the current region of the target location. So if you want to reproject all of RMapSrc (you may or may not), it is important that the extents of the current region in the target location are set so that they extend beyond the *reprojected extents* of RMapSrc. Sometimes it is easy to set this and sometimes it is more difficulty. This is where the v.in.region trick comes in.

If VMapSrc is a vector map in a source location, v.proj will read ALL of the data in VMapSrc and create a new map in the target location such that the coordinates of each node/point in VMapSrc is transformed to a node/point in the new map. Resolution is irrelevant for a vector map. VMapSrc will be reprojected in its entirety regardless of the region settings in the target location.

So, if you create a vector map VregionSrc in the source location (using v.in.region) that is the same size and shape as the raster map you want to reproject, and if you reproject VregionSrc into the target location, you can then run g.region vect=VregionSrc and set the current region in the target location to exactly match the raster area you want to reproject. This way, ALL the data in RMapSrc will be read and used to create a new map in the target location.

Workflow:
-Start in the source location.
1. Set the region to match the map you want to reproject in the source location: g.region rast=RMapSrc.
   (if you have multiple maps to match, you can do g.region rast=map1,map2,map3 or you can even do
   g.region vect=vmap1,vmap2,vmap3 if this encompasses the raster area you want to reproject.)
2. Create VregionSrc as a rectangular vector area to match the extents of the region: v.in.region output=VregionSrc

-Change to target location (quit GRASS and restart in target).
1. Reproject VregionSrc: v.proj input=VregionSrc location=SourceLocation mapset=SouceMapset output=VregionSrc
2. Set the current region in the target location to match the reprojected VregionSrc in the target region and
   set the resolution as desired: g.region vect=VregionSrc res=resolution.
3. Reproject RMapSrc, creating a new map in the target location at the desired resolution, with cells that match up
   with the cells of RMapSrc: r.proj input=RMapSrc location=SourceLocation mapset=SouceMapset output=NewReprojectedMap.

I hope this helps

Michael

2

Michael,

This email should be appended to the r.proj man page! Or put on the
wiki! Great explanation...

Rich,

I'm not understanding why you need / want different mapsets. I always
see mapsets as a group of maps related to each other that I would like
to separate in groups to keep things organized. For instance, I use
different mapsets for different climate simulation scenarios. You
could use different mapsets for different satellite images (years,
sensors etc). But keep in mind that you might not need several mapsets
and centanly, one mapset for each map is an overkill. Now, the
PERMANENT mapset is a special one that can be accessed from within any
other mapset. So that's the place to put maps like boundaries,
counties, roads etc, that will be accessed by all.

Another way of organizing is by creating different mapsets for
different models and in PERMANENT you can put all common data to the
models (DEM, streams, etc)

This page helped me a lot in understanding what mapsets are for:
http://grass.osgeo.org/grass64/manuals/html64_user/helptext.html

cheers
daniel

On Tue, Dec 22, 2009 at 9:51 PM, Michael Barton <michael.barton@asu.edu> wrote:

On Dec 22, 2009, at 1:57 PM, grass-user-request@lists.osgeo.org wrote:

Date: Tue, 22 Dec 2009 11:45:56 -0800 (PST)
From: Rich Shepard <rshepard@appl-ecosys.com>
Subject: [GRASS-user] Re-organizing Project
To: grass-users@lists.osgeo.org
Message-ID:
<alpine.LNX.2.00.0912221132470.14561@salmo.appl-ecosys.com>
Content-Type: TEXT/PLAIN; format=flowed; charset=US-ASCII

Since I'm having all sorts of issues getting usable data in one place so I
can start modeling, I should step back and reorganize what I'm doing. I've
just picked up that mapsets can be used for different issues or themes as
well as different users. This leads to several questions for which I don't
find answers in the Book or Web site; they may be in there but I've missed
them.

For the current project, all data is in Oregon. All source data files have
the same projection information, but each has different bounds. I'm thinking
that I should put all source data in a single location, called 'Oregon.' The
location's PERMANENT mapset, DEFAULT_WIND and PROJ_INFO will come from the
state boundary outline (an .e00 file). That will be placed in the 'statebnd'
mapset and should have the same WIND as the DEFAULT_WIND in PERMANENT.
County outlines will be placed in the 'cntybnd' mapset which should have the
same WIND, too.

Then I use v.proj and r.proj to copy each vector map and the raster DEM
map to their own mapsets within the Oregon location. The projection
information should all the same, but each will have different WIND
boundaries.

If I'm correct on the above, then it should be simpler and correct to
create a project mapset in the same location since it, too, has the same
projection, but bounds the smallest area.

With the above scheme, do I need to run g.region to set the bounds for
each mapset/map theme? I ask because some of the maps (e.g., stream/rivers,
hydrologic units, soil mapping units) may extend past the state or county
boundaries.

Re-learning after 5-6 years of not using GRASS at all will get me back to
speed so I can complete this project.

Thanks,

Rich

Rich,

Maybe I can help with suggesting a work flow. To start with, I want to explain a couple features about regions and projecting. Some of this has been explained before, but maybe this can put it all in one place.

Regions:
Regions are a set of 4 extents (6 for volumes) and resolution within a given geographic projection (I include latlon and xy here). The idea of regions is to define and limit where many GIS computational operations take place. Most importantly, RASTER data are READ only from within the current region (this is the default situation, but sometimes can be overridden for some operations). VECTOR data may or may not be read from only within the current region (this is NOT the default, but can be set for some operations).

There are various ways to set regions, but the most direct and reliable is to use g.region. You can set each extent and the resolution independently, or set the region to match one or more maps. To set a region to match the extents and resolution of a raster map named "map1", you run g.region rast=map1. All subsequent raster operations are constrained to only read data from within the extents (i.e., the NSEW boundaries) of map1. Data will be read at the resolution of map1 too.

You can set the region extents to match the maximum of multiple maps by giving g.region a list: g.region rast=map1,map2,map3 will set the region to the maximum extents of the 3 maps combined. The same thing goes for vector maps: g.region vect=vmap1,vmap2,vmap3. However, I'm not sure what happens with g.region rast=map1,map2,map3 vect=vmap1,vmap2,vmap3. I think that rast= will override vect= or vice versa.

Projecting (r.proj and v.proj):
If RMapSrc is a raster map in the source location, r.proj will read the the cells of RMapSrc **that lie within the extents of the current region of the target location** and create a new map, RMapTgt, **at the resolution of the current region** whose cell values are calculated from a transformation function that maps information from a map in one projection into another. No data are read from the cells of RMapSrc that lie outside the current region of the target location. So if you want to reproject all of RMapSrc (you may or may not), it is important that the extents of the current region in the target location are set so that they extend beyond the *reprojected extents* of RMapSrc. Sometimes it is easy to set this and sometimes it is more difficulty. This is where the v.in.region trick comes in.

If VMapSrc is a vector map in a source location, v.proj will read ALL of the data in VMapSrc and create a new map in the target location such that the coordinates of each node/point in VMapSrc is transformed to a node/point in the new map. Resolution is irrelevant for a vector map. VMapSrc will be reprojected in its entirety regardless of the region settings in the target location.

So, if you create a vector map VregionSrc in the source location (using v.in.region) that is the same size and shape as the raster map you want to reproject, and if you reproject VregionSrc into the target location, you can then run g.region vect=VregionSrc and set the current region in the target location to exactly match the raster area you want to reproject. This way, ALL the data in RMapSrc will be read and used to create a new map in the target location.

Workflow:
-Start in the source location.
1. Set the region to match the map you want to reproject in the source location: g.region rast=RMapSrc.
(if you have multiple maps to match, you can do g.region rast=map1,map2,map3 or you can even do
g.region vect=vmap1,vmap2,vmap3 if this encompasses the raster area you want to reproject.)
2. Create VregionSrc as a rectangular vector area to match the extents of the region: v.in.region output=VregionSrc

-Change to target location (quit GRASS and restart in target).
1. Reproject VregionSrc: v.proj input=VregionSrc location=SourceLocation mapset=SouceMapset output=VregionSrc
2. Set the current region in the target location to match the reprojected VregionSrc in the target region and
set the resolution as desired: g.region vect=VregionSrc res=resolution.
3. Reproject RMapSrc, creating a new map in the target location at the desired resolution, with cells that match up
with the cells of RMapSrc: r.proj input=RMapSrc location=SourceLocation mapset=SouceMapset output=NewReprojectedMap.

I hope this helps

Michael

_______________________________________________
grass-user mailing list
grass-user@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/grass-user

3

On Tue, 22 Dec 2009, Michael Barton wrote:

Maybe I can help with suggesting a work flow. To start with, I want to
explain a couple features about regions and projecting. Some of this has
been explained before, but maybe this can put it all in one place.

Michael,

   Thank you. I actually understand most of this (other than multiple maps
for setting the largest region), but did not clearly express myself. The key
question for me is why it's not working as it should.

There are various ways to set regions, but the most direct and reliable is
to use g.region.

   Without getting into my issues of reorganizing everything so it's in a
single location, with each map and/or issue a separate mapset, the largest
region should be the entire state. And, the state boundary map and counties
map _should_ have the same region even though the west boundary of the state
is doubled (from the .e00 import). I guess that's the effects of erosion by
the Pacific Ocean on the coast. It's not worth my time to clean that now.

   All the vector maps can happily reproject into the collective location.
It's the raster DEM that's being recalcitrant.

If RMapSrc is a raster map in the source location, r.proj will read the
the cells of RMapSrc **that lie within the extents of the current region
of the target location** and create a new map, RMapTgt, **at the
resolution of the current region** whose cell values are calculated from a
transformation function that maps information from a map in one projection
into another.

   Which suggests that I should have had no problems. The DEM coverage is for
the northwestern portion of the state; approximately from the top of the
Cascade Mountains on the east to the Pacific Ocean on the west, the Columbia
River on the north, and someplace about the southern end of the Willamette
Valley on that end. This source region is a sub-set of the target region
(the entire state), and any cells in the water can happily be ignored and
clipped off because I'm not doing bathymetry on this project.

   The region of the raster source map is a subset of the target region.
That's why I don't understand why GRASS kept throwing errors at me.

Workflow:
-Start in the source location.
1. Set the region to match the map you want to reproject in the source location: g.region rast=RMapSrc.

   That's the way I set the source region after importing the raster file.

2. Create VregionSrc as a rectangular vector area to match the extents of
the region: v.in.region output=VregionSrc

   I did this, too. It provides a box that is a subset of the target region.

-Change to target location (quit GRASS and restart in target).
1. Reproject VregionSrc: v.proj input=VregionSrc location=SourceLocation
mapset=SouceMapset output=VregionSrc

   I _think_ this step worked for me. I'll have to do it over again to be
sure.

2. Set the current region in the target location to match the reprojected
VregionSrc in the target region and set the resolution as desired:
g.region vect=VregionSrc res=resolution.

   I may not have done this.

3. Reproject RMapSrc, creating a new map in the target location at the
desired resolution, with cells that match up with the cells of RMapSrc:
r.proj input=RMapSrc location=SourceLocation mapset=SouceMapset
output=NewReprojectedMap.

   It's all quite different from when I worked with versions 4.1 through 4.3.
Heck, even the modules and steps in the 2nd edition of the GRASS book are
frequently no longer valid.

Thanks very much,

Rich

4

On Tue, 22 Dec 2009, Daniel Victoria wrote:

I'm not understanding why you need / want different mapsets. I always see
mapsets as a group of maps related to each other that I would like to
separate in groups to keep things organized. For instance, I use different
mapsets for different climate simulation scenarios. You could use
different mapsets for different satellite images (years, sensors etc). But
keep in mind that you might not need several mapsets and certainly, one
mapset for each map is an overkill. Now, the PERMANENT mapset is a special
one that can be accessed from within any other mapset. So that's the place
to put maps like boundaries, counties, roads etc, that will be accessed by
all.

Another way of organizing is by creating different mapsets for different
models and in PERMANENT you can put all common data to the models (DEM,
streams, etc)

This page helped me a lot in understanding what mapsets are for:
http://grass.osgeo.org/grass64/manuals/html64_user/helptext.html

   Yes, the latter URL repeated the information in the book. But, there seems
to be many different preferences for organizing data. I'm open to all
opinions on this.

   I have source data for the following:
   DEM 10m, northwest portion of the state (raster); projection: Oregon LCC
   streams/rivers in the the state at 24K; projection: Oregon LCC
   lakes/reservoirs in the state at 24K; projection: Oregon LCC
   dams in the state at 24K; projection: Oregon LCC
   hydrologic units (basins) in the state at 24K; projection: Oregon LCC
   soils (map units, line and point features) at 24K in one county;
     projection: Oregon LCC
   roads/highways in the state at 24K; projection: Oregon LCC
   state and county boundaries (not needed for this project, but to
     establish largest region if necessary)
   key point features specific to this project; projection: Oregon LCC

   Originally, each map was in a separate location. Now I want to put them
all in a single, statewide location. I don't care if they're all
imported/re-projected into the PERMANENT map set or individual ones.

   There's also the project-specific location which is a small subset of one
county. After getting all the source data correctly in their base
location/mapset(s), I want to define the region of this location or mapset
as slightly larger than the drainage basin I need to model. I plan to define
this by viewing the stream of interest with the HUC basins overlaid.
Running r.watershed on the DEM should match this area closely (I hope). This
will be the project's working directory. I was hoping to be ready to model
by now, but am high-centered on both organizing all source data and
moving/re-projecting it to a common place.

   There will be several (or more) projects in this and other states where I
will need to run GRASS models so defining an organizational structure and
being able to correctly move maps around is vital. I would like to establish
a standard now, use it for this project, and also use it for future
projects.

   Everyone's free to suggest organization. Interestingly, both the book (2nd
ed) and above URL focus more on mapsets per user, but show they can be used
per map theme.

Thanks,

Rich

5

Glad this is helpful. If there is a wiki place for it, feel free to post it.

Michael
____________________
C. Michael Barton
Director, Center for Social Dynamics & Complexity
Professor of Anthropology, School of Human Evolution & Social Change
Arizona State University

Phone: 480-965-6262
Fax: 480-965-7671
www: www.public.asu.edu/~cmbarton, http://csdc.asu.edu

On Dec 22, 2009, at 5:54 PM, Daniel Victoria wrote:

Michael,

This email should be appended to the r.proj man page! Or put on the
wiki! Great explanation...

6

Rich wrote:

  It's all quite different from when I worked with versions 4.1 through 4.3.
Heck, even the modules and steps in the 2nd edition of the GRASS book are
frequently no longer valid.

I am not sure what version the 2nd edition book is aimed at, but fwiw
we have maintained backwards compatibility with all of GRASS 6.x for the
last 5 years. all modules should still be there (at least 99% of them
anyway). And I'd say that most of the grass 5 modules are still there too.
A rewrite of the vector modules was the massive change between 5 and 6.

the gui and menu positions have changed a lot since then, and the d.*
modules don't work with the new GUI.

Hamish

7

On Tue, 22 Dec 2009, Daniel Victoria wrote:

I'm not understanding why you need / want different mapsets. I always see
mapsets as a group of maps related to each other that I would like to
separate in groups to keep things organized.

Another way of organizing is by creating different mapsets for different
models and in PERMANENT you can put all common data to the models (DEM,
streams, etc)

Daniel,

   This is yet another way of organizing data. Thinking about it in light of
what I want to do, it makes sense. I can put all the basic maps in
PERMANENT, then have a mapset for the project-specific area with its more
restricted bounds.

   I'll re-arrange everything once again. :slight_smile:

Rich