GDAL : Reprojecting netCDF file

Question

I am trying to convert netCDF files to EPSG:3857 for use with Mapbox by using GDAL. This would be .nc to .nc conversion. Not to raster. I am open to using GDAL or other methods to do this. This data must be reprojected before it goes to a console app - and this process is taking weeks to find a solution for - I figured it was simple.

I am working on colorizing satellite data. There are 3 .nc files (blue, red, and infrared) that when combined and processed create a color image. After the 3 files are downloaded (from Amazon AWS), a python console app does the processing and dumps a .jpg to the same folder. The source code for that application is Located here so you may validate the data. (It is slow as the files are super high resolution).

The code I have tried is :

gdalwarp -t_srs EPSG:3857 test.nc test-projected.nc

However, there have been several other variations tried and nothing works.

I am not a professional with this, but should I even be using gdalwarp to do this? I only want to change the projection - nothing else, so the python app can still work with the data. It must be able to create the .jpg using the reprojected files.

The following links are samples of the data that needs to be converted :

.nc file on AWS > Color Channel 1 (Blue 1km resolution)

.nc file on AWS > Color Channel 2 (Red, Higher 0.5km resolution & larger file size)

.nc file on AWS > Color Channel 3 (Infrared - serves as green)

Additionaly, someone else online has accomplished this using a similar projection via the pyproj module at https://github.com/blaylockbk/pyBKB_v2/tree/master/BB_GOES16. (Mine must be EPSG:3857 for use with Mapbox). If the python code were modified to do this all in one go, that would be great too. I am opening a bounty as the final hope.

I do not know python, so I have been attempting GDAL for the most part- however working python code added to my source code to achieve the expected result (or a working GDAL script) will earn the bounty.

Mammy

Sadly, the Github link is dead – Aggrieve

It should be corrected now. – Mammy

Answer 1

Here is my solution:

# -*- coding: utf-8 -*-
"""
Created on Mon Mar  4 17:39:45 2019

@author: Guy Serbin
"""

import os, sys, glob, argparse
from osgeo import gdal, osr
from scipy.misc import imresize

parser = argparse.ArgumentParser(description = 'Script to create CONUS true color image from GOES 16 L1b data.')
parser.add_argument('-i', '--indir', type = str, default = r'C:\Data\Freelancer\DavidHolcomb', help = 'Input directory name.')
parser.add_argument('-o', '--outdir', type = str, default = None, help = 'Output directory name.')
parser.add_argument('-p', '--proj', type = int, default = 3857, help = 'Output projection, must be EPSG number.')
args = parser.parse_args()

if not args.indir:
    print('ERROR: --indir not set. exiting.')
    sys.exit()
elif not os.path.isdir(args.indir):
    print('ERROR: --indir not set to a valid directory path. exiting.')
    sys.exit()

if not args.outdir:
    print('WARNING: --outdir not set. Output will be written to --indir.')
    args.outdir = args.indir

o_srs = osr.SpatialReference()
o_srs.ImportFromEPSG(args.proj)


# based upon code ripped from https://riptutorial.com/gdal/example/25859/read-a-netcdf-file---nc--with-python-gdal

# Path of netCDF file
netcdf_red = glob.glob(os.path.join(args.indir, 'OR_ABI-L1b-RadC-M3C02_G16_s*.nc'))[0]
netcdf_green = glob.glob(os.path.join(args.indir, 'OR_ABI-L1b-RadC-M3C03_G16_s*.nc'))[0]
netcdf_blue = glob.glob(os.path.join(args.indir, 'OR_ABI-L1b-RadC-M3C01_G16_s*.nc'))[0]
baselist = os.path.basename(netcdf_blue).split('_')

outputfilename = os.path.join(args.outdir, 'OR_ABI-L1b-RadC-M3TrueColor_1_G16_{}.tif'.format(baselist[3]))
print('Output file will be: {}'.format(outputfilename))
tempfile = os.path.join(args.outdir, 'temp.tif')

# Specify the layer name to read
layer_name = "Rad"

# Open netcdf file.nc with gdal
print('Opening red band file: {}'.format(netcdf_red))
dsR = gdal.Open("NETCDF:{0}:{1}".format(netcdf_red, layer_name))
print('Opening green band file: {}'.format(netcdf_green))
dsG = gdal.Open("NETCDF:{0}:{1}".format(netcdf_green, layer_name))
print('Opening blue band file: {}'.format(netcdf_blue))
dsB = gdal.Open("NETCDF:{0}:{1}".format(netcdf_blue, layer_name))
red_srs = osr.SpatialReference()
red_srs.ImportFromWkt(dsR.GetProjectionRef())
i_srs = osr.SpatialReference()
i_srs.ImportFromWkt(dsG.GetProjectionRef())
GeoT = dsG.GetGeoTransform()
print(i_srs.ExportToWkt())
red_transform = osr.CoordinateTransformation(red_srs, o_srs)
transform = osr.CoordinateTransformation(i_srs, o_srs)

# Read full data from netcdf

print('Reading red band into memory.')
red = dsR.ReadAsArray(0, 0, dsR.RasterXSize, dsR.RasterYSize)
print('Resizing red band to match green and blue bands.')
red = imresize(red, 50, interp = 'bicubic')
print('Reading green band into memory.')
green = dsG.ReadAsArray(0, 0, dsG.RasterXSize, dsG.RasterYSize)
print('Reading blue band into memory.')
blue = dsB.ReadAsArray(0, 0, dsB.RasterXSize, dsB.RasterYSize)
red[red < 0] = 0
green[green < 0] = 0
blue[blue < 0] = 0

# Stack data and output
print('Stacking data.')
driver = gdal.GetDriverByName('GTiff')
stack = driver.Create('/vsimem/stack.tif', dsB.RasterXSize, dsB.RasterYSize, 3, gdal.GDT_Int16)
stack.SetProjection(i_srs.ExportToWkt())
stack.SetGeoTransform(GeoT)
stack.GetRasterBand(1).WriteArray(red)
stack.GetRasterBand(2).WriteArray(green)
stack.GetRasterBand(3).WriteArray(blue)
print('Warping data to new projection.')
warped = gdal.Warp('/vsimem/warped.tif', stack, dstSRS = o_srs, outputType = gdal.GDT_Int16)

print('Writing output to disk.')

outRaster = gdal.Translate(outputfilename, '/vsimem/warped.tif')

outRaster = None
red = None
green = None
blue = None
tmp_ds = None
dsR = None
dsG = None
dsB = None

print('Processing complete.')

Runin

Dear Guy Serbin, thank you for this great solution of yours. One thing it got me wondering was how one can manage other dimensions in this WARPING convertion. For instance, there are several cases in which a NETCDF file has dimensions like TIME, if not even other ones (like altitude from the Ground). As it is now, the script above will mishandle these other dimensions within the NETCDF. Would you know a good approach to handle all dimensions within a NETCDF during a spatial reprojection? – Feather

Answer 2

You can use rioxarray to do this. An example of doing so is here: https://corteva.github.io/rioxarray/html/examples/reproject.html

Here is an example targeted for your use case:

import rioxarray

xds = rioxarray.open_rasterio("OR_ABI-L1b-RadC-M3C01_G16_s20190621802131_e20190621804504_c20190621804546.nc")

<xarray.Dataset>
Dimensions:      (band: 1, x: 5000, y: 3000)
Coordinates:
  * y            (y) float64 1.584e+06 1.585e+06 ... 4.588e+06 4.589e+06
  * x            (x) float64 -3.627e+06 -3.626e+06 ... 1.381e+06 1.382e+06
  * band         (band) int64 1
    spatial_ref  int64 0
Data variables:
    Rad          (band, y, x) int16 ...
    DQF          (band, y, x) int8 ...

xds.rio.crs

CRS.from_wkt('PROJCS["unnamed",GEOGCS["unknown",DATUM["unnamed",SPHEROID["Spheroid",6378137,298.2572221]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]]],PROJECTION["Geostationary_Satellite"],PARAMETER["central_meridian",-75],PARAMETER["satellite_height",35786023],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],EXTENSION["PROJ4","+proj=geos +lon_0=-75 +h=35786023 +x_0=0 +y_0=0 +ellps=GRS80 +units=m +no_defs +sweep=x"]]')

Then, reproject:

xds_3857 = xds.rio.reproject("epsg:3857")

<xarray.Dataset>
Dimensions:      (band: 1, x: 7693, y: 4242)
Coordinates:
  * x            (x) float64 -1.691e+07 -1.691e+07 ... -5.892e+06 -5.891e+06
  * y            (y) float64 7.714e+06 7.712e+06 ... 1.641e+06 1.64e+06
  * band         (band) int64 1
    spatial_ref  int64 0
Data variables:
    Rad          (band, y, x) int16 1023 1023 1023 1023 ... 1023 1023 1023 1023
    DQF          (band, y, x) int8 0 0 0 0 0 0 0 0 0 0 0 ... 0 0 0 0 0 0 0 0 0 0
Attributes:
    creation_date:  2019-09-25 01:02:54.590053

xds_3857.rio.crs

CRS.from_epsg(3857)

Write to netcdf:

xds_3857.to_netcdf("epsg3857.nc")

Reeve