wiki:grass
Last modified 4 years ago Last modified on 06/18/13 12:48:13

Creating a GRASS location

You only need to set this up once. Next time, you can just run grass65 to start

grass65 --text

LOCATION:  pastp
MAPSET:    PERMANENT

DATABASE   /scratch/usrname

<ESC> <ENTER> 

Would you like to create location <pastp> ? (y/n) [y] y

Do you have all this information? (y/n) [y] y

Please specify the coordinate system for location <pastp>

A   x,y
B   Latitude-Longitude
C   UTM
D   Other Projection
RETURN to cancel

> D

Other Projection coordinate system? (y/n) [y] y
Please enter a one line description for location <pastp>

> 

=====================================================

=====================================================
ok? (y/n) [n] y

Please specify projection name
Enter 'list' for the list of available projections
Hit RETURN to cancel request
> stp
Do you wish to specify a geodetic datum for this location?(y/n) [y] y

Please specify datum name
Enter 'list' for the list of available datums
or 'custom' if you wish to enter custom parameters
Hit RETURN to cancel request
>nad83

Now select Datum Transformation Parameters
Please think carefully about the area covered by your data
and the accuracy you require before making your selection.

Enter 'list' to see the list of available Parameter sets
Enter the corresponding number, or <RETURN> to cancel request
>

Specify State FIPS (numeric) code
Enter 'list' for the list of states with corresponding FIPS codes
Hit RETURN to cancel request
>42

You have chosen state PA, Correct(y/n) [y] y


Specify County FIPS (numeric) code for state PA
Enter 'list' for the list of counties in PA with corresponding FIPS codes
Hit RETURN to cancel request
>45
You have chosen DELAWARE county, correct(y/n) [y] y

Specify State Plane 1927 or 1983
Enter '27' or '83'
Hit RETURN to cancel request
>83
Specify the correct units to use:
Enter the corresponding number
1.	US Survey Foot (Default for State Plane 1927)
2.	International Foot
3.	Meter
>1

DEFINE THE DEFAULT REGION

Just hit <ESC><ENTER>
Do you accept this region? (y/n) [y] > y
LOCATION <pastp> created!

Hit RETURN -->

<ESC><ENTER>

GRASS 6.5.svn (pastp):~ > r.in.gdal /scratch/adanner/public/delco20.gtiff output=delco20 -o
WARNING: Datum <unknown> not recognized by GRASS and no parameters found
WARNING: Over-riding projection check

r.in.gdal complete. Raster map <delco20> created.
GRASS 6.5.svn (pastp):~ > g.region rast=delco20
GRASS 6.5.svn (pastp):~ > d.mon x0
using default visual which is TrueColor
ncolors: 16777216
Graphics driver [x0] started
GRASS 6.5.svn (pastp):~ > d.rast delco20
 100%

Reducing resolution/Export

g.region rast=delco20 res=100
r.mapcalc delco100=delco20
r.info delco100

 |   Rows:         1102                                                       |
 |   Columns:      1102                                                       |
 |   Total Cells:  1214404                                                    |
 |        Projection: State Plane                                             |
 |            N:     280000    S:     169800   Res:   100                     |
 |            E:    2680200    W:    2570000   Res:   100                     |
 |   Range of data:    min = -140.4787  max = 498.4281                        |

cd /scratch/adas/
r.out.gdal delco100 output=delco100.tif nodata=-9999.0

Flooding the DEM

cd /scratch/username/
mkdir tmp
ln -s /scratch/adanner/public/terrastream/

./terrastream/Hydrology-flooding --input-elev=delco100.tif --input-elev-type=gdal --flooded-elev=delco100flood.tif --flooded-elev-type=gdal --tmppath=/scratch/adas/tmp


Import back into grass

r.in.gdal delco100flood.tif output=delco100flood -o
g.region rast=delco100
r.mapcalc diff=delco100flood-delco100
r.colors diff color=wave
d.mon x0
d.rast diff val=1-500

Flow Directions

terrastream/Hydrology-flow-directions --input-elev=delco100flood.tif --input-elev-type=gdal --tmppath=/scratch/adas/tmp/ --method=sfd --output-dirs=delco100dirs.tif --output-dirs-type=gdal

r.in.gdal delco100dirs.tif output=delcodirs -o
r.mapcalc dirs="if(delcodirs==0,null(),delcodirs)"
d.mon x0
d.rast dirs

Flow Accumulation

/terrastream/Hydrology-flow-accumulation --input-dir-edges=delco100flood.tds --output-flow=delco100flow.tif --output-flow-type=gdal --output-flow-gdal-driver=GTiff --tmppath=/scratch/adas/tmp/
r.in.gdal delco100flow.tif output=flow -o
d.rast flow  val=1000-1000000
r.mapcalc dots="if(flow>1000,1,null())"
r.to.vect dots output=rivers
d.erase
d.rast delco100
d.vect rivers color=black