[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

Follow-up: My code for converting ASCII Grids to netCDF



Title: Follow-up: My code for converting ASCII Grids to netCDF

Hi folks,

I got a couple of really great responses to my earlier email asking for advice on converting ASCII raster grids to netCDF, and I wanted to share my solution with the group in case others find themselves needing similar functionality.

First, I tried to use the gdal_translate utility to create the netCDF file. This worked, however gdal_translate only has the very primitive ASCII raster metadata to work with and since it doesn't have a lot of configuration options to provide additional metadata the resulting netCDF is likewise void of useful metadata. In addition to that problem, gdal_translate only gives the beginning and end latitude and longitude points rather than an array of values for both. Some of my other code relies on having this lat/lon arrays, so I really needed another solution.

The code below is my current solution. Converting a 170 meg raster took 39 seconds, which I believe is comparable in speed to gdal_translate. One thing you may notice is I only used static methods in this code and so there's a little bit more parameter passing than some may be used to. Also, I used hashmaps to pass collections of parameters around, as noted in the code, these can be null if you like. The methods convertAsciiToNetcdfForVDM and convertAsciiDEMToNetcdfForVDM can be referred to for a specific example of how to use this class.

Comments and feedback appreciated, but since this now works for me, they aren't likely to be acted on ;-)

-Andrew

<code file="NetcdfHelper.java">
/**
 * Copyright: Copyright (c) 2006 Company: MBARI
 *
 * @author Andrew C. Chase
 */

package org.mbari.vdm;

import java.io.BufferedReader;
import java.io.File;
import java.io.FileReader;
import java.io.IOException;
import java.util.HashMap;
import java.util.Iterator;
import java.util.TimeZone;
import java.util.regex.Matcher;
import java.util.regex.Pattern;

import ucar.ma2.Array;
import ucar.ma2.DataType;
import ucar.ma2.Index;
import ucar.ma2.InvalidRangeException;
import ucar.nc2.Dimension;
import ucar.nc2.Group;
import ucar.nc2.NetcdfFileWriteable;

/* This library can be found at http://sourceforge.net/projects/timeandmoney/ */
import com.domainlanguage.timeutil.Clock;
import com.domainlanguage.timeutil.SystemClock;

public class NetcdfHelper {

    private static String TIME_ZONE_ID = "America/Los_Angeles";

    /**
     * @param gridName
     *            The grid file to convert
     * @param newNcFile
     *            The location to save the nc file to
     * @param gridDataTitle
     *            The name of the data that is in the ASCII grid. Null value
     *            okay
     * @param globalAttributes
     *            A collection of key value pairs (of type String:String) that
     *            will be added to the netcdf file as global attributes. Null
     *            value okay
     * @throws IllegalArgumentException
     * @throws IOException
     */
    public static void gridToNetcdf(File gridFile, File newNcFile, String gridDataTitle,
        HashMap globalAttributes, HashMap gridAttributes) throws IllegalArgumentException,
        IOException {
        if (gridFile == null || newNcFile == null || gridDataTitle == null) {
            throw new NullPointerException(
                "Both files and the gridDataTitle must all be non-null values");
        }
        NetcdfFileWriteable writeableNc = NetcdfFileWriteable.createNew(
            newNcFile.getAbsolutePath(), false);
        NetcdfHelper.addGlobalAttributes(writeableNc, globalAttributes);

        BufferedReader reader = new BufferedReader(new FileReader(gridFile));
        AsciiGridHeader header = AsciiGridHeader.parseHeader(reader);
        double[] latArray = createLatArray(header);
        double[] lonArray = createLonArray(header);

        writeableNc.addDimension("latitude", header.nrows);
        writeableNc.addDimension("longitude", header.ncols);
        Group rootGroup = writeableNc.getRootGroup();
        writeableNc.addVariable("latitude", DataType.DOUBLE, new Dimension[]{rootGroup
            .findDimension("latitude")});
        writeableNc.addVariableAttribute("latitude", "units", "degrees_north (+N/-S)");
        writeableNc.addVariable("longitude", DataType.DOUBLE, new Dimension[]{rootGroup
            .findDimension("longitude")});
        writeableNc.addVariableAttribute("longitude", "units", "degrees_east (+E/-W)");
        writeableNc.addVariable(gridDataTitle, DataType.FLOAT, new Dimension[]{
            rootGroup.findDimension("latitude"), rootGroup.findDimension("longitude")});
        if (gridAttributes != null) {
            addVariableAttributes(writeableNc, gridDataTitle, gridAttributes);
        }
        writeableNc.create();

        Array lonNcArray = createNcArray(lonArray);
        Array latNcArray = createNcArray(latArray);
        Array gridArray = getAsciiGridValues(reader, header);
        try {
            writeableNc.write("latitude", latNcArray);
            writeableNc.write("longitude", lonNcArray);
            writeableNc.write(gridDataTitle, gridArray);
        } catch (InvalidRangeException rangeE) {
            rangeE.printStackTrace();
            throw new RuntimeException(rangeE);
        }
        writeableNc.close();
    }

    private static void addGlobalAttributes(NetcdfFileWriteable writeableNc,
        HashMap globalAttributes) {
        Clock.setDefaultTimeZone(TimeZone.getTimeZone(TIME_ZONE_ID));
        Clock.setTimeSource(SystemClock.timeSource());
        if (globalAttributes == null || !globalAttributes.containsKey("creationDate")) {
            writeableNc.addGlobalAttribute("creationDate", Clock.today().toString("MM/dd/yyyy"));
        }
        if (globalAttributes == null || !globalAttributes.containsKey("lastModified")) {
            writeableNc.addGlobalAttribute("lastModified", Clock.today().toString("MM/dd/yyyy"));
        }
        if (globalAttributes != null) {
            for (Iterator iter = globalAttributes.keySet().iterator(); iter.hasNext();) {
                String key = (String) iter.next();
                String value = (String) globalAttributes.get(key);
                writeableNc.addGlobalAttribute(key, value);
            }
        }
    }

    private static void addVariableAttributes(NetcdfFileWriteable writeableNc, String variableName,
        HashMap variableAttributes) {
        for (Iterator iter = variableAttributes.keySet().iterator(); iter.hasNext();) {
            String attributeName = (String) iter.next();
            String attributeValue = (String) variableAttributes.get(attributeName);
            writeableNc.addVariableAttribute(variableName, attributeName, attributeValue);
        }
    }

    private static Array getAsciiGridValues(BufferedReader reader, AsciiGridHeader header)
        throws IOException {
        Array array = Array.factory(DataType.FLOAT, new int[]{header.nrows, header.ncols});
        Index arrayIndex = array.getIndex();
        String[] values;
        for (int i = 0; i < header.nrows; i++) {
            values = reader.readLine().split("\\s");
            for (int j = 0; j < header.ncols; j++) {
                arrayIndex.set(i, j);
                array.setDouble(arrayIndex, Double.parseDouble(values[j]));
            }
        }
        return array;
    }

    private static Array createNcArray(double[] array) {
        Array retval = Array.factory(DataType.DOUBLE, new int[]{array.length});
        Index index = retval.getIndex();
        for (int i = 0; i < array.length; i++) {
            index.set(i);
            retval.setDouble(index, array[i]);
        }
        return retval;
    }

    private static double[] createGridArray(int size, double start, double cellsize) {
        double[] retval = new double[size];
        for (int i = 0; i < size; i++) {
            retval[i] = start + (cellsize * i);
        }
        return retval;
    }

    private static double[] createLonArray(AsciiGridHeader header) {
        return createGridArray(header.ncols, header.xllcorner, header.cellsize);
    }

    private static double[] createLatArray(AsciiGridHeader header) {
        return createGridArray(header.nrows, header.yllcorner + (header.cellsize * header.nrows),
            -header.cellsize);
    }

    /**
     * This method is a convenience method for use within the VDM project, makes
     * assumptions about locations of files.
     *
     * @param gridName
     * @throws IOException
     * @throws IllegalArgumentException
     */
    public static void convertAsciiToNetcdfForVDM(String gridName, HashMap gAtts, String vName,
        HashMap vAtts) throws IllegalArgumentException, IOException {
        if (!gridName.endsWith("asc")) {
            throw new IllegalArgumentException("Must be an ascii raster grid with extension .asc");
        }
        File gridFile = new File(NetcdfHelper.class.getResource("/grids").getFile(), gridName);
        String newNcFileName = gridName.substring(0, gridName.lastIndexOf('.'));
        File newNcFile = new File(NetcdfHelper.class.getResource("/grids").getFile(), newNcFileName
            + ".nc");
        gridToNetcdf(gridFile, newNcFile, vName, gAtts, vAtts);
    }

    public static void convertAsciiDEMToNetcdfForVDM(String gridName)
        throws IllegalArgumentException, IOException {
        HashMap vHash = new HashMap();
        vHash.put("units", "meters");
        vHash.put("positive", "up");
        vHash.put("missing_value", "-99999");
        vHash.put("long_name", "Elevation");
        HashMap gHash = new HashMap();
        gHash.put("keywords", "dem bathymetry grid");
        gHash.put("title", "DEM Grid from " + gridName);
        gHash.put("creator_name", "Andrew Chase");
        gHash.put("institution", "MBARI");
        convertAsciiToNetcdfForVDM(gridName, gHash, "elevation", vHash);
    }

    public static void main(String[] args) {
        long start = System.currentTimeMillis();
        try {
            convertAsciiDEMToNetcdfForVDM("hrbath.asc");
        } catch (IOException e) {
            e.printStackTrace();
        }
        System.out.println("total time = " + ((System.currentTimeMillis() - start)/ 1000) + " seconds");
    }
}

class AsciiGridHeader {

    final int ncols;
    final int nrows;
    final double xllcorner;
    final double yllcorner;
    final double cellsize;
    final double NODATA_value;

    /*
     * This method should only be called by the parseHeader method.
     */
    private AsciiGridHeader(int ncols, int nrows, double xllcorner, double yllcorner,
        double cellsize, double nodata) {
        this.ncols = ncols;
        this.nrows = nrows;
        this.xllcorner = xllcorner;
        this.yllcorner = yllcorner;
        this.cellsize = cellsize;
        this.NODATA_value = nodata;
    }

    public static AsciiGridHeader parseHeader(BufferedReader gridReader) throws IOException {
        Pattern p = Pattern.compile("(\\w+)\\s+(-?\\d+.*)");
        Matcher m = p.matcher(gridReader.readLine());
        if (!m.matches()) {
            throw new RuntimeException("Header does not obey format contract");
        }
        int ncols = Integer.parseInt(m.group(2));
        m = p.matcher(gridReader.readLine());
        if (!m.matches()) {
            throw new RuntimeException("Header does not obey format contract");
        }
        int nrows = Integer.parseInt(m.group(2));
        m = p.matcher(gridReader.readLine());
        if (!m.matches()) {
            throw new RuntimeException("Header does not obey format contract");
        }
        double xllcorner = Double.parseDouble(m.group(2));
        m = p.matcher(gridReader.readLine());
        if (!m.matches()) {
            throw new RuntimeException("Header does not obey format contract");
        }
        double yllcorner = Double.parseDouble(m.group(2));
        m = p.matcher(gridReader.readLine());
        if (!m.matches()) {
            throw new RuntimeException("Header does not obey format contract");
        }
        double cellsize = Double.parseDouble(m.group(2));
        m = p.matcher(gridReader.readLine());
        if (!m.matches()) {
            throw new RuntimeException("Header does not obey format contract");
        }
        double nodata = Double.parseDouble(m.group(2));
        return new AsciiGridHeader(ncols, nrows, xllcorner, yllcorner, cellsize, nodata);
    }
}
</code>

====My Original Email====

> Hi Folks,
>
> I'm trying to do something relatively simple, but I haven't been able
> to find any examples of others doing it before me. Just thought I'd
> ask here before going off and writing the code myself.
>
> What I want to do is take a DEM stored in the  ESRI ASCII Raster
> format and convert it to a netcdf file. My ascii raster has the
> standard header format like so:
> ncols         5001
> nrows         4001
> xllcorner     -125
> yllcorner     34
> cellsize      0.001
> NODATA_value  -9999
>
> Anybody got some code to do this? :-)
>
> If not, no biggee, just thought I'd ask.


Andrew C. Chase
Software Engineer
Monterey Bay Aquarium Research Institute