[Gdal-dev] RE: [osg-users] Next gen blue marble imagery

[Norman Vine]

Burns, Christopher K writes:
> 
> I have downloaded the full image: world.200410.3x86400x43200.bin, and
> although it is simply a raw image gdal_translate does not seem to be
> accepting it.  Could someone suggest the proper command line to enter to
> convert this to a GeoTIFF?

attached find a python script that will build geotiffs from your file

note this is hardwired to "world.200410.3x86400x43200.bin" 
but this is easily modified 

see doc string in script for details

Norman

-------------- next part --------------
#!/usr/bin/env python

"""
Create tiled geotiff with overviews from the Next Gen Blue Marble datasets
we work with the full image and break it up into 8 tiles using the original
NASA Tiling and Naming scheme

requires the GDAL utilities gdaltranslate and gdaladdo
GDAL available from http://www.gdal.org
for .vrt format info see http://www.gdal.org/gdal_vrttut.html

data files available from  http://bmng.telascience.org/world_big/
note the files on above site are gzipped
this script expects the file to be decompressed

author  Norman Vine nhv at cape.com
19-Dec-2005 10:03:09 am
"""

import os

def do_cmd(cmd):    
    print cmd
    os.system(cmd)
    
def make_vrt(name):
    vrtStr =  '<VRTDataset rasterXSize=\"86400\" rasterYSize=\"43200\">\n'
    vrtStr += ' <SRS>WGS84</SRS>\n'
    vrtStr += ' <GeoTransform>-180.0,0.004166666666666,0,90.0,0,-0.004166666666666</GeoTransform>\n'

    for i in range(3):
        vrtStr += ' <VRTRasterBand dataType=\"Byte\" band=\"%d\" subClass="VRTRawRasterBand">\n'%(i+1)
        vrtStr += '  <SourceFilename relativetoVRT=\"1\">%s</SourceFilename>\n'%(name)
        vrtStr += '  <ImageOffset>%d</ImageOffset>\n'%(i)
        vrtStr += '  <PixelOffset>3</PixelOffset>\n'
        vrtStr += '  <LineOffset>259200</LineOffset>\n'
        vrtStr += ' </VRTRasterBand>\n'
    vrtStr += '</VRTDataset>\n'
    #print vrtStr
    #return

    base = name[:-4]
    vrtName = base+".vrt"

    vrtFile = open(vrtName,"w")
    vrtFile.write(vrtStr)
    vrtFile.close()

def make_tiffs(name):
    trans_extent = [ "-180.0 90.0 -90.0 0.0",
                     " -90.0 90.0   0.0 0.0",
                     "   0.0 90.0  90.0 0.0",
                     "  90.0 90.0 180.0 0.0",
                     "-180.0  0.0 -90.0 -90.0",
                     " -90.0  0.0   0.0 -90.0",
                     "   0.0  0.0  90.0 -90.0",
                     "  90.0  0.0 180.0 -90.0" ]
                    
    section = ["A1","B1","C1","D1","A2","B2","C2","D2"]

    base = name[:-4]
    vrtName = base+".vrt"
    
    trans_cmd = 'gdal_translate -of GTiff -co \"TILED=YES\" -co \"BLOCKXSIZE=128\" -co \"BLOCKYSIZE=128\" -projwin'
    oview_cmd = 'gdaladdo'
    for i in range(8):
        outfile = "%s.%s.tif"%(base,section[i]);
        do_cmd( '%s %s %s %s'%(trans_cmd,trans_extent[i],vrtName,outfile) )
        do_cmd( '%s %s 2 4 8 16 32 64 128 256'%(oview_cmd,outfile) )


def process(name):
    make_vrt(name)
    make_tiffs(name)

# hardwired test    
if __name__ == "__main__":
    process("world.200410.3x86400x43200.bin")