python - GDAL : Reprojecting netCDF file

Question

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

我正在为卫星数据着色。有 3 个 .nc 文件(蓝色、红色和红外线)在组合和处理时会创建彩色图像。在(从 Amazon AWS)下载 3 个文件后，python 控制台应用程序进行处理并将 .jpg 转储到同一文件夹。该应用程序的源代码是 Located here so you may validate the data . (因为文件是超高分辨率，所以速度很慢)。

我试过的代码是:

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

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

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

以下链接是需要转换的数据示例:

.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)

另外，网上的其他人已经通过 https://github.com/blaylockbk/pyBKB_v2/tree/master/BB_GOES16 的 pyproj 模块使用类似的投影来完成此操作。 . (我的必须是 EPSG:3857 才能与 Mapbox 一起使用)。如果修改 python 代码以一次性完成所有操作，那也很棒。我打开赏金作为最后的希望。



我不知道python，所以我大部分时间都在尝试GDAL-但是添加到我的源代码中以实现预期结果(或工作GDAL脚本)的工作python代码将获得赏金。

Answer 1

您可以使用 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")