Opened 7 years ago

Closed 7 years ago

#4432 closed defect (fixed)

gdal_rasterize produces error when burning vector dataset into netCDF dataset directly with -a_nodata nan

Reported by: hsumanto Owned by: etourigny
Priority: normal Milestone:
Component: default Version: unspecified
Severity: normal Keywords: netcdf, create, gdal_rasterize
Cc: etourigny

Description

I am trying to perform the below with the latest gdal in trunk:

gdal_rasterize -of netCDF -tr 500.0 500.0 -te 519223.62281715404 -4538725.344003787 1452723.622817154 -3983225.3440037873 -a_srs EPSG:3112 -a_nodata nan -ot float32 -a value -l SHAPE_TEST SHAPE_TEST.shp Output_Invalid.nc

and I am getting these errors:

ERROR 1: netCDF scanline fetch failed: NetCDF: Operation not allowed in define mode
ERROR 1: Output.nc, band 1: IReadBlock failed at X offset 0, Y offset 0
ERROR 1: GetBlockRef? failed at X block offset 0, Y block offset 0

It seems that it doesn't like -a_nodata nan because I could burn vector dataset into netCDF without error if I remove -a_nodata nan. Still wondering the reason why -a_nodata nan is causing the error.

Burn vector dataset into netCDF (OK without -a_nodata nan)

gdal_rasterize -of netCDF -tr 500.0 500.0 -te 519223.62281715404 -4538725.344003787 1452723.622817154 -3983225.3440037873 -a_srs EPSG:3112 -ot float32 -a value -l SHAPE_TEST SHAPE_TEST.shp Output_OK.nc

Burn vector dataset into Geotiff Raster (OK even with -a_nodata nan)

gdal_rasterize -tr 500.0 500.0 -te 519223.62281715404 -4538725.344003787 1452723.622817154 -3983225.3440037873 -a_srs EPSG:3112 -a_nodata nan -ot float32 -a value -l SHAPE_TEST SHAPE_TEST.shp Output_Invalid.tif

Attachments (4)

SHAPE_TEST.shp (575.1 KB) - added by hsumanto 7 years ago.
vector dataset file
SHAPE_TEST.dbf (513.5 KB) - added by hsumanto 7 years ago.
vector dataset file
SHAPE_TEST.prj (470 bytes) - added by hsumanto 7 years ago.
vector dataset file
SHAPE_TEST.shx (164.4 KB) - added by hsumanto 7 years ago.
vector dataset file

Download all attachments as: .zip

Change History (11)

Changed 7 years ago by hsumanto

Attachment: SHAPE_TEST.shp added

vector dataset file

Changed 7 years ago by hsumanto

Attachment: SHAPE_TEST.dbf added

vector dataset file

Changed 7 years ago by hsumanto

Attachment: SHAPE_TEST.prj added

vector dataset file

Changed 7 years ago by hsumanto

Attachment: SHAPE_TEST.shx added

vector dataset file

comment:1 Changed 7 years ago by Even Rouault

Cc: etourigny added

comment:2 Changed 7 years ago by etourigny

Keywords: create added
Owner: changed from warmerdam to etourigny
Status: newassigned

It seems that gdal_rasterize sets nodata values after setting projection info, which I had not tested before in the new Create() method (which gdal_rasterize uses).

IReadBlock() is missing a call to get out of define mode, so the fix is to add a call to SetDefinemode?(False), just like in IWriteBlock(). I don't think this should be detrimental to performance, as the nodata value was set to default previously and does not need to be allocated, only modified.

Try this simple fix in IReadBlock():

            Taken += start[panBandZPos[i]] * Sum;
        }
    }

+    /* make sure we are in data mode */
+    ( ( netCDFDataset * ) poDS )->SetDefineMode( FALSE );

+    /* read data according to type */
    if( eDataType == GDT_Byte ) 
    {

Can you try with other nodata values also?

comment:3 Changed 7 years ago by etourigny

Oh and in the furure - please send shapefiles in a single zip files, it makes it easier to download and also dramatically smaller - thanks.

comment:4 Changed 7 years ago by hsumanto

Thanks Etienne for the quick response. I will send them as a single zip next time.

I try the fix above, it fixes the issue:

WITH nan as nodata value
gdal_rasterize -of netCDF -tr 500.0 500.0 -te 519223.62281715404 -4538725.344003787 1452723.622817154 -3983225.3440037873 -a_srs EPSG:3112 -a_nodata nan -ot float32 -a value -l SHAPE_TEST SHAPE_TEST.shp Output.nc
0...10...20...30...40...50...60...70...80...90...100 - done.

WITH -9999 as nodata value
gdal_rasterize -of netCDF -tr 500.0 500.0 -te 519223.62281715404 -4538725.344003787 1452723.622817154 -3983225.3440037873 -a_srs EPSG:3112 -a_nodata -9999 -ot float32 -a value -l SHAPE_TEST SHAPE_TEST.shp Output2.nc
0...10...20...30...40...50...60...70...80...90...100 - done.

Check gdalinfo on both files, they seem OK.

Then I just try to check the content by converting the netcdf into ASCII

gdal_translate -of AAIGrid Output2.nc Output2.asc
Input file size is 1867, 1111
0ERROR 1: netcdf error #-38 : NetCDF: Operation not allowed in data mode .
at (netcdfdataset.cpp,SetDefineMode,1339)

...10...20...30...40...50...60...70...80...90...100 - done.

gdal_translate -of AAIGrid Output2.nc Output2.asc
Input file size is 1867, 1111
0ERROR 1: netcdf error #-38 : NetCDF: Operation not allowed in data mode .
at (netcdfdataset.cpp,SetDefineMode,1339)

...10...20...30...40...50...60...70...80...90...100 - done.

The ASCII datasets are still being produced but I noticed the similar warning about SetDefineMode? do appear again here.

comment:5 Changed 7 years ago by etourigny

You should also get the same error with gdalinfo -stats , as it also calls IReadBlock(). The call to nc_enddef/nc_redef (in SetDefineMode?) should not happen when only reading a file, not writing.

Here is the proper fix in SetDefineMode?() (in addition to the previous one):

int netCDFDataset::SetDefineMode( int bNewDefineMode )
{
    /* if dataset is in read-only mode, do nothing */
    if( GetAccess() == GA_ReadOnly ) return CE_None;

    /* if already in new define mode, do nothing */
    if ( bDefineMode == bNewDefineMode ) return CE_None;
...
}

I will commit this in the next few days, please confirm it works well, thanks.

comment:6 Changed 7 years ago by hsumanto

Yeah, that fix resolves the issues which I have encountered so far. Thanks, Etienne.

comment:7 Changed 7 years ago by etourigny

Keywords: gdal_rasterize added
Resolution: fixed
Status: assignedclosed

fixed in 1.9 (r23760) and trunk (r23759)

Note: See TracTickets for help on using tickets.