Reproject NEODAAS NetCDF

reproject_neodaas_netcdf.py

#!/usr/bin/env python

"""
Create a single layer tif in WGS84 projection from
NEODAAS (https://www.neodaas.ac.uk/) netCDF files in mercator projection.

Converts the location array to GCPs and uses these to warp the image with
GDAL.

Dan Clewley (dac@pml.ac.uk).
2018-11-12

"""

import argparse
import os
import shutil
import tempfile

import netCDF4
from osgeo import gdal

# The geolocation grid in the netCDF file contains a location for
# each pixel. This is often more than GDAL will handle as GCPs so set
# a resampling factor to reduce.
GEOLOCATION_GRID_RESAMPLING = 50

def create_single_var_tiff(input_netcdf, output_tif, variable):
    """
    Export netCDF as a single variable GeoTiff
    """
    # Create single variable tif from input netCDF
    in_ds = gdal.Open('NETCDF:"{}":{}'.format(input_netcdf, variable), gdal.GA_ReadOnly)
    gdal.Translate(output_tif, in_ds)
    in_ds = None

if __name__ == "__main__":

    parser = argparse.ArgumentParser(description="Reproject NetCDF")
    parser.add_argument("innetcdf", nargs=1, type=str,
                        help="Input netCDF file")
    parser.add_argument("outfile", nargs=1, type=str,
                        help="Output file (.tif)")
    parser.add_argument("--variable", required=True,
                        help="Variable from netCDF to extract (e.g., CHL_OC4ME)")
    args = parser.parse_args()

    # Make temporary directory
    temp_dir = tempfile.mkdtemp(prefix="netcdf_convert_")
    temp_tif = os.path.join(temp_dir, os.path.basename(args.innetcdf[0]))

    # Export a single layer as a geotiff
    create_single_var_tiff(args.innetcdf[0], temp_tif, args.variable)

    # Open new geotiff to add GCPs to
    data = gdal.Open(temp_tif, gdal.GA_Update)

    # Open input netCDF
    data_nc = netCDF4.Dataset(args.innetcdf[0])

    gcp_list = []

    # Go through longitude and latitude arrays and convert to GCPs
    for x in range(1, data_nc.variables['longitude'].size, GEOLOCATION_GRID_RESAMPLING):
        for y in range(1, data_nc.variables['latitude'].size, GEOLOCATION_GRID_RESAMPLING):
            gcp_list.append(gdal.GCP(float(data_nc.variables['longitude'][x]),
                                     float(data_nc.variables['latitude'][y]), 0,
                                     x, y))

    # Add GCPs to tif. Set projection is None to ignore it
    data.SetGCPs(gcp_list, None)

    # Warp file based on GCPs
    gdal.Warp(args.outfile[0], data, dstSRS='EPSG:4326', format='GTiff')

    # Close files
    data = None
    data_nc = None

    # Remove temporary directory
    shutil.rmtree(temp_dir)