Extract Water Use data from .NC file using python

Question

I am trying to use global wateruse data from the following website https://zenodo.org/record/897933#.Yj4vbufMKUl. I managed to open the zipped file, domestic water use.7z, but noticed that the both files cons_dom.nc and withd_dom.ncare in .nc format.

I am fairly new with these types of files and I tried downloading the data to .tif using the following code. However it seems that I keep on getting a wrong result. The data should be a coarse cover of the world, but all I am seeing is a long rectangle in the middle of the screen when opening the .tif in arcpro.

Also, I am sorry for not sharing an image, but I do not know how to plot the tifs directly into 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 for one of the inputs 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 for the output 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

@mmomtchev just edited the question, both `gdalinfo` for input and output are added — Jay Antonio Oliver, Mar 28 '22 at 15:49
The input file is not properly georeferenced and lacks a `geoTransform` — mmomtchev, Mar 28 '22 at 16:20
@mmomtchev these are the 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)` I do not seem to be able to find the projection — Jay Antonio Oliver, Mar 28 '22 at 17:19
Ask whoever provided you this file about the projection and its geographic corner coordinates - these are the raster set coordinates — mmomtchev, Mar 28 '22 at 17:24

score 0 · Answer 1 · answered Mar 29 '22 at 00:18

0

I downloaded and examined your dataset.

This is not spatial information at all: read carefully its description:

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

Your raster dataset has 67420x480 values.

This is a temporal dataset.

You can still use GDAL to read it - but everything it can do for you is to present you with the 67420x480 values - it is up to you to interpret that data.

answered Mar 29 '22 at 00:18

mmomtchev

2,497
1
8
23

Oh, well thank you for the clarification. So in that case, I should be able to convert this data into a CSV, correct? – Jay Antonio Oliver Mar 29 '22 at 00:32
Yes, CSV would be a good format for it – mmomtchev Mar 29 '22 at 08:31

Extract Water Use data from .NC file using python

1 Answers1