Geologic Application of RadarSat S2 Mode Data in Northern Nova Scotia

The following image is a photo representation of a larger poster that I made along with Blair Sangster at COGS in March of 1999.

The project also included a detailed paper and presentation that was presented at the Center of Geographic Sciences Auditorium, in Lawrencetown Nova Scotia.

The project used RADARSAT S2 beam mode SAR imagery and ERDAS Imagine together to provide 3D models representing various geological terrain features.

PCI EASI Script [PIX_to_SHP_batch_export.eas]

This simple EASI script used with the MODEL command in PCI will batch convert vector files stored in PCIDSK (PIX) format into Shape File (SHP) format with the exact same file name as the input files. This was originally created for PCI Geomatica v9.1 but was last tested and working with no problems in PCI Geomatica v10.0


!—————————————————————————-
!—————————————————————————-
!
! Batch Export PIX Vector to SHP Script
! [PIX_to_SHP_batch_export.eas]
!
! This script will export vector segments from PIX files located in a
! given directory into Shape files and place the new files into the same
! directory. The script assumes that no SHP files with the same names
! already exist and all the input files are setup the same with the
! vector segments stored the same.
!
!—————————————————————————-
!—————————————————————————-
! Define variables
!—————————————————————————-

!to store location of input & output files

local string in_files

!for the file format and extension types

local string type, ext

!file names

local string bn, fn

!to store directory listing of input files

local mstring dirlist

!to store vector segment number

local integer vec

local integer i

local $Z

!—————————————————————————-
! Clear the EASI window and then show the header information
!—————————————————————————-

PRINT @(1 ,1,CLREOS)

print “———————————————————————–”
print @reverse,” ‘Batch export PIX vector to SHP’ EASI Script “,@alloff
print “”
print “”
print “”
print “———————————————————————–”
print “”
print ” This script will export vectors from PIX files in a given directory ”
print ” into SHP format using the same file names as the input files. ”
print “”
print ” All input vector files are expected to be setup the same with all ”
print ” vectors to be exported, stored using the same segment number.”
print “”
print “———————————————————————–”
print “”

!—————————————————————————-
! Collect input from user
!—————————————————————————-

print “Enter the directory that contains the PIX files to export to SHP:”

input “>” in_files

print “Enter the vector segment number that the vector is stored in:”

input “>” vec

!——————————————————————-
! Get the contents of the directory
!——————————————————————-

type =”pix

dirlist = getdirectory(in_files)

let $Z = “\

for i = 1 to f$len(dirlist)

!——————————————————————-
! Extract parts of the filenames
!——————————————————————-

fn = in_files + $Z + dirlist[i]

ext = getfileextension(fn)

bn = getfilebasename(fn)

if (ext ~= type) then

print “”
print “Exporting:”, bn, “from PIX to SHP”
print “”

!—————————————————————-
! Set up the parameters and execute the FEXPORT command
!—————————————————————-

fili = in_files + $Z + dirlist[i]
filo = in_files + $Z + bn
dbiw =
dbic =
dbib =
dbvs = vec
dblut =
dbpct =
ftype =”SHP
foptions =

R Fexport

PRINT @(1 ,1,CLREOS)

endif

endfor

PRINT @(1 ,1,CLREOS)

print “——————————————————————-”
print “”
print @reverse,” ‘Batch export PIX vector to SHP’ EASI Script Finished “,@alloff
print “”
print “——————————————————————-”

return

!—————————————————————————-
!—————————————————————————-

Multiple File Subset Script – PCI EASI Script

This simple EASI script used with the MODEL command in PCI will batch subset a directory of input files based on a defined input window. It assumes that all the input files will all be located within the same specified directory, all the files will be of the same format with three channels and that the output directory does not already contain any files.

This was originally created for PCI Geomatica v9.1 but was last tested and working with no problems in PCI Geomatica v10.0 – Just copy the code below into a blank text file and edit as needed.

!—————————————————————————-
!—————————————————————————-
! Multiple File Subset Script [mul_subset.eas]
!
! This script was originally written for a user who wanted
! to have all of the fiducial marks and details removed from
! their air photos.
!
! This script assumes that all the input files will all be located
! within a given directory, all the files will be of the same format
! with three channels and that the output directory does not contain
! any files.
!
!—————————————————————————-

! Define variables
!—————————————————————————-

!for input and output directory local string in_files, out_files
!for directory listing of the input directories

local mstring dirlist

!for the file format and extension types

local string type, ext

!file names

local string bn, fn

!Parameters for the subset window

local integer x, y, width, height

local integer i

local $Z

local string confirm

!—————————————————————————-
! Clear the EASI window and then show the header information
!—————————————————————————-

PRINT @(1 ,1,CLREOS)

print “——————————————————————-“

print @reverse,” Multiple Subset “,@alloff

print “”
print “This script assumes that all the input files will all be located”
print “within a given directory, the files will be of the same format,”
print “have three image channels, will be clipped to the same extents and”
print “that the output directory will not contain any files.”

print “”

print “——————————————————————-“

!—————————————————————————-
! Collect input from user
!—————————————————————————-

print “”

print “Enter the directory that contains the input files:”

input “>” in_files

print “”

print “Enter the directory for the output files:”

input “>” out_files

print “”

print “Enter the file format of the files (3-letter file extension) : “

input “>” type

print “”
print “Enter the X coordinate for the dbiw parameter : “

input “>” x

print “”
print “Enter the Y coordinate for the dbiw parameter : “

input “>” y

print “”
print “Enter the width for the dbiw parameter : “

input “>” width

print “”
print “Enter the height for the dbiw parameter : “

input “>” height

print “”

PRINT @(1 ,1,CLREOS)

!—————————————————————————-
! Confirm with user to ensure that the parameters are correct
! If they are correct then continue with the script and if they are not
! then run the script over again.
!—————————————————————————-

print “—————————————————————————–“
print “”
print “The input directory you specified was:”

print ” “, in_files

print “”

print “The output directory you specified was:”

print ” “, out_files

print “”

print “The file format you specified was:”

print ” “, type

print “”

print “The initial coordinate of the clip you specified was:”

print ” “, x, ” “, y

print “”

print “Your image size you specified was:”
print ” “, width, ” “, height
print “”
print “—————————————————————————–“

print “Are these parameters correct? (Y/N)”

print “”

input “>” confirm

if (confirm ~= “y” or confirm ~= “Y” ) then

!——————————————————————-
! Get the contents of the directory
!——————————————————————-

dirlist = getdirectory(in_files)

let $Z = “\

for i = 1 to f$len(dirlist)

!——————————————————————-
! Extract parts of the filenames
!——————————————————————-

fn = in_files + $Z + dirlist[i]

ext = getfileextension(fn)

bn = getfilebasename(fn)

if (ext ~= type) then

print “”

print “Clipping:”, bn, “.”, ext

print “”

!——————————————————————-
! Set up the parameters and execute the FEXPORT command
!——————————————————————-

fili = in_files + $Z + dirlist[i]
filo = out_files + $Z + bn +”_clipped”
dbiw = x, y, width, height
dbic = 1,2,3
dbib =
dbvs =
dblut =
dbpct =
ftype = type
foptions =

R Fexport

PRINT @(1 ,1,CLREOS)

endif

endfor

else

run “mul_subset.eas

endif

PRINT @(1 ,1,CLREOS)

print “——————————————————————-“
print “”
print “The clipped files are stored in the following directory:”
print ” “, out_files
print “”
print “”

print @reverse,” Multiple Subset EASI Script Finished “,@alloff

print “”
print “——————————————————————-“

return

!—————————————————————————-
!—————————————————————————-

 

Script to generate a PIX file large enough to contain all input DEM files

Below is a simple EASI script used with the MODEL command in PCI will generate an output PIX file with extents large enough to contain all of the input files PIX files and then mosaic each of the input files into the output pix file. It was designed to mosaic several DEM files into one large DEM file.

This was originally created for PCI Geomatica v9.1 but was last tested and working with no problems in PCI Geomatica v10.0 – Just copy the code below into a blank text file and edit as needed.

!—————————————————————————-

! Multiple DEM Mosaic Script
! [mul_DEM_mosaic.eas]

! Ted MacKinnon – tmackinnon.com
!
! This script will generate an output PIX file with extents large enough
! to contain all of the input files and then mosaic each of the input
! files into the new output pix file.
!
! This script was designed to mosaic several DEM files into one large
! DEM file. Each input DEM file must have the same projection,
! resolution and at each DEM should have at least 150m overlap.
!
! The user is expected that the user will enter (or verify) the parameters
! of ‘CIMPRO’ at the bottom of this script prior to running the script.
! Parmaters such as BXPXSZ and DBIC have been hard coded.

!—————————————————————————-
!—————————
! Define variables
!—————————

!for input & output pix files and pixel size

local string inputFile, outputFile
local string pixel

!to obtain list of input files
local mstring inputList
local integer inputIndex

!for extracting georeferencing information

local GeoInfo geoInfo
local integer geoFile

!to contain bounding rectangle of all input files

local double boundULX, boundULY, boundLRX, boundLRY

!—————————————————————————-
! Clear the EASI window and then show the header information
!—————————————————————————-

PRINT @(1 ,1,CLREOS)

print “———————————————————————–

print @reverse,” ‘Mosaic multiple DEM files’ EASI Script “,@alloff
print “”
print “”
print “”
print “———————————————————————–“
print “”
print ” This script will generate an new output PIX file with extents large “
print ” enough to contain all of the PIX files in this directory and then “
print ” mosaic each of the PIX files into the new output file.”
print “”
print ” The output file will contain the following projection: “
print “”

! The projection will be printed to the screen as a general reminder

print ” UTM 11 S E000″

print “”
print “———————————————————————–“
print “”

!—————————————————————————-
! Collect input from user
!—————————————————————————-

print “”
print “Enter the Output file name:”
input “>” outputFile
print “”

print “Enter the pixel size for: “,outputfile

input “>” pixel

print “”
print “”

PRINT @(1 ,1,CLREOS)

!—————————————————————————-
! Create list of input files
!—————————————————————————-

sys “dir *.pix /b > pixlist.txt”

inputList = Text$Import(“pixlist.txt”)

!—————————————————————————
! Read georeferencing of each input file
!—————————————————————————-

for inputIndex = 1 to F$LEN(inputList)

print inputList(inputIndex)

geoFile = DBOpen(inputList(inputIndex), “r”)
call DBReadGeoInfo(geoFile, geoInfo)
call DBClose(geoFile)

print “Georeferencing: “, geoInfo.Units
print “Image extents: “,geoInfo.ULX, “, “, geoInfo.ULY, ” “, geoInfo.LRX, “, “, geoInfo.LRY

PRINT @(1 ,1,CLREOS)

!————————————————————————-
! initialize bounds using first file
!————————————————————————-

if inputIndex = 1 then

boundULX = geoInfo.ULX
boundULY = geoInfo.ULY
boundLRX = geoInfo.LRX
boundLRY = geoInfo.LRY

else

!———————————————————————–
! assumes UTM projection (ULY > LRY)
!———————————————————————–

if geoInfo.ULX < boundULX then
boundULX = geoInfo.ULX

endif

if geoInfo.ULY > boundULY then

boundULY = geoInfo.ULY

endif

if geoInfo.LRX > boundLRX then

boundLRX = geoInfo.LRX

endif

if geoInfo.LRY < boundLRY then
boundLRY = geoInfo.LRY

endif

endif

print “”

endfor

print outputFile, ” file extents: “, boundULX, ” “, boundULY, ” “, boundLRX, ” “, boundLRY

!—————————————————————————-
! create the output PIX file containing bounds
!—————————————————————————-

FILE = outputFile
TEX1 =

!Specify the # of channels here

DBNC = 1,0,0,0
DBLAYOUT = “PIXEL”

!Specify the projection info here

PROJECT = “UTM”
ZONE = 11
ROW = “S”
ELLIPS = “0”
LLBOUND = “N”

ULX = F$STRING(boundULX)
ULY = F$STRING(boundULY)
LRX = F$STRING(boundLRX)
LRY = F$STRING(boundLRY)
BXPXSZ = pixel
BYPXSZ = pixel
REPORT = “TERM”

run CIMPRO

!—————————————————————————-
! Mosaic each PIX file from the listing into the outputfile
!—————————————————————————-

inputList = Text$Import(“pixlist.txt”)

for inputIndex = 1 to F$LEN(inputList)

inputFile = inputList(inputIndex)

!————————————————————————–
! Clear the EASI window and show progress
!————————————————————————–

PRINT @(1 ,1,CLREOS)

print “”
print “Mosaicking “, inputFile, ” into “, outputFile
print “”

FILI = inputFile
DBIC = 1
DBVS =
DBLUT =
FILO = outputFile
DBOC = 1
BLEND =
BACKVAL = 0

run MOSAIC

endfor
!—————————————————————————-
!—————————————————————————-

Script to clip the unnecessary excess portions of an air photo mosaic

Below is a simple EASI script that when used with the MODEL command in PCI will clip the unnecessary excess portions of an air photo mosaic created from Ortho Engine to an irregular buffered shape around a set study area. This was originally created for PCI Geomatica v9.1 but was last tested and working with no problems in PCI Geomatica v10.0 – Just copy the code below into a blank text file and edit as needed.

!——————————————————————————————–

! Irregular_Polygon_Clip [bitmap_clip.eas]

! Ted MacKinnon – tmackinnon.com

! This simple script used with the MODEL command in PCI EASI

! will clip the unnecessary excess portions of an airphoto mosaic

! to an irregular buffered shape around the study area.

!

! The working file ‘working-file.pix’ has the existing

! mosaic image located in the first three channels,

! an existing irregular shaped polygon bitmap and

! three empty 8 bit channels.

!

! %%2 is the bitmap mask of the irregular shape

! %1, %2, %3 are the RGB existing image channels

! %4, %5, %6 will be the new modeled RGB image channels

! the RGB value of 255, 255, 255 will set the background to white

!

! Simply change the file name to use this script with a different

! file and ensure that the channels and segments are setup the same

!——————————————————————————————–

!——————————–

! Set up and run the model

!——————————–

MODEL ON “working-file.pix” OVER dbiw

if %%2 = 1 then

%4 = %1;

%5 = %2;

%6 = %3;

else

%4 = 255;

%5 = 255;

%6 = 255;

endif;

ENDMODEL

!———————————————-

! Export the resultant channels to a new file

!———————————————-

FILI = “working-file.pix

FILO = “mosaic_clipped.tif

DBIW =

DBIC = 4,5,6

DBIB =

DBVS =

DBLUT =

DBPCT =

FTYPE = “TIF

FOPTIONS = “”

RUN FEXPORT

!——————————————————————————————–

 

Color Shaded Relief Model – Gatineau, Quebec

Color Shaded Relief Model of the Gatineau foot hills of Quebec

 

The above image is one of several color shaded relief (CSR) models I created of Gatineau foot hills of Quebec. The city of Ottawa is featured in the national capital region near the bottom center of the image across the river just south of Gatineau Park. This was one of the color shaded relief models submitted for use with an online interactive web GIS mapping project for the municipal of Collines-de-l’Outaouais, Quebec. The image below is the actual one that is currently featured in the web GIS project.

The Gatineau Hills are part of a geological formation in Canada which represent the foothills of the Laurentian Mountains which stretch east through Quebec, beginning north of Montreal and joining up with others into Vermont and New Hampshire. The geology of Gatineau Park, which encompass these foothills, is related to the Eardly Escarpment, which is a fault line that lies along the southern edge of the hills. This escarpment makes the part an attractive location for rock climbers and hikers, offering a beautiful view of the relatively flat fields below, which extend to the Ottawa River.

Color Shaded Relief Model of the Gatineau foot hills of Quebec

 

 

Color Shaded Relief Model of Vancouver Island

Color Shaded Relief Model of Vancouver Island, British Columbia

The above image is a Color shaded relief (CSR) model created using a DEM of Vancouver Island, British Columbia using PCI Geomatica software. The change in relief ranges from 0 m (blue) up to over 2200 m (red colors). Vancouver Island is a large island in the province of British Columbia off the Pacific coast of Canada, it is the largest island (32,134 square kilometers) on the western side of both North and South America.

This image has been published in various books, University Course material (such as UVic) and papers such as the Vancouver Sun & the Montreal Gazette

More examples of Color shaded relief (CSR) models …

Color Shaded Relief Models – Pointe Du Chene, New Brunswick

Color Shaded Relief Models of Pointe Du Chene, New BrunswickColor Shaded Relief Models of Pointe Du Chene, New Brunswick with flood
The above two images were created for my LIDAR flood modeling graduate research project. The first image is before the flood scenario; featuring a color shaded relief perspective view pointing south east from the Northumberland Strait landwards across the Pointe Du Chene wharf. The second image is of the same color shaded relief perspective view but features a 2.55 m flood level super imposed on top of it.

The 2.55 m flood level was an actual recorded storm surge water level that effected this area during a winter storm on January 2001. The two images below show the same flood level and area but from an overhead aerial view. The first image is with an orthophoto and the second image is with the color shaded relief.

Pointe Du Chene, New Brunswick with floodPointe Du Chene, New Brunswick with flood