[GRASS-user] Working with USA Public Land Survey System

[Hamish]

> Henry House wrote:
> > I have a set of surveyed-section boundaries* for northern California
> > from the United States Public Land Survey System and would like to
> > use GRASS to find geometric fractions of these roughly square
> > sections. For example, I would like to divide a particular section
> > into equal quarters by bisecting it east-west and north-south and
> > then save the northwest quarter for further analysis. Can anyone
> > suggest a method to do this acurately?
> >
> > * Used for legal descriptions of parcels given in terms of townships
> > and ranges.

Dylan Beaudette wrote:
> I think that the approach suggested by Hamish will work, but only when
> using a  projection where the PLSS sections area oriented _exactly_ NS
> and EW. For  example, looking at the PLSS data in the AEA projection
> that we use for many  things, the sections are _not_ oriented as such:
> 
> http://169.237.35.250/~dylan/temp/plss1.png
> 
> However, there is an 'angle' operator in the v.mkgrid module, which
> may be  able to account for this- although I cannot think of an
> _automated_ approach  to deriving this angle.
> 
> I have extracted a single section, marked red in the above image, like
> this:
> 
> v.extract in=pls out=pls1 list=80800
> 
> where 80800 is the cat, or feature id, of the section of interest.
> Zooming to  the extents of this section:
> 
> g.region vect=pls1
> 
> http://169.237.35.250/~dylan/temp/plss2.png
> 
> looking at the options for v.mkgrid:
> 
> Parameters:
>        map   name of vector map
>       grid   number of ROWS and COLUMNS in grid
>   position   Where to place the grid:
>         region - current region
>         coor - use 'coor' and 'box' options
>              options: region,coor
>              default: region
>       coor   lower left EASTING and NORTHING coordinates of map
>        box   WIDTH and HEIGHT of boxes in grid
>      angle   angle of rotation (in degrees counter-clockwise)
>              default: 0
> 
> .... it looks like we might be able to define the grid in terms of a 
> lower-left point, and the width and height of grid cells, and an angle
> of  rotation.
> 
> this will require some math, but might easily be accomplished with
> some right  triangle-based trig.
> 
> fitting the region to _exactly_ that of the bounding box of our
> extracted  section: 
>  g.region `v.info -g pls1  | awk ' {gsub("north","n");
>  gsub("south","s") ; 
> gsub("east","e") ; gsub("west","w"); gsub("top","t") ;
> gsub("bottom","b")}  $1 !~ /b|t/ {print $0}' | tr "\n" " "` -p
> 
> extracting vertices of selected section: can do this with v.out.ascii
> and awk,  or with v.to.points


If true N-S E-W grid is wanted, why not just run v.mkgrid in a lat-lon
location and reproject that to the target location with v.proj?

IIRC, v.mkgrid breaks the line ~ 3 times per box, so results can look
quite nice even when "curved". I think this is how Markus made this
screenshot?

http://grass.ibiblio.org/grass60/screenshots/images/etopo2_grass_laea_9_48N_0E.jpg



Hamish