[GRASS-user] Re-organizing Project

[Michael Barton]

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

> Date: Tue, 22 Dec 2009 11:45:56 -0800 (PST)
> From: Rich Shepard <rshepard at appl-ecosys.com>
> Subject: [GRASS-user] Re-organizing Project
> To: grass-users at lists.osgeo.org
> Message-ID:
>        <alpine.LNX.2.00.0912221132470.14561 at 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