Tuesday, February 21, 2006

What_nick's PHP script for dishing out BILs...

Tisham (What_Nick, see his blog at http://whatnick.blogspot.com) of Apogee (http://www.apogee.com.au) sent me his PHP script for dishing out BILs behind a web server or a geo-trawler. It will take a source TIF, use GDAL to cut up the source elevation data into WW tiles (and corresponding directory structure), and then GZip the files for delivery.


$szMapCacheDir="G:/CFSImagery/wwcache/DEM";
/* create the main cache directory if necessary */
if (!@is_dir($szMapCacheDir))
makeDirs($szMapCacheDir);

/* get the various request parameters
* also need to make sure inputs are clean, especially those used to
* build paths and filenames
*/
$X = isset( $_REQUEST['X'] ) ? intval($_REQUEST['X']) : 0;
$Y = isset( $_REQUEST['Y'] ) ? intval($_REQUEST['Y']) : 0;
$L = isset( $_REQUEST['L'] ) ? intval($_REQUEST['L']) : 0;
$T = isset( $_REQUEST['T'] ) ? intval($_REQUEST['T']) : 103;
$szExt = ".7z";

$szLevelcache = $szMapCacheDir."/$L";
$szYcache = sprintf($szLevelcache."/%04d",$Y);
if (!@is_dir($szYcache))
makeDirs($szYcache);
$szCacheFile = sprintf($szYcache."/%04d_%04d".$szExt,$Y,$X);
$szIntFile = sprintf($szYcache."/%04d_%04d.bil",$Y,$X);;
/* Hit Test the cache tile */
if (!file_exists($szCacheFile) || $bForce){
$lzts = 1.0;
//Our Layer Zero Tile Size

$lat1 = ($Y*$lzts*pow(0.5,$L))-90;
$lon1 = ($X*$lzts*pow(0.5,$L))-180;
//Lat2 and Lon2 as figured from tile size and level
$lat2 = $lat1 + $lzts*pow(0.5,$L);
$lon2 = $lon1 + $lzts*pow(0.5,$L);
if($T==103){
if(($lat1>-33)||($lon1<138)||($lat2<-36)||($lon2>140)){
header("HTTP/1.0 404 Not Found");
exit();
}
else{
$gdalwarp = "gdalwarp.exe -te $lon1 $lat1 $lon2 $lat2 -ot Int16 -ts 150 150 -of ENVI ".
"G:/CFSImagery/latlong/dem/hillsdem.tif ".$szIntFile;
exec($gdalwarp);
$za7 = "7za.exe a ".$szCacheFile." ".$szIntFile;
exec($za7);
}
}
if($T==104){
if(($lat1>28)||($lon1<9)||($lat2<27)||($lon2>11)){
header("HTTP/1.0 404 Not Found");
exit();
}
else{
$gdalwarp = "gdalwarp.exe -te $lon1 $lat1 $lon2 $lat2 -ot Int16 -ts 250 250 -of ENVI ".
"G:/CFSImagery/latlong/dem/LibyaDuneslatlonint ".$szIntFile;
exec($gdalwarp);
$za7 = "7za.exe a ".$szCacheFile." ".$szIntFile;
exec($za7);
}
}
}
$h = fopen($szCacheFile, "r");
header("Content-Type: "."application/x-7z-compressed");
header("Content-Length: " . filesize($szCacheFile));
header("Expires: " . date( "D, d M Y H:i:s GMT", time() + 31536000 ));
header("Cache-Control: max-age=31536000, must-revalidate" );
fpassthru($h);
fclose($h);

Viewing BIL files...

You've got a BIL, but want to look at the data. Here's how:
  1. Make sure its not 7zip'd -- if it is, then unzip it. Usually you do not have to worry about this, as WW unzips them for you.
  2. Make a copy of the file, and rename it to .raw
  3. Figure out the dimensions of the file: you can do this by math or by looking at the .XML description. For instance, the SRTM30 has the following description inside C:\Program Files\NASA\World Wind 1.3\Config\Earth.xml

SamplesPerTile: 150
DataFormat: Int16
FileExtension: BIL
CompressonType: 7z

This says that each tile is 150 x 150 pixels, with 16 bit integers. To use math, just take the file size of the uncompressed BIL (e.g. 45000 bytes for SRTM30 data), then divide by 2 (2 bytes per pixel), then square root (assuming all BILs are square). (45000 / 2) ^ 0.5 = 150 pixels.

Now, choose your poison to view the BILs. For quickness, I use IrfanView:

  1. Load the .RAW file copy into Irfan View.
  2. Set the RAW Parameters:
    The image width and height are calculated in step 3 above
    The file header is 0 bytes
    The BIL is 16 bits per pixel (or different, as spec'd in the XML file)
  3. The color will give you only a rough idea, as it is using three channels to display a monochrome color.

To view it in Photoshop:

  1. Open the .RAW file in Photoshop
  2. Set the RAW Parameters:
    The image width and height are calculated in step 3 above
    The number of channels is 1
    The Depth is 16 bits per pixel (or different, as spec'd in the XML file)
    The Byte order is IBM.
    The file header is 0 bytes
  3. Pull up the levels adjustment (Ctrl-L)
    Move the gamma (center grey) tab to see the intermediate height values.

More information:

http://www.worldwindcentral.com/wiki/Add-on:Colored_SRTM_elevation_data_add-on
http://forum.worldwind.arc.nasa.gov/lofiversion/index.php?t5508.html
http://support.intergraph.com/Geospatial/Downloads/Tools.asp?ID=70&SORT=Title


Note: the referenced TIF2BIL tool from Integraph above requires zifl00.dll -- which comes with INGR's I/RasC and ImageViewer. It also may be available in their raster utilities (thanks Mazop!). Also:from Pangloss: actually it's bil2tif not tif2bil; ISRUStandalone contains zirfl01.dll not zirfl00.dll.