LandScript - Tutorial

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20

# Script to convert binary image file into a UTM LandSerf file.
# v1.1 Jo Wood 2nd October, 2007.
# -------------------------------------------------------------
version(1.0);

# Store the base directory in a variable to allow easy changing.
baseDir = "/data/us/lincoln/landscriptTutorial/";

# Open the .bil file and store the raster as a variable.
dem = open(baseDir & "lincolnNED.bil");

# Add the projection metadata to the DEM.
edit(dem,"projection","LatLong");
edit(dem,"ellipsoid","wgs 84");

# Create a new Universal Transverse Mercator (UTM) reprojected DEM.
dem = reproject(dem_,"UTM","true",30,30);

# Save the edited raster as a LandSerf file.
save(dem, baseDir & "step1.srf");

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29

# Script to edit the metadata of a spatial object.
# v1.1 Jo Wood 2nd October, 2007.
# -------------------------------------------------------------
version(1.0);

# Store the base directory in a variable to allow easy changing.
baseDir = "/data/us/lincoln/landscriptTutorial/";

# Open the file from the previous step and store as a variable.
dem = open(baseDir & "step1.srf");

# Give the raster a new colour table.
colouredit(dem,"land3");

# Guess at the water level by finding modal elevation
maxFreq = info(dem,"modef");
if (maxFreq > 1000);
{
  # If there are more than 1000 cells of same height, colour them blue
  waterLevel = info(dem,"mode");
  colouredit(dem,"addRule",waterLevel & " 141 141 166 (D)");
}

# Update title and notes.
edit(dem,"title","Lincoln 30m DEM");
edit(dem,"notes","Elevation data from USGS National Elevation Dataset (NED) available from http://seamless.usgs.gov");

# Save the edited raster as a LandSerf file.
save(dem, baseDir & "step2.srf");

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35

# Script to process a landcover raster.
# v1.2 Jo Wood 2nd October, 2007.
# -------------------------------------------------------------
version(1.0);

# Store the base directory in a variable to allow easy changing.
baseDir = "/data/us/lincoln/landscriptTutorial/";

# Open the landcover file, reproject it to UTM.
landcover = open(baseDir & "lincolnLandcover.bil");
edit(landcover,"projection","LatLong");
edit(landcover,"ellipsoid","wgs 84");
landcover = reproject(landcover_,"UTM","false",30,30);

# Update the landcover metadata.
edit(landcover,"title","Lincoln landcover");
colouredit(landcover,"file",baseDir & "landcover.ctb");
edit(landcover,"attributes",baseDir & "landcover.atr");
edit(landcover,"attCol",3);
edit(landcover,"type","other");

# Find average height of all forested areas.
dem = open(baseDir & "step2.srf");
forestHeight = new(landcover);
forestHeight = ifelse((landcover >= 40) and (landcover <=49),dem,null());
echo("Average height of forest is " & round(info(forestHeight,"mean")) & "m.");

# Find average height of scrubland areas.
dem = open(baseDir & "step2.srf");
scrubHeight = new(landcover);
scrubHeight = ifelse(landcover == 51,dem,null());
echo("Average height of scrubland is " & round(info(scrubHeight,"mean")) & "m.");

# Save the edited landcover as a LandSerf file.
save(landcover, baseDir & "step3.srf");

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49

# Script to perform some simple morphometric analysis on a DEM.
# v1.1 Jo Wood 2nd October, 2007.
# -------------------------------------------------------------
version(1.0);

# Store the base directory in a variable to allow easy changing.
baseDir = "/data/us/lincoln/landscriptTutorial/";

# Open the DEM and landcover raster from the previous steps.
dem = open(baseDir & "step2.srf");
landcover = open(baseDir & "step3.srf");

# Create a shaded relief image saved as a file.
edit(dem,"vExag",3);
draw(baseDir & "relief.jpg",dem,landcover,"null","relief");

# Output some relevant DEM metadata
showMetadata(dem,"elevation");

# Calculate and output slope of DEM
slope = surfparam(dem,"slope");
draw(baseDir & "slope.jpg",dem,slope,"null","relief");
showMetadata(slope,"slope");

# Create images of the other surface parameters.
draw(baseDir & "aspect.jpg",dem,surfparam(dem,"aspect"),"null","relief");
draw(baseDir & "profc.jpg", dem,surfparam(dem,"profc"),"null","relief");
draw(baseDir & "planc.jpg", dem,surfparam(dem,"planc"),"null","relief");

# Finally create a surface feature map and save as a LandSerf file.
features = surfparam(dem,"feature");
save(features, baseDir & "LincolnFeatures.srf");


# Function to display some spatial metadata.
# Params: raster Raster from which to extract metadata.
#         name Name of raster.
function showMetadata(raster,name)
{
  echo();
  echo(name);
  echo("East-west extent: "& (info(raster,"E")-info(raster,"W"))/1000 & "km.");
  echo("North-south extent: "& (info(raster,"N")-info(raster,"S"))/1000 & "km.");
  echo(name&" range from "& info(raster,"min")&" to "&info(raster,"max")&".");
  echo("Mean "&name&": "& info(raster,"mean"));
  echo("Modal "&name&" value: "& info(raster,"mode"));
  echo("Standard deviation of "&name&": "& info(raster,"stdev"));
  echo("Skewness of "&name&": "& info(raster,"skew"));
}

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28

# Script to perform more advanced morphometric analysis on a DEM.
# v1.1 Jo Wood 2nd October, 2007.
# -------------------------------------------------------------
version(1.0);

# Store the base directory in a variable to allow easy changing.
baseDir = "/Users/jwo/data/us/lincoln/landscriptTutorial/";

# Open the DEM from the earlier step.
dem = open(baseDir & "step2.srf");

# Calculate 11x11 cell smothed elevation and draw output.
winSize = 11*info(dem,"xRes");
echo("Smoothing elevation using "&winSize&"m x "&winSize&"m window");
draw(baseDir & "elev11.jpg", surfparam(dem,"elev",11),"null","null","relief");

# Calculate 65x65 cell smothed elevation and draw and save output.
winSize = 65*info(dem,"xRes")/1000;
echo("Smoothing elevation using "&winSize&"km x "&winSize&"km window");
elev65 = surfparam(dem,"elev",65);
draw(baseDir & "elev65.jpg",elev65,"null","null","relief");
save(elev65,baseDir&"elev65.srf");

# Create a 5x5 feature network and metric surface network
feat5 = surfparam(dem,"netfeature",5,1);
msn = surfnetwork(dem,feat5);
draw(baseDir & "featnetwork.jpg",dem,"null",msn,"relief");
save(msn,baseDir&"msn.vec");