[Gdal-dev] RE: [osg-users] Next gen blue marble imagery
Norman Vine
nhv at cape.com
Mon Dec 19 10:42:34 EST 2005
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")
More information about the Gdal-dev
mailing list