Created
July 21, 2016 15:55
-
-
Save robintw/fccd3fcb084dbe78c51d5b1833760f80 to your computer and use it in GitHub Desktop.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
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