GDAL: перепроектирование файла netCDF

Question

Я пытаюсь преобразовать файлы netCDF в EPSG: 3857 для использования с Mapbox с помощью GDAL.Это будет преобразование .nc в .nc.Не в растр.Я открыт для использования GDAL или других методов для этого.Эти данные должны быть перепроецированы до того, как они попадут в консольное приложение - и этот процесс занимает недели, чтобы найти решение - я подумал, что это было просто.

Я работаю над раскрашиванием спутниковых данных.Существует 3 файла .nc (синий, красный и инфракрасный), которые при объединении и обработке создают цветное изображение.После загрузки трех файлов (из Amazon AWS) консольное приложение python выполняет обработку и выгружает файл .jpg в ту же папку.Исходный код для этого приложения Расположен здесь, чтобы вы могли проверить данные .(Это медленный файл, поскольку файлы имеют сверхвысокое разрешение.)

Код, который я пробовал:

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

Однако было предпринято несколько других вариантов, и ничего не работает.

Я не профессионал с этим, но я должен даже использовать gdalwarp, чтобы сделать это?Я только хочу изменить проекцию - больше ничего, поэтому приложение python все еще может работать с данными.Он должен иметь возможность создавать .jpg с помощью перепроецированных файлов.

Следующие ссылки представляют собой примеры данных, которые необходимо преобразовать:

.nc файл в AWS> ColorКанал 1 (синий, разрешение 1 км)

.nc файл на AWS> Цвет Канал 2 (красный, более высокое разрешение 0,5 км и больший размер файла)

.nc-файл на AWS> Цветовой канал 3 (инфракрасный - служит зеленым)

Кроме того, кто-то еще в сети сделал это, используя аналогичную проекцию через модуль pyproj на https://github.com/blaylockbk/pyBKB_v2/tree/master/BB_GOES16. (Mine mustбыть EPSG: 3857 для использования с Mapbox).Если бы код Python был изменен, чтобы сделать все это за один раз, это тоже было бы здорово.Я открываю щедрость как последнюю надежду.

Я не знаю Python, поэтому я пытался GDAL по большей части - однако рабочий код Python добавленк моему исходному коду для достижения ожидаемого результата (или работающий скрипт GDAL) получит награду.

Answer 1

Вот мое решение:

# -*- 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.')