我正在嘗試使用 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。
2 回答
慕婉清6462132
TA貢獻1804條經驗 獲得超2個贊
這是我的解決方案：
# -*- 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.')
斯蒂芬大帝
TA貢獻1827條經驗 獲得超8個贊
您可以使用 rioxarray 來執行此操作。這樣做的一個例子是：https : //corteva.github.io/rioxarray/html/examples/reproject.html
以下是針對您的用例的示例：
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"]]')
然后，重新投影：
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)
寫入 netcdf：
xds_3857.to_netcdf("epsg3857.nc")
添加回答
舉報
0/150
提交
取消