[GRASS-user] Fixing Collared DRGs

[Tom Russo]

On Mon, Mar 31, 2008 at 03:00:41PM -0700, we recorded a bogon-computron collision of the <dhirsch226 at gmail.com> flavor, containing:
> Hi folks- I'm new to GRASS, but am eager to learn a GIS that doesn't require
> me to use Windows.  I can manage to do the basics with GRASS (although I'm
> having trouble with the python GUI, but that's another story).  What I'm
> trying to do now is remove white collars from some GeoTiff topo maps I got
> for an area in Montana (my region of interest overlaps the corners of four
> quadrangles).  These have white borders, and the border region is the same
> white as the white of the map itself (color index 1).  When I downloaded the
> maps, there was a link to some scripts (for ArcMap) to remove the collars,
> and I imagine that there is some equivalent way to do it for GRASS.
> 
> I'm running GRASS 6.3.0RC5 (using the X11 tcltk interface) on a PowerMac G4.
> Thanks,

It is a little tedious but not especially difficult, and easily scripted.  
Here's how I do it:

Go to www.geocomm.com and download their 24kgrid.shp "quad index" file.  It
contains the lat/lon rectangles for all the USGS 7.5' DRGs in the whole US.

Extract the set of outlines for your state to a shapefile:

  ogr2ogr -where "STATE='MT'" montana_quads.shp 24kgrid.shp

Now you need to translate those outlines to the UTM zone and datum that your
maps are in.  Some states use quad maps all in NAD27, some have a mixture of
NAD83 and NAD27 maps.  The 24kgrid shapefile has no associated .prj file,
and just contains the lat/lon boundaries of the quads (inside the neatlines).

If your GRASS location is UTM zone 12, NAD27 (and your DRGs are also NAD27, 
UTM zone 12), you could do this:

  ogr2ogr -t_srs EPSG:26712 -s_srs EPSG:4267 montana_quads_UTM12_NAD27.shp montana_quads.shp

This reprojects the (assumed) NAD27 lat/lon rectangles into UTM Zone 12/NAD27 
polygons.  These will be used to clip the collars off of DRGs.

Then import this shapefile into your GRASS location with v.in.ogr:

  v.in.ogr dsn=montana_quads_UTM12_NAD27.shp output=quads27

Once you've done that, you're ready to collar-clip all the NAD27 DRGs you want. 
Say you have a quad for Abbot Lake, MT that's called o48112d2.tif.  You could 
then do this:

  r.in.gdal input=o48112d2.tif output=o48112d2_withcollar
  v.extract where="MRC='48112-D2'" input=quads27 output=o48112d2_bounds
     (this selects only the single quad outline applicable)
  g.region rast=o48112d2_withcollar   # this just to set resolution properly
  g.region vect=o48112d2_bounds       # to set the region only to neatline
  v.to.rast input=o48112d2_bounds output=o48112d2_bounds use=val value=1
     (this one generates a raster that's 1 inside the neatline, null outside)
  echo "o48112d2_withoutcollar=if(!isnull(o48112d2_bounds),o48112d2_withcollar,null())" | r.mapcalc
    (this populates the _withoutcollar map only where the mask is non-null)
    (you could also do it with 
        r.mask input=o48112d2_bounds
        echo "o48112d2_withoutcollar=o48112d2_withcollar"| r.mapcalc
        r.mask -r
    )
  r.colors map=o48112d2_withoutcollar rast=o48112d2_withcollar
    (to copy colormap from original)

If the location is in a different datum than the DRG, I add a couple steps.  For
example, if my location is in NAD27 and the DRG is in NAD83, I'd do this:

  ogr2ogr -t_srs EPSG:26712 -s_srs EPSG:4269 montana_quads_UTM12_NAD83.shp montana_quads.shp

this creates a NAD27/UTM Zone 12 shapefile of the outlines of NAD83 quad maps.
I'd then import this as, say, "quads83" with v.in.ogr, and the steps would be:

  gdalwarp -t_srs EPSG:26712 o48112d2.tif o48112d2_NAD27.tif 
  r.in.gdal input=o48112d2_NAD27.tif output=o48112d2_withcollar
  v.extract where="MRC='48112-D2'" input=quads83 output=o48112d2_bounds
     (this selects only the single quad outline applicable)
  g.region rast=o48112d2_withcollar   # this just to set resolution properly
  g.region vect=o48112d2_bounds       # to set the region only to neatline
  v.to.rast input=o48112d2_bounds output=o48112d2_bounds use=val value=1
     (this one generates a raster that's 1 inside the neatline, null outside)
  echo "o48112d2_withoutcollar=if(!isnull(o48112d2_bounds),o48112d2_withcollar,n
ull())" | r.mapcalc
    (this populates the _withoutcollar map only where the mask is non-null)
    (you could also do it with
        r.mask input=o48112d2_bounds
        echo "o48112d2_withoutcollar=o48112d2_withcollar"| r.mapcalc
        r.mask -r
    )
  r.colors map=o48112d2_withoutcollar rast=o48112d2_withcollar
    (to copy colormap from original)

Only the first three steps are different.

-- 
Tom Russo    KM5VY   SAR502   DM64ux          http://www.swcp.com/~russo/
Tijeras, NM  QRPL#1592 K2#398  SOC#236 AHTB#1 http://kevan.org/brain.cgi?DDTNM
"And, isn't sanity really just a one-trick pony anyway? I mean all you get is
 one trick, rational thinking, but when you're good and crazy, oooh, oooh,
 oooh, the sky is the limit!"  --- The Tick