Opened 16 years ago

Closed 16 years ago

#2128 closed defect (invalid)

Bug in reading georeferencing of Arc Info ASCII Grid

Reported by: moellney Owned by: warmerdam
Priority: normal Milestone:
Component: default Version: 1.5.0
Severity: normal Keywords:
Cc:

Description

The ARC/Info file:

ncols 1
nrows 1
xllcorner 0.0000000000
yllcorner 0.0000000000
cellsize 0.0083333333
ycellsize 0.0083333333
nodata_value -32767
-32767 

Results in an adfGeoTransform:

0.0
0.0083333333
0.0
0.0
0.0
-0.0083333333

but shouldn't it be

-0.0083333333 / 2.0
0.0083333333
0.0
0.0083333333 / 2.0
0.0
-0.0083333333

to reflect the definition of the adfGeoTransform reference point to be the upper left corner of the area of the upper left data point?

So I think in aaigriddataset.cpp (http://trac.osgeo.org/gdal/browser/branches/1.5/gdal/frmts/aaigrid/aaigriddataset.cpp#L444) the calculation is wrong.

444   if ((i = CSLFindString( papszTokens, "xllcorner" )) >= 0 &&
445 	        (j = CSLFindString( papszTokens, "yllcorner" )) >= 0 )
446 	    {
447 	        poDS->adfGeoTransform[0] = atof( papszTokens[i + 1] );
448 	        poDS->adfGeoTransform[1] = dfCellDX;
449 	        poDS->adfGeoTransform[2] = 0.0;
450 	        poDS->adfGeoTransform[3] = atof( papszTokens[j + 1] )
451 	            + poDS->nRasterYSize * dfCellDY;
452 	        poDS->adfGeoTransform[4] = 0.0;
453 	        poDS->adfGeoTransform[5] = - dfCellDY;
454 	    }
455 	    else if ((i = CSLFindString( papszTokens, "xllcenter" )) >= 0 &&
456 	             (j = CSLFindString( papszTokens, "yllcenter" )) >= 0 )
457 	    {
458 	        poDS->SetMetadataItem( GDALMD_AREA_OR_POINT, GDALMD_AOP_POINT );
459 	
460 	        poDS->adfGeoTransform[0] = atof(papszTokens[i + 1]) - 0.5 * dfCellDX;
461 	        poDS->adfGeoTransform[1] = dfCellDX;
462 	        poDS->adfGeoTransform[2] = 0.0;
463 	        poDS->adfGeoTransform[3] = atof( papszTokens[j + 1] )
464 	            - 0.5 * dfCellDY
465 	            + poDS->nRasterYSize * dfCellDY;
466 	        poDS->adfGeoTransform[4] = 0.0;
467 	        poDS->adfGeoTransform[5] = - dfCellDY;
468 	    }

I think both calculations are wrong. The (X|Y)LLCORNER defines the position of the dataset in the lower right corner, so we have to add the "half-pixel"-size offset here.

The (X|Y)LLCENTER seems to define the position of the center of the lower left tile of the dataset (defined by the 4 points in the lower left corner of the dataset). So we have to add "full-pixel"-size.

Michael

Attachments (1)

nine_values.img (7.6 KB ) - added by moellney 16 years ago.
Nine DTM values around 0°, 0° generated by geoengine.nima.mil

Download all attachments as: .zip

Change History (4)

comment:1 by warmerdam, 16 years ago

Status: newassigned

Michael,

I have been operating on the assumption that the (X|Y)LLCORNER defines the lower left corner of the lower left pixel. As such, we would want adfGeoTransform[0] to be XLLCORNER. I have been operating on the assumption that XLLCENTER is the center of the bottom left pixel, hence a half pixel offset.

If you have reason to believe these assumptions are wrong, then can you provide some explanation? Is this based on some document? Careful export of data to this format using ArcGIS?

in reply to:  1 comment:2 by moellney, 16 years ago

Replying to warmerdam:

Michael,

I have been operating on the assumption that the (X|Y)LLCORNER defines the lower left corner of the lower left pixel. As such, we would want adfGeoTransform[0] to be XLLCORNER. I have been operating on the assumption that XLLCENTER is the center of the bottom left pixel, hence a half pixel offset.

I did not find 'hard' evidence in the internet. And I googled for xllcorner and xllcenter in the same document and found e.g.:
http://geotools.codehaus.org/ArcInfo+ASCII+Grid+format (last section on page seems to very clear on this topic)
http://www.anuva.de/service_arcforum.php?action=vthread&forum=2&topic=2866
http://www.csc.noaa.gov/crs/tcm/instructions.html#binaryrasters

So I now think you are right.

If you have reason to believe these assumptions are wrong, then can you provide some explanation? Is this based on some document? Careful export of data to this format using ArcGIS?

The example I gave in my issue was generated by geoengine.nima.mil. Out of the respect from the data they provide I assumed that they should be coding them right:

I requested just one DTM value for the sake of simpicity of the file: elevation at 0° latitude and 0° longitude.

Requesting 3 values in latitude and 3 in longitude gives a sheet of 9 values:

ncols 3
nrows 3
xllcorner -0.0083333333
yllcorner -0.0083333333
cellsize 0.0083333333
ycellsize 0.0083333333
nodata_value -32767
-32767 -32767 -32767 
-32767 -32767 -32767 
-32767 -32767 -32767 

Again you can see, that the (x|y)llcorner of the dataset is just one cellsize away from the center of data (0°, 0°) that I requested, and not 1.5 cellsizes away.

You can read out the DTM data in different file formats on this site, so I requested the same data in Erdas IMG format, too and attached it to this issue. GDAL delivers different georeferencings for the same data depending on the format used. The referencing for IMG is right... for ESRI ASCII GRID the referencing is wrong.

However the fault might be on the geoengine.nima.mil, setting (x|y)llcorner instead of (x|y)llcenter in their data....

I will contact them for more details.

by moellney, 16 years ago

Attachment: nine_values.img added

Nine DTM values around 0°, 0° generated by geoengine.nima.mil

comment:3 by moellney, 16 years ago

Resolution: invalid
Status: assignedclosed

I'm in contact with geoengine.nima.mil and one support worker confirmed me, that there seems to be a bug in the ARC ASCII grid export.

So I close this request as resolved invalid.

Note: See TracTickets for help on using tickets.