ASC to JSON conversion. Where are the attributes?

Question

I am converting MetOffice Nimrod .dat files to GeoJSON via the asc format.

The MetOffice dat format data is a proprietary type and I am using a script found here that converts from dat to asc.

Using QGIS I can see a reasonable looking output from this asc file.

From the dat file I use gdal to convert to json. I also do some geospatial clipping and convert to the epsg:4326 projection. I am confident from testing other files that these sections are correct.

I iterate over the asc files and use gdal, as is show below.

for file in glob.glob("*.asc"):
        print("Converting file " + file + " to json")
        src_ds = gdal.Open(file)
        if src_ds.RasterCount !=1:
            print "Raster issue - more than one band detected"
            sys.exit(1)
        if src_ds is None:
            print 'Unable to open %s' % src_filename
            sys.exit(1)
        try:
            srcband = src_ds.GetRasterBand(1)
        except RuntimeError, e:
            print 'Band ( %i ) not found' % band_num
            print e
            sys.exit(1)
        dst_layername = file[:-4]
        drv = ogr.GetDriverByName("GeoJSON")
        dst_ds = drv.CreateDataSource( dst_layername + ".json" )
        dst_layer = dst_ds.CreateLayer(dst_layername, srs = None )
        gdal.Polygonize( srcband, None, dst_layer, -1, [], callback=None )

This outputs a reasonable looking GeoJSON file, as shown below. I am happy that this is geospatially correct.

However, my issue lies in the attributes for each of these polygons. The MetOffice NIMROD data contains precipitation measurements.

If I open the json file in a text editor, I see the following:

{
"type": "FeatureCollection",
"crs": { "type": "name", "properties": { "name": "urn:ogc:def:crs:OGC:1.3:CRS84" } },

"features": [
{ "type": "Feature", "properties": { }, "geometry": { "type": "Polygon", "coordinates": [ [ [ -0.1872638449, 51.653290137600003 ], [ -0.1728164805, 51.653066222 ], [ -0.1735389107, 51.635092885699997 ], [ -0.1879805698, 51.635316657899999 ], [ -0.1872638449, 51.653290137600003 ] ] ] } },
{ "type": "Feature", "properties": { }, "geometry": { "type": "Polygon", "coordinates": [ [ [ -0.1735389107, 51.635092885699997 ], [ -0.1590974392, 51.634867338200003 ], [ -0.1594612233, 51.625880718799998 ], [ -0.1738998445, 51.626106194099997 ], [ -0.1735389107, 51.635092885699997 ] ] ] } },
{ "type": "Feature", "properties": { }, "geometry": { "type": "Polygon", "coordinates": [ [ [ -0.1594612233, 51.625880718799998 ], [ -0.1450227909, 51.625653468899998 ], [ -0.1453892349, 51.616666906500001 ], [ -0.1598248185, 51.616894083799998 ], [ -0.1594612233, 51.625880718799998 ] ] ] } },
{ "type": "Feature", "properties": { }, "geometry": { "type": "Polygon", "coordinates": [ [ [ 0.1687182185, 51.539341131199997 ], [ 0.1831263429, 51.53907385 ], [ 0.1826964281, 51.5300889699 ], [ 0.168291136, 51.530356165599997 ], [ 0.1687182185, 51.539341131199997 ] ] ] } },
{ "type": "Feature", "properties": { }, "geometry": { "type": "Polygon", "coordinates": [ [ [ -0.4960296, 51.2999988 ], [ -0.4960296, 51.6806202 ], [ 0.275826, 51.6806202 ], [ 0.275826, 51.2999988 ], [ -0.4960296, 51.2999988 ] ] ] } }
]
}

There is no property or attribute metric here which gives a number to the precipitation for this polygon.

Where am I going wrong? Is there a further gdal argument for bringing the file attributes across from the asc file too?

Answer 1

I think the issue is you need to explicitly add a field to your output layer to store the rainfall value for each raster pixel when it gets converted to a polygon. Using your script, I added an extra two lines just before the last line:

import glob
import sys
import gdal
import osgeo.ogr

for file in glob.glob("*.asc"):
        print("Converting file " + file + " to json")
        src_ds = gdal.Open(file)
        if src_ds.RasterCount != 1:
            print "Raster issue - more than one band detected"
            sys.exit(1)
        if src_ds is None:
            print 'Unable to open %s' % file
            sys.exit(1)
        try:
            srcband = src_ds.GetRasterBand(1)
        except RuntimeError, e:
            print 'Band ( %i ) not found' % srcband
            print e
            sys.exit(1)
        dst_layername = file[:-4]
        drv = osgeo.ogr.GetDriverByName("GeoJSON")
        dst_ds = drv.CreateDataSource(dst_layername + ".json")
        dst_layer = dst_ds.CreateLayer(dst_layername, srs=None)

        precipField = osgeo.ogr.FieldDefn('PRECIP', osgeo.ogr.OFTInteger)
        dst_layer.CreateField(precipField)

        gdal.Polygonize(srcband, None, dst_layer, 0, [], callback=None)

Then in my JSON file I get:

{
"type": "FeatureCollection",
"features": [
{ "type": "Feature", "properties": { "PRECIP": 12 }, "geometry": { "type": "Polygon", "coordinates": [ [ [ 126500.0, 699500.0 ], [ 126500.0, 698500.0 ], [ 131500.0, 698500.0 ], [ 131500.0, 699500.0 ], [ 126500.0, 699500.0 ] ] ] } },
{ "type": "Feature", "properties": { "PRECIP": 27 }, "geometry": { "type": "Polygon", "coordinates": [ [ [ 131500.0, 699500.0 ], [ 131500.0, 698500.0 ], [ 136500.0, 698500.0 ], [ 136500.0, 699500.0 ], [ 131500.0, 699500.0 ] ] ] } },
{ "type": "Feature", "properties": { "PRECIP": 8 }, "geometry": { "type": "Polygon", "coordinates": [ [ [ 136500.0, 699500.0 ], [ 136500.0, 698500.0 ], [ 141500.0, 698500.0 ], [ 141500.0, 699500.0 ], [ 136500.0, 699500.0 ] ] ] } },
{ "type": "Feature", "properties": { "PRECIP": 3 }, "geometry": { "type": "Polygon", "coordinates": [ [ [ 161500.0, 699500.0 ], [ 161500.0, 698500.0 ], [ 166500.0, 698500.0 ], [ 166500.0, 699500.0 ], [ 161500.0, 699500.0 ] ] ] } },

...

This answer here helped me: gdal polygonize in python error, blank polygon created

Answer 2

The problem is that GeoJSON is a vector format and asc is a raster format. You could try exporting to a commonly used (and smaller) raster format like GeoTif.