0

I'm trying to loop through and read a stack of ESRI grids in a specific projection, convert them to GDAL and then write them out to a NetCDF, but I've never performed this from scratch and I'm running into problems understanding some conventions/processes.

My approach is:

  1. glob my ascii files and perform a subprocess.call of gdalwarp to reproject/reformat/rename my source grids

  2. For each group of variables (e.g. precip, tasmin, tasmax), open a new NetCDF container, define some global attributes, and loop through the files.

  3. Open my files to be added with gdal.Open(), create variables on the first run, then read in my files as gdal arrays.

My issue is that I can't seem to succesfully write out my array to NetCDF. AFAIK, the lats, lons, and variable values are properly being written to objects, but when these objects are read in to the NetCDF variables, I keep getting consecutive values on my axes (e.g. lon = [-180, -180 .. -180 -180]). I'm not sure what I might be missing. Code below:

#!/bin/python3

import os
import glob
import subprocess
import netCDF4 as nc
import numpy as np
from osgeo import gdal
import datetime

#from pdb import set_trace as stop

# Enable GDAL/OGR exceptions
gdal.UseExceptions()

path_source = "~/Documents/Script_Examples/Practice/"
path_destination = "~/Documents/Script_Examples/Practice/"

CNAprj = "+proj=lcc +lat_1=49 +lat_2=77 +lat_0=0 +lon_0=-95 +x_0=0 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs"
WGS84prj = "+proj=longlat +datum=WGS84 +no_defs"
output = "GTiff"

file_location = os.path.join(path_source, "**", "*.asc")
files = glob.glob(file_location, recursive=True)

tifs = []
for f in files:
    print(f)
    tif = '.'.join([f.split(".")[0], "rpj.tif"])
    subprocess.call(['gdalwarp', f, tif, "-s_srs", CNAprj, "-t_srs", WGS84prj, "-of", output])
    tifs.append(tif)
tifs.sort()
print(tifs)

# Dictionaries and lists of climate variables handled
mon_variables = {"PPT": "pr", "Tave": "tasmean", "Tmax": "tasmax", "Tmin": "tasmin"}
bio_variables = ['AHM', 'bFFP', 'CMD', 'DD5', 'DD18', 'DD_0', 'DD_18', 'e_FFP', 'EMT', 'Eref', 'EXT',\
                 'FFP', 'MAP', 'MAR', 'MAT', 'MCMT', 'MSP', 'MWMT', 'NFFD', 'PAS', 'PPT_sm', 'PPT_wt',\
                 'RH', 'SHM', 'Tave_sm', 'Tave_wt', 'TD']

Conventions = "CF-1.7"
institution = "University of Alberta"
contact = "Andreas Hamann (andreas.hamann@ualberta.ca)"
comment = "30-year climatologies based on CMIP5 model ensembles and downscaled using the PRISM technique"
title = "ClimateNA"
project = "Historical and projected climate data for North America"
references = "Wang, T., Hamann, A. Spittlehouse, D.L. and Carroll, C. 2016. Locally downscaled and spatially customizable climate data for historical and future periods for North America. PLoS One 11: e0156720"

var_metadata = {
    "pr": {
        "long_name": "Precipitation",
        "standard_name": "precipitation",
        "units": "mm",
        "coordinates": "lon lat time",
        "cell_methods": "time: mean"
    },
    "tasmean": {
        "long_name": "Mean Air Temperature",
        "standard_name": "mean_air_temperature",
        "units": "deg C",
        "coordinates": "lon lat time",
        "cell_methods": "time: mean"
    },
    "tasmin": {
        "long_name": "Minimum Air Temperature",
        "standard_name": "min_air_temperature",
        "units": "deg C",
        "coordinates": "lon lat time",
        "cell_methods": "time: mean"
    },
    "tasmax": {
        "long_name": "Maximum Air Temperature",
        "standard_name": "max_air_temperature",
        "units": "deg C",
        "coordinates": "lon lat time",
        "cell_methods": "time: mean"
    }
}

for var in mon_variables.keys():
    locale = ''.join(path_source.split('/')[-2:-1])
    fname = var + '-' + locale + ".nc"
    print(fname)
    outf = os.path.join(path_destination, fname)
    ncf = nc.Dataset(outf, 'w', format="NETCDF4", clobber=True)
    now = datetime.datetime.now().strftime("%Y-%m-%dT%H:%M:%S")

    # Write global attributes
    ncf.Conventions = Conventions
    ncf.title = title
    ncf.institution = institution
    ncf.contact = contact
    ncf.comment = comment
    ncf.project = project
    ncf.references = references
    ncf.history = str(now) + ": Original ascii grids converted to geoTIFF, temporally merged, and metedata appended to netCDF container."

    index = 0

    for tif in tifs:

        # Where files start with PPT, Tave, Tmax, Tmin
        if tif.split('/')[-1].startswith(var):
            print(tif)
            ds = gdal.Open(tif)
            if ds is None:
                print("GTIff", tif, "failed to open!")

            a = ds.ReadAsArray()
            ama = np.ma.masked_array(a, np.equal(a, -9999))

            nlat, nlon = np.shape(ama)

            b = ds.GetGeoTransform()  # bbox, interval

            print(b)

            lon = np.arange(nlon) * b[1] + b[0]
            lat = np.arange(nlat) * b[5] + b[3]


            print("Lats:", nlat, lat)
            print("Lons:", nlon, lon)

            try:
                # First run of a variable to set up the NetCDF variables

                ncf.createDimension("time", None)
                ncf.createDimension("lat", nlat)
                ncf.createDimension("lon", nlon)

                lono = ncf.createVariable("lon", "i4", ("lon",), zlib=True)
                lono.axis = 'X'
                lono.units = "degrees_east"
                lono.long_name = "longitude"
                lono.standard_name = "longitude"
                lono[:] = lon
                print(lono[:])

                lato = ncf.createVariable("lat", "i4", ("lat",), zlib=True)
                lato.axis = 'Y'
                lato.units = "degrees_north"
                lato.long_name = "latitude"
                lato.standard_name = "latitude"
                lato[:] = lat
                print(lato[:])

                time = ncf.createVariable("time", "i4", ("time",), zlib=True)
                time.axis = 'T'
                time.units = "Month"
                time.long_name = "time"
                time.standard_name = "time"
                time[:] = 0

                variable = ncf.createVariable(mon_variables[var], "f4", ("lon", "lat", "time"), zlib=True, fill_value=-9999)
                variable.units = var_metadata[mon_variables[var]]["units"]
                variable.long_name = var_metadata[mon_variables[var]]["long_name"]
                variable.standard_name = var_metadata[mon_variables[var]]["standard_name"]
                variable.coordinates = var_metadata[mon_variables[var]]["coordinates"]
                variable.cell_methods = var_metadata[mon_variables[var]]["cell_methods"]

                variable[:] = np.empty((15399, 3236, 3), dtype=np.float)

            except RuntimeError:
                # Increment the index for each new file read

                index = index+1
                print(index)

            finally:
                # Write the array to the NetCDF as [lat, lon, time]

                variable[:, :, index] = a
                print(variable)

    ncf.close()

edit: I've looked at similar questions such as this one and have tried my best to replicate their methods but I still seem to run into problems. It could be that I'm forgetting a step.

  • If your main objective is to convert ESRI/ArcGIS grid (ASCII) files to netCDF, no programming is required. Just use gdal_translate. – Ralph Tee Apr 25 '18 at 17:38
  • I hadn't considered that, but that would be good fix for a situation less complicated than mine. I'm trying to restructure a database worth of ascii images in format and in file structure, and make them compliant with Climate and Forecast metadata standards (ie: cfconventions.org/Data/cf-conventions/cf-conventions-1.7/…). I have made some progress here so I'll update soon with what I've discovered. – Trevor J. Smith Apr 25 '18 at 18:09

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.