GDAL：重新投影 netCDF 文件

Question

我正在尝试使用 GDAL 将 netCDF 文件转换为 EPSG:3857 以与 Mapbox 一起使用。这将是 .nc 到 .nc 的转换。不光栅。我愿意使用 GDAL 或其他方法来做到这一点。这些数据在进入控制台应用程序之前必须重新投影 - 这个过程需要数周时间才能找到解决方案 - 我认为这很简单。

我正在为卫星数据着色。有 3 个 .nc 文件（蓝色、红色和红外线），它们在组合和处理后会创建一个彩色图像。下载 3 个文件（从 Amazon AWS）后，python 控制台应用程序进行处理并将 .jpg 转储到同一文件夹。该应用程序的源代码位于此处，因此您可以验证数据。（由于文件是超高分辨率的，所以速度很慢）。

我试过的代码是：

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

但是，已经尝试了其他几种变体，但没有任何效果。

我不是这方面的专业人士，但我什至应该使用 gdalwarp 来做到这一点吗？我只想改变投影 - 没有别的，所以 python 应用程序仍然可以处理数据。它必须能够使用重新投影的文件创建 .jpg。

Answer 1

这是我的解决方案：# -*- coding: utf-8 -*-"""Created on Mon Mar&nbsp; 4 17:39:45 2019@author: Guy Serbin"""import os, sys, glob, argparsefrom osgeo import gdal, osrfrom scipy.misc import imresizeparser = 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:&nbsp; &nbsp; print('ERROR: --indir not set. exiting.')&nbsp; &nbsp; sys.exit()elif not os.path.isdir(args.indir):&nbsp; &nbsp; print('ERROR: --indir not set to a valid directory path. exiting.')&nbsp; &nbsp; sys.exit()if not args.outdir:&nbsp; &nbsp; print('WARNING: --outdir not set. Output will be written to --indir.')&nbsp; &nbsp; args.outdir = args.indiro_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 filenetcdf_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 readlayer_name = "Rad"# Open netcdf file.nc with gdalprint('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 netcdfprint('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] = 0green[green < 0] = 0blue[blue < 0] = 0# Stack data and outputprint('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 = Nonered = Nonegreen = Noneblue = Nonetmp_ds = NonedsR = NonedsG = NonedsB = Noneprint('Processing complete.')

Answer 2

您可以使用 rioxarray 来执行此操作。这样做的一个例子是：https : //corteva.github.io/rioxarray/html/examples/reproject.html以下是针对您的用例的示例：import rioxarrayxds = rioxarray.open_rasterio("OR_ABI-L1b-RadC-M3C01_G16_s20190621802131_e20190621804504_c20190621804546.nc")<xarray.Dataset>Dimensions:&nbsp; &nbsp; &nbsp; (band: 1, x: 5000, y: 3000)Coordinates:&nbsp; * y&nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; (y) float64 1.584e+06 1.585e+06 ... 4.588e+06 4.589e+06&nbsp; * x&nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; (x) float64 -3.627e+06 -3.626e+06 ... 1.381e+06 1.382e+06&nbsp; * band&nbsp; &nbsp; &nbsp; &nbsp; &nbsp;(band) int64 1&nbsp; &nbsp; spatial_ref&nbsp; int64 0Data variables:&nbsp; &nbsp; Rad&nbsp; &nbsp; &nbsp; &nbsp; &nbsp; (band, y, x) int16 ...&nbsp; &nbsp; DQF&nbsp; &nbsp; &nbsp; &nbsp; &nbsp; (band, y, x) int8 ...xds.rio.crsCRS.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"]]')然后，重新投影：xds_3857 = xds.rio.reproject("epsg:3857")<xarray.Dataset>Dimensions:&nbsp; &nbsp; &nbsp; (band: 1, x: 7693, y: 4242)Coordinates:&nbsp; * x&nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; (x) float64 -1.691e+07 -1.691e+07 ... -5.892e+06 -5.891e+06&nbsp; * y&nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; (y) float64 7.714e+06 7.712e+06 ... 1.641e+06 1.64e+06&nbsp; * band&nbsp; &nbsp; &nbsp; &nbsp; &nbsp;(band) int64 1&nbsp; &nbsp; spatial_ref&nbsp; int64 0Data variables:&nbsp; &nbsp; Rad&nbsp; &nbsp; &nbsp; &nbsp; &nbsp; (band, y, x) int16 1023 1023 1023 1023 ... 1023 1023 1023 1023&nbsp; &nbsp; DQF&nbsp; &nbsp; &nbsp; &nbsp; &nbsp; (band, y, x) int8 0 0 0 0 0 0 0 0 0 0 0 ... 0 0 0 0 0 0 0 0 0 0Attributes:&nbsp; &nbsp; creation_date:&nbsp; 2019-09-25 01:02:54.590053xds_3857.rio.crsCRS.from_epsg(3857)写入 netcdf：xds_3857.to_netcdf("epsg3857.nc")