[Gdal-dev] Fwd: Re: Support for ASTER AST_07 products
Gerald Buckmaster
gbuckmaster at cox.net
Fri Jun 13 02:29:27 EDT 2003
Hmm...maybe I should have been subscribed sooner...
Anyone else having trouble with ASTER Level 2 Surface Reflectance products?
---------- Forwarded Message ----------
Subject: Re: Support for ASTER AST_07 products
Date: Thursday 12 June 2003 07:43 pm
From: Gerald Buckmaster <gbuckmaster at cox.net>
To: openev-discuss at lists.sourceforge.net
Cc: Andrey Kiselev <dron at ak4719.spb.edu>
Andrey,
First, I am not a C/C++ programmer, nor very smart about compiling code so
mostly stabbing in the dark here.
After reviewing hdf4dataset.h, hdf4imagedataset.cpp, and hdf4dataset.cpp,
it appears to me AST_07 ASTER Level 2 Surface Reflectance products are
specifically NOT supported yet, unlike ASTER_L1A, ASTER_L1B, and others. It
also appears one reason I am able to open and extract bands is because the
poDS->iDataType and pszDataType is being handled as "UNKNOWN". I would be
hard-pressed to extract bands using gdal_translate without knowing the
SubdatasetTypes were "UNKNOWN".
I double-checked my gdalinfo against L1B and L1A and GeoTIFFs produced with
gdal_translate from HDF. No problems...just faulting on AST_07 products.
gdal_translate also does not appear to be extracting the projection or
coordinates from the AST_07 metadata.
I've included a vbscript that I and others wrote to handle various forms of
ASTER, including AST_07. Don't know if it will help, but it may give you
some additional ideas. plaster was really just a cli frontend to other
executables (in Windoze) to ease the extraction, georeferencing and merging
of ASTER image bands from hdf files. If I remember right, AST_07 bandnames
are Band1, Band2, Band3N, etc.
HTH and thanks for all your hard work.
Gerald
On Thursday 12 June 2003 08:54 am, Gerald Buckmaster wrote:
> Andrey,
>
> I uploaded the gdalinfo output for a ASTER L1B dataset that was acting
> strangely. it is at ftp.remotesensing.org/incoming/aster-l1b.txt
>
> I am actually using the nightly build from the 10th of June.
>
> The bzip2 brought the AST_07 dataset that includes SWIR and TIR bands down
> to 13mb. Much better than 78mb! Hope this doesn't clog your pipes too
> long. Unfortunately, I'm not technically skilled enough to extract only the
> TIR bands into a separate HDF...and that would probably negate
> troubleshooting anyway.
>
> You can find the AST_07 also at ftp.remotesensing.org/incoming/
>
> Thanks for the quick response!
>
> Gerald
>
> On Thursday 12 June 2003 07:14 am, you wrote:
> > Gerald,
> >
> > On Thu, Jun 12, 2003 at 06:39:49AM -0700, Gerald Buckmaster wrote:
> > > I noticed gdal_translate georeferences ASTER L1B HDF-EOS products very
> > > nicely, though there tends to be a initial display bug in
> > > OpenEV...something to do with what corner to start counting pixels
> > > from.
> >
> > It is very strange and should not happen. OpenEV works good for me, but
> > it is possible there is a bug in HDF4 code. Could you show me gdalinfo
> > output for those ASTER L1B dataset?
> >
> > > My real need, and I believe most public consumers of ASTER need, is the
> > > ability to do the same with ASTER Level 2 surface reflectance data set
> > > (AST07) products...gdal_translate works, but without interrogating the
> > > metadata and georeferencing. In fact, gdalinfo
> > > gives me a Segmentation fault runtime on AST07 HDF files.
> >
> > Very bad. What GDAL version do you have?
> >
> > Regarding AST_07 datasets: it will be good if you will upload sample
> > file into incoming. Try to compress it with bzip2, it should get some
> > effect. 90m dataset will be enough for me, I shall manage to order other
> > samples myself.
> >
> > Regards,
> > Andrey
-------------------------------------------------------
-------------- next part --------------
' PLASTER.VBS
' VBScript to process ASTER Satellite Imagery
' Written by Gerald Buckmaster, www.imagery-analyst.com
' Version 1.6 - 08 April 2002
' Version History
' v1.8 - 20 Aug 02 - Recoded to handle AST_07 SurfaceReflectance product
' v1.6 - 08 Apr 02 - Loose code optimized by Emilio, merged with previous enhancements of v1.5.2
' Fixed 'projection initialization failure by adding +ellps=clrk66 to proj.exe commandline
' Eliminated need to for ensuring batchproj.bat is included; script writes, executes, and deletes batchproj.bat
' v1.5.2 - 06 Apr 02 - Fixed script to determine processing level via coremetadata.0 vs. filename
' Combined Subs Geocode_VNIR_Bands(), Geocode_SWIR_Bands(), Geocode_TIR_Bands()
' into 1 subroutine called Geocode_Bands().
' Re-wrote subroutine Geocode_Image(BandNumber, SENSOR) to handle any band.
' Added Remove_A_File() to clean out the non-georeferenced tiff files.
' Added Command Line Argument processing - no longer need to edit this script
' v1.5.1 - 02 Apr 02 - Combined improvements below by G. Nelson and E. Mayorga
' v1.5 - 01 Apr 02 - Modified by Gerald Nelson
' Moved user-editable constants closer to top of script, added comments and a
' and a detailed seperate readmefirst.txt file.
' v1.4 - 15 Feb 02 - Added support for processing L1B Expedited Product. Fixed "double-negative"
' MAPORIENTATIONANGLE bug.
' 1.3.1- 06 Feb 02 - Modified by Emilio Mayorga
' Added extraction of TIR bands. Made some minor changes also; mostly cosmetic,
' but also fixed one apparent mistake with the name of the subroutine
' Render_GeoTiff_CartesianType_PointFile
' v1.3 - 16 Jan 02 - Adjusted PullMapOrientation() to use a default -10 degrees in case the
' Metadata keyword MAPORIENTATIONANGLE or the seemingly rarer SCENEORIENTATIONANGLE
' is not found. -10 degrees seems to be a fair average orientation angle.
' v1.2 - 15 Jan 02 - Able to extract map orientation from Level-1B and conduct transformation
' This rotates the image north up. We are now able to geocode in either
' geographic lat-longs or utm cartesian coordinates
' v1.1 - 08 Jan 02 - Modified to process level 2 ASTER hdf files: AST_06V and AST_06S
' v1.0 - 06 Jan 02 - Initial working release of PLASTER.VBS works with level 1b ASTER hdf files
' This script accomplishes the following:
' Examines metadata from ASTER HDF files.
' Determines Level of Processing and Band Availability.
' Extracts and writes geographic coordinates to an ASCII file.
' Extracts map orientation angle.
' Provides a choice between rendering cartesian or lat-long coordinate system
' Calculates hemisphere and UTM zone.
' Converts the geographic coordinates to cartesian coordinates.
' Renders a geotiff modeltransformation definition file.
' Renders a geotiff four-corner tiepoint definition file.
' Extracts individual image bands to TIFF format.
' Geocodes each TIFF with definition file.
' This script assumes the following files are in your system with appropriate paths:
' vbscript file - this file
' cscript.exe - Windows Scripting Host Command Line Interface
' Copyright (C) Microsoft Corporation 1996-2001.
' hdp.exe - a HDF command-line utility
' ftp://ftp.ncsa.uiuc.edu/HDF/HDF/HDF_Current/zip/utilities/hdp.exe
' Copyright 1988-1999 The Board of Trustees of the University of Illinois
' proj.exe - a geographic projection command-line utility
' ftp://ftp.remotesensing.org/pub/proj/from_kai/DOS-PJ433.zip
' hdfbrowse.exe - a HDF command-line utility that extracts tif images from HDF
' http://www.scanex.ru/downloads/eostools/hdfbrowse.zip
' Copyright (C) 2001 R&D center ScanEx.
' geotifcp.exe - a utility to convert TIF to GeoTIFF
' http://www.possys.com/utilities.html
' gdal_merge.py - python script by Frank Warmerdam packaged with OpenEV
' Needed only if you want a single multiband GeoTiff file by not using
' the command line flag -s
' http://openev.sourceforge.net/
' The OpenEV version must be 1.41+ in order to ensure the correct gdal_merge.py is present.
Const OPENEV_HOME = "C:\openev_fw\"
Const ForReading = 1, ForWriting = 2, ForAppending = 3
Const pi = 3.1415926535897932
Const DegToMeters = 111111.111 ' meters -> degrees conversion factor
Const VNIR = 0, SWIR = 1, TIR = 2 ' Band Group Indices
Const BandGroupTot = 3
Const BandTot = 15
dim strOutputFile
dim strMetaFile
dim strShellCmd
dim UTMZone
dim Hemi_1
dim Hemi_2
Dim cartesian_file
Dim Cart_Coord_UL
Dim Cart_Coord_UR
Dim Cart_Coord_LL
Dim Cart_Coord_LR
Dim proj_def_file
Dim PIX_COORDS (3)
Dim PIX_RES
Dim Geog_Coord_UL
Dim Geog_Coord_UR
Dim Geog_Coord_LL
Dim Geog_Coord_LR
Dim strCoordFile
Dim SDSString
Dim MapRotRads
Dim Processing_Level
Dim BandGroupIdx, BandIdx
Dim PixResScl, GeoTiffCoordPar
Dim bBandGroupProcess(3)
Dim BandGroupName
Dim BandGroupRows, BandGroupCols
Dim BandGroupBandStart, BandGroupBandEnd
Dim PixRES_meters
Dim BandName
Dim BaseFilePath
Dim Rendered_File_Count
Dim GDAL_Input(14)
Dim strHDFname
Dim drvPath
Dim drvDestPath
Dim GeographicType
Dim ModelTransformation
Dim MergeBands
Dim objArgs
Dim TmpFileCnt
Dim TempFileList(10)
Dim BNDFLAG
Dim BSTPROJ
' ASTER Sensor and Dataset defaults
BandGroupName = Array("VNIR","SWIR","TIR")
BandGroupRows = Array(4200,2100,700)
BandGroupCols = Array(4980,2490,830)
PixRES_meters = Array(15,30,90)
BandGroupBandStart = Array(0,4,10)
BandGroupBandEnd = Array(2,9,14)
BandName = Array("1","2","3N","3B","4","5","6","7","8","9","10","11","12","13","14")
TmpFileCnt = 0
' PROCESS COMMAND-LINE ARGUMENTS
' cscript plaster16.vbs -fhdfname -iinputpath -ooutputpath [-g] [-t] [-s] [-evo|eso|eto|evs|est|etv]
' 0 1 2 3 4 5 6 7 8
' -f (name of hdf file)
' -i (directory of input hdf file)
' -o (directory for output files)
' -g (use geographic coordinate system vs. UTM cartesian system - the assumed default)
' -t (use ModelTiepointTag vs. ModelTransformationTag - the assumed default)
' -s (separate bands extraction (no merge) vs. merge of bands selected for processing - the assumed default)
' All bands are extracted if flags below are not utilized
' -evo (extract VNIR only)*
' -eso (extract SWIR only)
' -eto (extract TIR only)
' -evs (extract VNIR & SWIR)*
' -est (extract SWIR & TIR)
' -etv (extract TIR & VNIR)*
' * Band 3B is not processed
' Default Settings
GeographicType = False ' use optional -g flag to use geographic coordinate system
ModelTransformation = True ' use optional -t flag to use ModelTiePointTag
MergeBands = True ' use optional -s flag to prevent merging of extracted bands into a single GeoTIFF
' use optional flags -evo, -eso, -eto, -evs, -est, -etv to control which bands are extracted
BNDFLAG = "VST"
BSTPROJ = "VNIR"
bBandGroupProcess(VNIR) = True
bBandGroupProcess(SWIR) = True
bBandGroupProcess(TIR) = True
'Parse out the Command Line Arguments
Set objArgs = Wscript.Arguments
For I = 0 to objArgs.Count - 1
Select Case left(objArgs.Item(I),2)
Case "-f"
strHDFname = mid(objArgs.Item(I),3)
Case "-i"
drvPath = mid(objArgs.Item(I),3)
If right(drvPath,1) <> "\" then drvPath = drvPath & "\"
Case "-o"
drvDestPath = mid(objArgs.Item(I),3)
If right(drvDestPath,1) <> "\" then drvDestPath = drvDestPath & "\"
Case "-g"
GeographicType = True
Case "-t"
ModelTransformation = False
Case "-s"
MergeBands = False
Case "-e"
Select Case left(objArgs.Item(I),4)
Case "-evo"
BNDFLAG = "V"
BSTPROJ = "VNIR"
bBandGroupProcess(SWIR) = False
bBandGroupProcess(TIR) = False
Case "-eso"
BNDFLAG = "S"
BSTPROJ = "SWIR"
bBandGroupProcess(VNIR) = False
bBandGroupProcess(TIR) = False
Case "-eto"
BNDFLAG = "T"
BSTPROJ = "TIR"
bBandGroupProcess(VNIR) = False
bBandGroupProcess(SWIR) = False
Case "-evs"
BNDFLAG = "VS"
BSTPROJ = "VNIR"
bBandGroupProcess(TIR) = False
Case "-est"
BNDFLAG = "ST"
BSTPROJ = "SWIR"
bBandGroupProcess(VNIR) = False
Case "-etv"
BNDFLAG = "VT"
BSTPROJ = "VNIR"
bBandGroupProcess(SWIR) = False
End Select
End Select
Next
' Check for missing mandatory arguments
If strHDFname ="" Then
Wscript.Echo "Missing the -f argument...you need to provide the name of the HDF file."
WriteCommandSyntax
Wscript.Quit
End If
If drvPath ="" Then
Wscript.Echo "Missing the -i argument...you need to provide the entire path to the HDF file."
WriteCommandSyntax
Wscript.Quit
End If
If drvDestPath ="" Then
Wscript.Echo "Missing the -o argument...you need to an output path for the resulting GeoTIFF files."
WriteCommandSyntax
Wscript.Quit
End If
Rendered_File_Count = 0 ' initialize the count
BaseFilePath = drvDestPath & strHDFname & "_"
' Determine Level of Processing represented in HDF file
strMeta = "coremetadata.0"
Extract_MetaData (strMeta)
Processing_Level = PullShortName (strMetaFile, "SHORTNAME")
' No matter what bands the user asked for, some bands will not be available
' At the same time, we cannot dictate what bands the user wants processed
' This filter could still be improved by parsing metadata for actual bands present
Select Case Processing_Level
Case "AST_07" 'Added 2002-08-19 for processing Surface Reflectance Product
'Situation: AST_07 product arrives in two hdf files, separating
' VNIR and SWIR, but AST_07 is used in each file.
strMeta = "productmetadata.0"
Extract_MetaData (strMeta)
Processed_Bands = PullShortName (strMetaFile, "PROCESSEDBANDS")
if Processed_Bands = "XXXXXXXX040506070809XXXXXXXXXX" then
BSTPROJ = "SWIR"
bBandGroupProcess(VNIR) = False
bBandGroupProcess(TIR) = False
end if
if Processed_Bands = "01023NXXXXXXXXXXXXXXXXXXXXXXXX" then
BSTPROJ = "VNIR"
bBandGroupProcess(SWIR) = False
bBandGroupProcess(TIR) = False
end if
Case "AST_06V"
BSTPROJ = "VNIR"
bBandGroupProcess(SWIR) = False
bBandGroupProcess(TIR) = False
Case "AST_06S"
BSTPROJ = "SWIR"
bBandGroupProcess(VNIR) = False
bBandGroupProcess(TIR) = False
Case "ASTL1B"
' Every band could actually be available for processing
End Select
' Extract coordinates from HDF file - writes a text file for further processing
strMeta = "productmetadata.0"
Extract_MetaData (strMeta)
' Find and stores geographic latitudes and longitudes to arrays
Read_Geographic_Coords Geog_Coord_UL, "UPPERLEFT"
Read_Geographic_Coords Geog_Coord_UR, "UPPERRIGHT"
Read_Geographic_Coords Geog_Coord_LL, "LOWERLEFT"
Read_Geographic_Coords Geog_Coord_LR, "LOWERRIGHT"
' Write Latitude-Longitude coordinates to file
strCoordFile = BaseFilePath & "coords.txt"
WriteCoordFile strCoordFile
' Pull Map Orientation data & determine proper radian
MapRotRads = (pi/180) * PullMapOrientation(strMetaFile)
If GeographicType = True Then 'user wants images with geographic coordinate system (latitude longitude pairs)
PixResScl = 1 / DegToMeters
GeoTiffCoordPar = 0
Else ' User wants images projected to utm cartesians (eastings & northings)
PixResScl = 1
GeoTiffCoordPar = 1
If Trim(Geog_Coord_UL(1)) > 0 then ' Determine hemisphere of image location - either north or south
Hemi_1 = "north"
Hemi_2 = "N"
Else
Hemi_1 = "south"
Hemi_2 = "S"
End if
' Convert Latitude-Longitude pairs to UTM Cartesians
UTMZone = Calculate_UTM_Zone(Geog_Coord_LL(0), Geog_Coord_UR(0))
Convert_LatLon_2_Cartesians Hemi_1, UTMZone, strCoordFile ' added missing Hemi_1 08 April 02
Read_Cartesian_Coords
End If
'Create temporary GeoTIFF projection definition file for each Band Group requested
If ModelTransformation = True Then ' We are using rotation factor
For BandGroupIdx = 0 To BandGroupTot - 1
If bBandGroupProcess(BandGroupIdx) = True Then
PIX_COORDS(0) = PixResScl * PixRES_meters(BandGroupIdx) * sin(MapRotRads)
PIX_RES = PixResScl * PixRES_meters(BandGroupIdx) * cos(MapRotRads)
Render_GeoTiff_File BandGroupName(BandGroupIdx)
End If
Next
Else ' We are using tiepoints for corner coordinates
For BandGroupIdx = 0 To BandGroupTot - 1
If bBandGroupProcess(BandGroupIdx) = True Then
PIX_COORDS(0) = GeoTiffCoordPar & " " & GeoTiffCoordPar & " 0"
PIX_COORDS(1) = BandGroupCols(BandGroupIdx) & " " & GeoTiffCoordPar & " 0"
PIX_COORDS(2) = GeoTiffCoordPar & " -" & BandGroupRows(BandGroupIdx) & " 0"
PIX_COORDS(3) = BandGroupCols(BandGroupIdx) & " -" & BandGroupRows(BandGroupIdx) & " 0"
PIX_RES = PixResScl * PixRES_meters(BandGroupIdx)
Render_GeoTiff_File BandGroupName(BandGroupIdx)
End If
Next
End If
' Here's the magic, courtesy of hdfbrowse and geotifcp
For BandGroupIdx = 0 To BandGroupTot - 1
If bBandGroupProcess(BandGroupIdx) = True Then
Extract_ASTER_Bands BandGroupIdx
Geocode_Bands BandGroupIdx
End If
Next
' More magic, if you want a single, multiband GeoTIFF
If MergeBands = True Then
Render_GDAL_Merge_Batch_File
strShellCmd = drvDestPath & strHDFname & ".bat"
RunShell strShellCmd
End If
' Clean Up unnecessary files
' Remove all the individual band files geotiff files
For X = 1 To Rendered_File_Count
Remove_A_File GDAL_Input(X)
Next
Remove_A_File drvDestPath & strHDFname & "_TEMPGDAL.tif"
Remove_A_File drvDestPath & strHDFname & "_cartesians.txt"
Remove_A_File drvDestPath & strHDFname & "_coords.txt"
Remove_A_File drvDestPath & strHDFname & "_" & BSTPROJ & "_projinfo.txt"
Remove_A_File drvDestPath & strHDFname & ".bat"
Remove_A_File drvDestPath & strHDFname & "_batchproj.bat"
Remove_A_File drvDestPath & strHDFname & "_coremetadata.0.met"
Remove_A_File drvDestPath & strHDFname & "_productmetadata.0.met"
' End of Main Program
' =============================================================================
' SUBROUTINES & FUNCTIONS TO CALL FROM MAIN PROGRAM
Sub Extract_MetaData (strMetaName)
' Added to allow dumping of both coremetadata.0 and productmetadata.0 - modified 06 April 2002 - buckmaster
' Requires presence of hdp.exe in $PATH
strOutputFile = BaseFilePath & strMetaName & ".met"
strShellCmd = "hdp dumpvd -n " & strMetaName & " -d -o " & strOutputFile & " -x " & drvPath & strHDFname
RunShell strShellCmd
' Clean data of extraneous characters including spaces
strMetaFile = ReadNewMetaFile(strOutputFile)
strMetaFile = StripSpaces(strMetaFile)
End Sub
Sub Extract_ASTER_Bands(CurrBGIdx)
' Requires presence of hdfbrowse.exe in $PATH
For BandIdx = BandGroupBandStart(CurrBGIdx) To BandGroupBandEnd(CurrBGIdx)
Select Case Processing_Level
Case "AST_07" 'Added 20020819
if Processed_Bands = "01023NXXXXXXXXXXXXXXXXXXXXXXXX" then
SDSString = chr(34) & "SurfaceReflectanceVNIR!Data Fields!Band" & BandName(BandIdx) & chr(34)
end if
if Processed_Bands = "XXXXXXXX040506070809XXXXXXXXXX" then
SDSString = chr(34) & "SurfaceReflectanceSWIR!Data Fields!Band" & BandName(BandIdx) & chr(34)
end if
Case "ASTL1B"
SDSString = chr(34) & BandGroupName(CurrBGIdx) & "!" & BandGroupName(CurrBGIdx) & "_Swath!Data Fields!ImageData" & BandName(BandIdx) & chr(34)
Case "AST_06V"
SDSString = chr(34) & "DecorrelationStretchVNIR!Data Fields!Band" & BandName(BandIdx) & chr(34)
Case "AST_06S"
SDSString = chr(34) & "DecorrelationStretchSWIR!Data Fields!Band" & BandName(BandIdx) & chr(34)
End Select
' Extract image band from ANY level of hdf
strShellCmd = "hdfbrowse.exe -f" & drvPath & strHDFname & " -o" & drvDestPath & strHDFname & "_" & BandName(BandIdx) & ".tiff -s" & SDSString & " -b1 -d1 -x"
RunShell(strShellCmd)
Next
End Sub
Sub Remove_A_File (oldfile) ' Added 06 Apr 02 - buckmaster
' Facilitate the cleanup of unneeded raster and text files temporarily generated
Set fso = CreateObject("Scripting.FileSystemObject")
Set a = fso.GetFile(oldfile)
a.Delete
End Sub
Sub Geocode_Bands(CurrBGIdx)
Dim InTifFile
Dim OutTifFile
Dim GeoProjFile
For BandIdx = BandGroupBandStart(CurrBGIdx) To BandGroupBandEnd(CurrBGIdx)
InTifFile = BaseFilePath & BandName(BandIdx) & ".tiff"
OutTifFile = BaseFilePath & BandName(BandIdx) & "_geo.tif"
GeoProjFile = BaseFilePath & BandGroupName(CurrBGIdx) & "_projinfo.txt"
' Geocode single existing image with rendered geotiff projection file
strShellCmd = "geotifcp -g " & GeoProjFile & " " & InTifFile & " " & OutTifFile
RunShell strShellCmd
Remove_A_File InTifFile
Rendered_File_Count = Rendered_File_Count + 1
' Populate an array with a file geotiff to be merged with Render_GDAL_Merge_Batch_File
GDAL_Input(Rendered_File_Count) = OutTifFile
Next
End Sub
Sub Render_GDAL_Merge_Batch_File()
' Requires presence of gdal_merge.py, thus the presence of OpenEV
' A command is written here to re-geotifcp the merged tif
' for some reason, gdal_merge.py corrupts the geotiff transformation tag
Dim strMergeFiles 'list of files to merge
strMergeFiles = "" 'initially null
For X = 1 To Rendered_File_Count
InputCheck = IsNull(GDAL_Input(X))
' If Array item is not empty, add geotiff filename to merge list
If InputCheck = False Then strMergeFiles = strMergeFiles & " " & GDAL_Input(X)
Next
Set fso = CreateObject("Scripting.FileSystemObject")
Set temp_batch = fso.CreateTextFile(drvDestPath & strHDFname & ".bat", True)
' Write coords to file
temp_batch.WriteLine "@echo off"
temp_batch.WriteLine "set OPENEV_HOME=" & OPENEV_HOME
temp_batch.WriteLine "set OLD_PATH=%PATH%"
temp_batch.WriteLine "set path=%OPENEV_HOME%\bin;%OPENEV_HOME%\python;%PATH%"
temp_batch.WriteLine "set PYTHONPATH=%OPENEV_HOME%\pymod;%OPENEV_HOME%\numpy"
temp_batch.WriteLine "pythonw %OPENEV_HOME%\pymod\gdal_merge.py -o " & BaseFilePath & "TEMPGDAL.tif -v -separate " & strMergeFiles
' Need to change "VNIR_projinfo.txt " to a variable that is set when determining band choices.
temp_batch.WriteLine "geotifcp -g " & BaseFilePath & BSTPROJ & "_projinfo.txt " & BaseFilePath & "TEMPGDAL.tif " & BaseFilePath & BNDFLAG & ".tif"
temp_batch.WriteLine "set PATH=%OLD_PATH%"
temp_batch.Close
End Sub
Sub Read_Geographic_Coords(Geog_Coords_Corner, Corner)
' Find and stores geographic latitudes and longitudes to arrays
Dim GeogrCoords
GeogrCoords = Split(PullCoord(strMetaFile, Corner), ",", -1, 1)
Geog_Coords_Corner = Array(GeogrCoords(1), GeogrCoords(0))
End Sub
Sub Read_Cartesian_Coords()
'Read in the cartesians from the file written by PROJ earlier
Dim fso, f
Dim mytab
cartesian_file = BaseFilePath & "cartesians.txt"
Set fso = CreateObject("Scripting.FileSystemObject")
Set f = fso.OpenTextFile(cartesian_file, 1)
mytab = Chr(9)
' coordinate order: upperleft - 0, upperright - 1, lowerleft - 2, lowerright -3
tempcoord = trim(f.ReadLine)
Cart_Coord_UL = Split(tempcoord,mytab, -1, 1)
chopped_length = len(Cart_Coord_UL(1))-1
Cart_Coord_UL(1) = mid(Cart_Coord_UL(1),1,chopped_length)
tempcoord = trim(f.ReadLine)
Cart_Coord_UR = Split(tempcoord,mytab, -1, 1)
chopped_length = len(Cart_Coord_UR(1))-1
Cart_Coord_UR(1) = mid(Cart_Coord_UR(1),1,chopped_length)
tempcoord = trim(f.ReadLine)
Cart_Coord_LL = Split(tempcoord,mytab, -1, 1)
chopped_length = len(Cart_Coord_LL(1))-1
Cart_Coord_LL(1) = mid(Cart_Coord_LL(1),1,chopped_length)
tempcoord = trim(f.ReadLine)
Cart_Coord_LR = Split(tempcoord,mytab, -1, 1)
chopped_length = len(Cart_Coord_LR(1))-1
Cart_Coord_LR(1) = mid(Cart_Coord_LR(1),1,chopped_length)
f.Close
End Sub
Sub Render_GeoTiff_File(CurrBGName)
' Will render for any row/column set, pixel set and cartesian UTM coord set
' But you still need to provide those sets...see how to set up variable for
' calling Render_GeoTiff_ProjFile in the code immediately following the call
' to Read_Cartesian_Coords(). Future Feature will allow SELECTIVE PROJECTIONS
Dim i, Coord_UL(2), Coord_UR(2), Coord_LL(2), Coord_LR(2)
Dim ModelTypeGeoKey, LastLineTypeGeoKey
If GeographicType = True Then
For i = 0 To 1
Coord_UL(i) = Geog_Coord_UL(i)
Coord_UR(i) = Geog_Coord_UR(i)
Coord_LL(i) = Geog_Coord_LL(i)
Coord_LR(i) = Geog_Coord_LR(i)
Next
' Setting actual GeoTIFF Tags - completely dependent upon the specification
ModelTypeGeoKey = "ModelTypeGeographic"
LastLineTypeGeoKey = "GeographicTypeGeoKey (Short,1): GCS_WGS_84"
Else
For i = 0 To 1
Coord_UL(i) = Cart_Coord_UL(i)
Coord_UR(i) = Cart_Coord_UR(i)
Coord_LL(i) = Cart_Coord_LL(i)
Coord_LR(i) = Cart_Coord_LR(i)
Next
' Setting actual GeoTIFF Tags - completely dependent upon the specification
ModelTypeGeoKey = "ModelTypeProjected"
LastLineTypeGeoKey = "ProjectedCSTypeGeoKey (Short,1): PCS_WGS84_UTM_zone_" & UTMZone & Hemi_2
End If
' Open File for writing
Set fso = CreateObject("Scripting.FileSystemObject")
Set proj_def_file = fso.CreateTextFile(drvDestPath & strHDFname & "_" & CurrBGName & "_projinfo.txt", True)
' Write coords to file
proj_def_file.WriteLine "Geotiff_Information:"
proj_def_file.WriteLine "Version: 1"
proj_def_file.WriteLine "Key_Revision: 1.0"
proj_def_file.WriteLine "Tagged_Information:"
If ModelTransformation = True Then
proj_def_file.WriteLine "ModelTransformationTag(4,4):"
proj_def_file.WriteLine PIX_RES & " " & PIX_COORDS(0) & " 0 " & Coord_UL(0)
proj_def_file.WriteLine PIX_COORDS(0) & " -" & PIX_RES & " 0 " & Coord_UL(1)
proj_def_file.WriteLine "0 0 0 0"
proj_def_file.WriteLine "0 0 0 1"
Else
proj_def_file.WriteLine "ModelTiepointTag (8,3):"
proj_def_file.WriteLine PIX_COORDS(0)
proj_def_file.WriteLine Coord_UL(0) & " " & Coord_UL(1) & " 0"
proj_def_file.WriteLine PIX_COORDS(1)
proj_def_file.WriteLine Coord_UR(0) & " " & Coord_UR(1) & " 0"
proj_def_file.WriteLine PIX_COORDS(2)
proj_def_file.WriteLine Coord_LL(0) & " " & Coord_LL(1) & " 0"
proj_def_file.WriteLine PIX_COORDS(3)
proj_def_file.WriteLine Coord_LR(0) & " " & Coord_LR(1) & " 0"
proj_def_file.WriteLine "ModelPixelScaleTag (1,3):"
proj_def_file.WriteLine PIX_RES & " " & PIX_RES & " 0"
End If
proj_def_file.WriteLine "End_Of_Tags."
proj_def_file.WriteLine "Keyed_Information:"
proj_def_file.WriteLine "GTModelTypeGeoKey (Short,1): " & ModelTypeGeoKey
proj_def_file.WriteLine "GTRasterTypeGeoKey (Short,1): RasterPixelIsArea"
proj_def_file.WriteLine LastLineTypeGeoKey
proj_def_file.WriteLine "End_Of_Keys."
proj_def_file.WriteLine "End_Of_Geotiff."
proj_def_file.Close
End Sub
Sub Convert_LatLon_2_Cartesians(Hemi_1, UTMZone, strCoordFile)
' Convert Latitude-Longitude pairs to UTM Cartesians
'Write batch file
Set fso = CreateObject("Scripting.FileSystemObject")
Set temp_batch = fso.CreateTextFile(BaseFilePath & "batchproj.bat", True)
' Write coords to file
temp_batch.WriteLine "@echo off"
temp_batch.WriteLine "proj +proj=utm +" & Hemi_1 & " +zone=" & UTMZone & " -r +ellps=clrk66 " & strCoordFile & " > " & BaseFilePath & "cartesians.txt"
temp_batch.Close
'Execute batch file
strShellCmd = BaseFilePath & "batchproj.bat"
RunShell strShellCmd
End Sub
Sub WriteCommandSyntax()
Wscript.Echo
Wscript.Echo "cscript plaster16.vbs -fhdfname -iinpath -ooutpath [-g] [-t] [-s] [-evo|eso|eto|evs|est|etv]"
Wscript.Echo
Wscript.Echo " -f (name of hdf file)"
Wscript.Echo " -i (directory of input hdf file)"
Wscript.Echo " -o (directory for output files)"
Wscript.Echo
Wscript.Echo " -g (use geographic coordinate system vs. UTM cartesian system - the assumed default)"
Wscript.Echo " -t (use ModelTiepointTag vs. ModelTransformationTag - the assumed default)"
Wscript.Echo " -s (separate bands extraction (no merge) vs. merge of bands selected for processing - the assumed default)"
Wscript.Echo
Wscript.Echo "All bands are extracted if no flags below are stated"
Wscript.Echo " -evo (extract VNIR only)*"
Wscript.Echo " -eso (extract SWIR only)"
Wscript.Echo " -eto (extract TIR only)"
Wscript.Echo " -evs (extract VNIR & SWIR)*"
Wscript.Echo " -est (extract SWIR & TIR)"
Wscript.Echo " -etv (extract TIR & VNIR)*"
Wscript.Echo
Wscript.Echo "* Band 3B is not processed in this version."
End Sub
Function RunShell(strShellCmd)
' Handles most of the calls to external DOS executables to include batch files
Set WshShell = Wscript.CreateObject("Wscript.Shell")
Return = WshShell.Run (strShellCmd, 1, TRUE)
End Function
Function ReadNewMetaFile(strOutputFile)
' Reads an entire file into a string for further parsing
Dim fso, f
Set fso = CreateObject("Scripting.FileSystemObject")
Set f = fso.OpenTextFile(strOutputFile, ForReading)
ReadNewMetaFile = f.ReadAll
End Function
Function WriteCoordFile(strFile)
' Write out the geographic coordinates to a file. I know it probably better
' to pipe it to other functions, but PROJ 4.4.4 is what I use to convert to
' UTM cartesians, and it doesn't like working well with anything but files
' for input...it's the ole STDIN/STDOUT process that's biting me.
Dim fso
Dim txtfile
Dim s
Set fso = CreateObject("Scripting.FileSystemObject")
Set txtfile = fso.CreateTextFile(strFile, True)
' Write coords to file
txtfile.WriteLine Trim(Geog_Coord_UL(1)) & " " & Trim(Geog_Coord_UL(0))
txtfile.WriteLine Trim(Geog_Coord_UR(1)) & " " & Trim(Geog_Coord_UR(0))
txtfile.WriteLine Trim(Geog_Coord_LL(1)) & " " & Trim(Geog_Coord_LL(0))
txtfile.WriteLine Trim(Geog_Coord_LR(1)) & " " & Trim(Geog_Coord_LR(0))
txtfile.Close
End Function
Function PullCoord(strMetaFile,strCoord)
' Got to pull geographic coordinates (degree-decimal format) before
Dim intStart, intEnd
intStart = Instr(strMetaFile,strCoord)
intStart = Instr(intStart,strMetaFile,"(")
intEnd = Instr(intStart,strMetaFile,")")
PullCoord = Mid(strMetaFile, intStart+1, (intEnd-intStart)-1)
End Function
Function PullProcessedBands(strMetaFile,strCoord)
' This is for determining what bands to process...very useful for AST_06S files
' This will assist SELECTIVE band processing - a future feature
Dim intStart, intEnd
intStart = Instr(strMetaFile,strCoord)
intStart = Instr(intStart,strMetaFile,chr(34))
intEnd = Instr(intStart+1,strMetaFile,chr(34))
PullProcessedBands = Trim(mid(strMetaFile, intStart, intEnd))
End Function
Function PullShortName(strMetaFile,strCoord) '- added 06 April 2002 - buckmaster
' This is for determining what type of processing the ASTER HDF-EOS experienced
Dim intStart, intEnd
intStart = Instr(strMetaFile,strCoord)
intStart = Instr(intStart,strMetaFile,chr(34))
intEnd = Instr(intStart+1,strMetaFile,chr(34))
PullShortName = Trim(mid(strMetaFile, intStart+1, intEnd-intStart-1))
End Function
Function PullMapOrientation(strMetaFile)
' This is for pulling orientation angle...differs with product :^(
Dim intStart, intEnd
Dim tempstring, OrientationDir
intStart = Instr(strMetaFile,"MAPORIENTATIONANGLE")
If intStart > 0 Then
OrientationDir = -1
Else
intStart = Instr(strMetaFile,"SCENEORIENTATIONANGLE")
If intStart > 0 Then
OrientationDir = 1
Else ' Use a default 10 degrees of angle
PullMapOrientation = -10.000
msgbox "No orientation metadata found...using default -10 degrees"
'Return ' Type mismatch: 'Return'
End If
End If
' Extract and Process the orientation angle
intStart = Instr(intStart,strMetaFile,"VALUE") ' Invalid procedure call or argument: 'Instr'
intStart = Instr(intStart,strMetaFile,"=")
intEnd = Instr(intStart+1,strMetaFile,"END_OBJECT")
tempstring = Trim(mid(strMetaFile, intStart+1, intEnd-intStart-5))
PullMapOrientation = OrientationDir * tempstring
End Function
Function StripSpaces(strMetaFile)
' Clean extraneous spaces from metadata - need to accomplish before parsing for any data
StripSpaces = Replace(strMetaFile, " ", "")
End Function
Function Calculate_UTM_Zone(W_Lon, E_Lon)
' Subroutine to Determine UTMZone
Dim Center_Lon
If W_Lon < 0 and E_Lon < 0 then
Center_Lon = E_Lon + ((W_Lon - E_Lon)/2)
End If
If W_Lon > 0 and E_Lon > 0 then
Center_Lon = W_Lon + ((E_Lon - W_Lon)/2)
End If
If Center_Lon < 0 and Center_Lon > -6 then Calculate_UTM_Zone = "30"
If Center_Lon < -6 and Center_Lon > -12 then Calculate_UTM_Zone = "29"
If Center_Lon < -12 and Center_Lon > -18 then Calculate_UTM_Zone = "28"
If Center_Lon < -18 and Center_Lon > -24 then Calculate_UTM_Zone = "27"
If Center_Lon < -24 and Center_Lon > -30 then Calculate_UTM_Zone = "26"
If Center_Lon < -30 and Center_Lon > -36 then Calculate_UTM_Zone = "25"
If Center_Lon < -36 and Center_Lon > -42 then Calculate_UTM_Zone = "24"
If Center_Lon < -42 and Center_Lon > -48 then Calculate_UTM_Zone = "23"
If Center_Lon < -48 and Center_Lon > -54 then Calculate_UTM_Zone = "22"
If Center_Lon < -54 and Center_Lon > -60 then Calculate_UTM_Zone = "21"
If Center_Lon < -60 and Center_Lon > -66 then Calculate_UTM_Zone = "20"
If Center_Lon < -66 and Center_Lon > -72 then Calculate_UTM_Zone = "19"
If Center_Lon < -72 and Center_Lon > -78 then Calculate_UTM_Zone = "18"
If Center_Lon < -78 and Center_Lon > -84 then Calculate_UTM_Zone = "17"
If Center_Lon < -84 and Center_Lon > -90 then Calculate_UTM_Zone = "16"
If Center_Lon < -90 and Center_Lon > -96 then Calculate_UTM_Zone = "15"
If Center_Lon < -96 and Center_Lon > -102 then Calculate_UTM_Zone = "14"
If Center_Lon < -102 and Center_Lon > -108 then Calculate_UTM_Zone = "13"
If Center_Lon < -108 and Center_Lon > -114 then Calculate_UTM_Zone = "12"
If Center_Lon < -114 and Center_Lon > -120 then Calculate_UTM_Zone = "11"
If Center_Lon < -120 and Center_Lon > -126 then Calculate_UTM_Zone = "10"
If Center_Lon < -126 and Center_Lon > -132 then Calculate_UTM_Zone = "09"
If Center_Lon < -132 and Center_Lon > -138 then Calculate_UTM_Zone = "08"
If Center_Lon < -138 and Center_Lon > -144 then Calculate_UTM_Zone = "07"
If Center_Lon < -144 and Center_Lon > -150 then Calculate_UTM_Zone = "06"
If Center_Lon < -150 and Center_Lon > -156 then Calculate_UTM_Zone = "05"
If Center_Lon < -156 and Center_Lon > -162 then Calculate_UTM_Zone = "04"
If Center_Lon < -162 and Center_Lon > -168 then Calculate_UTM_Zone = "03"
If Center_Lon < -168 and Center_Lon > -174 then Calculate_UTM_Zone = "02"
If Center_Lon < -174 and Center_Lon > -180 then Calculate_UTM_Zone = "01"
If Center_Lon > 0 and Center_Lon < 6 then Calculate_UTM_Zone = "31"
If Center_Lon > 6 and Center_Lon < 12 then Calculate_UTM_Zone = "32"
If Center_Lon > 12 and Center_Lon < 18 then Calculate_UTM_Zone = "33"
If Center_Lon > 18 and Center_Lon < 24 then Calculate_UTM_Zone = "34"
If Center_Lon > 24 and Center_Lon < 30 then Calculate_UTM_Zone = "35"
If Center_Lon > 30 and Center_Lon < 36 then Calculate_UTM_Zone = "36"
If Center_Lon > 36 and Center_Lon < 42 then Calculate_UTM_Zone = "37"
If Center_Lon > 42 and Center_Lon < 48 then Calculate_UTM_Zone = "38"
If Center_Lon > 48 and Center_Lon < 54 then Calculate_UTM_Zone = "39"
If Center_Lon > 54 and Center_Lon < 60 then Calculate_UTM_Zone = "40"
If Center_Lon > 60 and Center_Lon < 66 then Calculate_UTM_Zone = "41"
If Center_Lon > 66 and Center_Lon < 72 then Calculate_UTM_Zone = "42"
If Center_Lon > 72 and Center_Lon < 78 then Calculate_UTM_Zone = "43"
If Center_Lon > 78 and Center_Lon < 84 then Calculate_UTM_Zone = "44"
If Center_Lon > 84 and Center_Lon < 90 then Calculate_UTM_Zone = "45"
If Center_Lon > 90 and Center_Lon < 96 then Calculate_UTM_Zone = "46"
If Center_Lon > 96 and Center_Lon < 102 then Calculate_UTM_Zone = "47"
If Center_Lon > 102 and Center_Lon < 108 then Calculate_UTM_Zone = "48"
If Center_Lon > 108 and Center_Lon < 114 then Calculate_UTM_Zone = "49"
If Center_Lon > 114 and Center_Lon < 120 then Calculate_UTM_Zone = "50"
If Center_Lon > 120 and Center_Lon < 126 then Calculate_UTM_Zone = "51"
If Center_Lon > 126 and Center_Lon < 132 then Calculate_UTM_Zone = "52"
If Center_Lon > 132 and Center_Lon < 138 then Calculate_UTM_Zone = "53"
If Center_Lon > 138 and Center_Lon < 144 then Calculate_UTM_Zone = "54"
If Center_Lon > 144 and Center_Lon < 150 then Calculate_UTM_Zone = "55"
If Center_Lon > 150 and Center_Lon < 156 then Calculate_UTM_Zone = "56"
If Center_Lon > 156 and Center_Lon < 162 then Calculate_UTM_Zone = "57"
If Center_Lon > 162 and Center_Lon < 168 then Calculate_UTM_Zone = "58"
If Center_Lon > 168 and Center_Lon < 174 then Calculate_UTM_Zone = "59"
If Center_Lon > 174 and Center_Lon < 180 then Calculate_UTM_Zone = "60"
End Function
' 2. "cv" - The Latitude/Longitude reporting function was
' changed for HDF-EOS files. If the starting coordinate
' for the mapping between the image array and the
' Latitude/Longitude arrays is a coordinate that is outside
' the image, then the program will attempt to determine the
' mapping by examining the LatticePoint data field.
' (Previously, it assumed that the starting coordinate was
' the upper left corner of the image array, which gave
' incorrect results for ASTER Level 1A files.)
' Deprecated Command Line Argument Parser
' If left(objArgs.Item(I),2) = "-f" then strHDFname = mid(objArgs.Item(I),3)
'
' If left(objArgs.Item(I),2) = "-i" then drvPath = mid(objArgs.Item(I),3)
' ' Add a terminating delimiter if not present
' If right(drvPath,1) <> "\" then drvPath = drvPath & "\"
'
' If left(objArgs.Item(I),2) = "-o" then drvDestPath = mid(objArgs.Item(I),3)
' ' Add a terminating delimiter if not present
' If right(drvDestPath,1) <> "\" then drvDestPath = drvDestPath & "\"
'
' If left(objArgs.Item(I),2) = "-g" then GeographicType = True
' If left(objArgs.Item(I),2) = "-t" then ModelTransformation = False
' If left(objArgs.Item(I),2) = "-s" then MergeBands = False
'
' If left(objArgs.Item(I),3) = "-vo" then
' bBandGroupProcess(VNIR) = True
' bBandGroupProcess(SWIR) = False
' bBandGroupProcess(TIR) = False
' end if
More information about the Gdal-dev
mailing list