Take me to pcigeomatics.com

PCI Geomatics Help Center

How can we help you today?

mul_DEM_mosaic

PCI Geomatics -

!----------------------------------------------------------------------------! 
!----------------------------------------------------------------------------!
!                                                                            !
!       	Multiple DEM Mosaic Script [mul_DEM_mosaic.eas]              !
!                                                                            !
!                                                                            !
!  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

!---------------------------------------------------------------------------- 
!---------------------------------------------------------------------------- 
Have more questions? Submit a request

Comments

Powered by Zendesk