Skip to content

Instantly share code, notes, and snippets.

@robintw
Created July 21, 2016 15:55
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save robintw/fccd3fcb084dbe78c51d5b1833760f80 to your computer and use it in GitHub Desktop.
Save robintw/fccd3fcb084dbe78c51d5b1833760f80 to your computer and use it in GitHub Desktop.
import os
import gdal
import glob
# # Translate Lat and Lon files
# - Find out information about the lat lon file
# - Using the gdal_translate tool convert the lat and lon hdr files into vrt format
# gdal.Translate(r"h00v01_lat.vrt",
# r'HDF4_EOS:EOS_GRID:"MAIACLatlon.h00v01.hdf":latlon:lat', format='VRT')
#
# gdal.Translate(r"h00v01_lon.vrt",
# r'HDF4_EOS:EOS_GRID:"MAIACLatlon.h00v01.hdf":latlon:lon', format='VRT')
# # Create geoloc.vrt
geoloc_xml = """<VRTDataset rasterXSize="1200" rasterYSize="1200">
<Metadata domain="GEOLOCATION">
<mdi key="X_DATASET">h00v01_lon.vrt</mdi>
<mdi key="X_BAND">1</mdi>
<mdi key="Y_DATASET">h00v01_lat.vrt</mdi>
<mdi key="Y_BAND">1</mdi>
<mdi key="PIXEL_OFFSET">0</mdi>
<mdi key="LINE_OFFSET">0</mdi>
<mdi key="PIXEL_STEP">1</mdi>
<mdi key="LINE_STEP">1</mdi>
</Metadata>
<VRTRasterBand dataType="Float32" band="1">
<ColorInterp>Gray</ColorInterp>
<SimpleSource>
<SourceFilename relativeToVRT="0">HDF4_EOS:EOS_GRID:"{fname}":grid1km:Optical_Depth_Land</SourceFilename>
<SourceBand>1</SourceBand>
<SourceProperties RasterXSize="1200" RasterYSize="1200" DataType="Float32" BlockXSize="1200" BlockYSize="1200" />
<SrcRect xOff="0" yOff="0" xSize="1200" ySize="1200" />
<DstRect xOff="0" yOff="0" xSize="1200" ySize="1200" />
</SimpleSource>
</VRTRasterBand>
</VRTDataset>
"""
def geolocate_and_warp(input_filename, output_filename, output_crs='EPSG:27700'):
"""Takes a MAIAC AOT file and performs geolocation using
the associated MAIAC Lat Lon files, and then warps to
the specified co-ordinate system.
Arguments:
- `input_filename`:
The input MAIAC AOTfile, in HDF format.
- `output_filename:
The file to store the output in. Will be a GeoTIFF file
so should have a .tif extension
- `output_crs`: (optional, default: EPSG:27700 for OSGB)
The output Co-ordinate Reference System (CRS), specified
as a string such as 'EPSG:27700'.
Notes:
- This currently assumes that you've created lat and lon VRTs
already, and has only been tested on h00v01 data for Europe
"""
with open(r'geoloc_test.vrt', 'w') as myfile:
myfile.write(geoloc_xml.format(fname=input_filename))
gdal.Warp(output_filename,
r'geoloc_test.vrt',
geoloc=True, srcSRS='EPSG:4326', dstSRS=output_crs,
srcNodata=-28672, dstNodata=-28672)
def preprocess_MAIAC(input_folder, output_folder, glob_string):
# Get the list of all files in the input folder that
# match the glob string
filenames = glob.glob(os.path.join(input_folder, glob_string))
for filename in filenames:
print('Processing {fname}'.format(fname=filename))
# Set the output 'basename' (the base for all the other filenames)
# based on the input filename and output folder - so that the folder
# structure is replicated in output_folder.
output_basename = os.path.join(output_folder,
os.path.relpath(filename,
start=input_folder))
# Make sure the output folder exists before we try to write to it
output_dir = os.path.dirname(output_basename)
if not os.path.exists(output_dir):
os.makedirs(output_dir)
# Actually do the geolocation and warping!
geolocate_and_warp(filename, output_basename + "_proj.tif")
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment