使用 python 从 .NC 文件中提取水使用数据

Question

我正在尝试使用以下网站的全球用水数据 https://zenodo.org/record/897933#.Yj4vbufMKUl。 我设法打开压缩文件“ Domestic Water use.7z”，但注意到文件 cons_dom.nc 和 withd_dom.nc 都在 .nc格式。

我对这些类型的文件相当陌生，我尝试使用以下代码将数据下载到 .tif 。然而，我似乎一直得到错误的结果。数据应该是世界的粗略覆盖，但在 arcpro 中打开 .tif 时，我看到的只是屏幕中间的一个长矩形。

另外，我很抱歉没有分享图像，但我不知道如何将 tifs 直接绘制到 python 中。

from osgeo import gdal
import glob
import os
import numpy as np
from osgeo import osr

NCfilenames = sorted(glob.glob("*.NC"))

for ncfile in NCfilenames:
    #print(ncfile)
    basencfilename = os.path.basename(ncfile)
    # print(basencfilename)
    domestic_name = basencfilename
    print(domestic_name)
    CleanBaseStr = os.path.splitext(basencfilename)[0]

    OutFileName = CleanBaseStr + "_domestic_wateruse.tif"
    print(OutFileName)

    NDVI_ds = gdal.Open(domestic_name, gdal.GA_ReadOnly)
    # print(NDVI_ds)
    width = NDVI_ds.RasterXSize
    height = NDVI_ds.RasterYSize
    NDVI_band = NDVI_ds.GetRasterBand(1)
    NDVI_arr = NDVI_band.ReadAsArray()

    gt = NDVI_ds.GetGeoTransform()

    wkt = NDVI_ds.GetProjection()
    driver = gdal.GetDriverByName("GTiff")
    out_ds = driver.Create(OutFileName, NDVI_band.XSize, NDVI_band.YSize, 1, gdal.GDT_Int16)

    # #writing output raster
    out_ds.GetRasterBand(1).WriteArray(NDVI_arr)
    out_ds.GetRasterBand(1).SetNoDataValue(-9999)

     # setting extent of output raster
     # top left x, w-e pixel resolution, rotation, top left y, rotation, n-s pixel resolution
    out_ds.SetGeoTransform(gt)

    # setting spatial reference of output raster
    srs = osr.SpatialReference()
    srs.ImportFromWkt(wkt)
    out_ds.SetProjection(srs.ExportToWkt())

print('Processing Done')

gdalinfo 用于输入之一 cons_dom.nc

Driver: netCDF/Network Common Data Format
Files: cons_dom.nc
Size is 67420, 480
Metadata:
  cons_dom#description= global gridded domestic water consumption results: 67420 grids, 480 Months
  cons_dom#units=mm/month
  month#units=months since 1971-1
Corner Coordinates:
Upper Left  (    0.0,    0.0)
Lower Left  (    0.0,  480.0)
Upper Right (67420.0,    0.0)
Lower Right (67420.0,  480.0)
Center      (33710.0,  240.0)
Band 1 Block=67420x1 Type=Float64, ColorInterp=Undefined
  NoData Value=9.96920996838686905e+36
  Unit Type: mm/month
  Metadata:
    description= global gridded domestic water consumption results: 67420 grids, 480 Months
    NETCDF_VARNAME=cons_dom
    units=mm/month

gdalinfo 用于输出 cons_dom_domestic_wateruse.tif

Driver: GTiff/GeoTIFF
Files: cons_dom_domestic_wateruse.tif
Size is 67420, 480
Origin = (0.000000000000000,0.000000000000000)
Pixel Size = (1.000000000000000,1.000000000000000)
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (   0.0000000,   0.0000000)
Lower Left  (       0.000,     480.000)
Upper Right (   67420.000,       0.000)
Lower Right (   67420.000,     480.000)
Center      (   33710.000,     240.000)
Band 1 Block=67420x1 Type=Int16, ColorInterp=Gray
  NoData Value=-9999

Answer 1

我下载并检查了你的数据集。

这根本不是空间信息：仔细阅读其描述：

cons_dom#description=global gridded domestic water consumption results: 67420 grids, 480 Months

您的栅格数据集有 67420x480 值。

这是一个时间数据集。

您仍然可以使用 GDAL 来读取它 - 但它能为您做的就是向您呈现 67420x480 值 - 由您来解释该数据。