Opened 19 years ago

Last modified 19 years ago

#751 closed defect (fixed)

Typo in hfadataset.cpp

Reported by: yangz@… Owned by: warmerdam
Priority: high Milestone:
Component: GDAL_Raster Version: unspecified
Severity: normal Keywords:
Cc:

Description

there are two typos in hfadataset.cpp. 
The word "Azimuthal" was typed as "Azitmuthal" on line 1241 and 1250. As a
result, if you set a projection to to "Lambert Azimuthal Equal-area" or
"Azimuthal Equidistant", the resulting image will not be valid in projection.

Change History (6)

comment:1 by yangz@…, 19 years ago

After correcting the typo in the source code, there seems to be another related
problem with this: the resulting projection of an image is not the same as you
asked for. The following is an example:

if create a new image and set the projection to be 

PROJCS["Lambert Azimuthal
Equal-area",GEOGCS["Undefined",DATUM["Undefined",SPHEROID["Sphere of Radius
6370997m",6370997,"1.#INF"]],PRIMEM["Greenwich",0],
UNIT["degree",0.0174532925199433]],
PROJECTION["Lambert_Azimuthal_Equal_Area"],
PARAMETER["latitude_of_center",45],
PARAMETER["longitude_of_center",-100],
PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["meters",1]]

but the projection for the created new image comes out as 

PROJCS["Lambert Azimuthal
Equal-area",GEOGCS["Undefined",DATUM["Undefined",SPHEROID["Sphere of Radius
6370997m",6370997,447.2538558200429],TOWGS84[0,0,0,0,0,0,0]],
PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433]],
PROJECTION["Lambert_Azimuthal_Equal_Area"],
PARAMETER["latitude_of_center",45],PARAMETER["longitude_of_center",-100],
PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["meters",1]]

As a result, the image can not be displayed with the original image properly
within Imagine.

comment:2 by warmerdam, 19 years ago

I have patched the spelling of azimuthal as suggested.


I'm not sure why you try to set a spheroid with an inverse flattening
of "1.#INF".  I believe the inverse flattening should be infinite
for spheres, but that we use the special value "0" as a marker for this case.

How did you come up with the input WKT definition?

comment:3 by yangz@…, 19 years ago

The input WKT definition is the result of getProjectionRef() from an existing
Imagine file, and it seems like that if you define the projection for an imagine
file with "Lambert Azimuthal Equal-area" using the GUI interface, the resulting
WKT string will be the input WKT string I have been trying to use.

If you need, I can send you a test image.

comment:4 by warmerdam, 19 years ago

Yang, 

Yes, I would appreciate it if you could make a smallish example available.
I think there is a bug in the way I originally form the WKT for images
on a sphere. 

comment:5 by yangz@…, 19 years ago

Here is an example to show the projection problem. I will send the images to you
in another email.

#include "gdal_priv.h"
#include "cpl_string.h"
#include "ogr_spatialref.h"

void CreateNewImg() {
	GDALAllRegister();
	GDALDataset *poDstDS, *poSrcDS, *poReopenDS;
	
	/* Project for the src.img is
	
	   PROJCS["Lambert Azimuthal Equal-area",GEOGCS["Undefined",
	   DATUM["Undefined",SPHEROID["Sphere of Radius 6370997m",6370997,"1.#INF"]],
	   PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433]],
	   PROJECTION["Lambert_Azimuthal_Equal_Area"],
	   PARAMETER["latitude_of_center",45],PARAMETER["longitude_of_center",-100],
	   PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["meters",1]].
	
	*/
	const char* pszSrcFileName = "H:\\testimg\\src.img";
	const char* pszNewFileName = "H:\\testimg\\new.img";
	
	// open the source image to get the projection information
	poSrcDS = (GDALDataset *)GDALOpen(pszSrcFileName, GA_ReadOnly);
	const char * pszPrj = poSrcDS->GetProjectionRef();
	printf("Source Projection is:\n\t %s.\n", pszPrj);

	// create a new image with the same projection information as 
	// the source image
	double adfGeoTransform[6] = {-945000.0, 1000, 0, -161000.0, 0, -1000};
	OGRSpatialReference oSRS;
	GDALRasterBand *poBand;
	GByte abyRaster[512*512];
	poDstDS =
GetGDALDriverManager()->GetDriverByName("HFA")->Create(pszNewFileName, 512, 512,
1, GDT_Byte, NULL);
	poDstDS->SetGeoTransform( adfGeoTransform );
	poDstDS->SetProjection( pszPrj );

	poBand = poDstDS->GetRasterBand(1);
	poBand->RasterIO( GF_Write, 0, 0, 512, 512,	abyRaster, 512, 512, GDT_Byte, 0, 0
);    

	//No surprise. Now the projection information is still correct.
	printf("New Image Projection is:\n\t %s.\n", poDstDS->GetProjectionRef());

	//Reopen the same file in shared mode, the projection is still correct. 
	poReopenDS = (GDALDataset *)GDALOpenShared(pszNewFileName, GA_ReadOnly);
	printf("Open New Image (shared mode) Projection is:\n\t %s.\n",
poReopenDS->GetProjectionRef());

	//close the dataset
	delete poDstDS;
	delete poSrcDS;

	// reopen the newly created image, the projection is wrong now!
	// Something is wrong in when writes the WKT projection to the image.
	poReopenDS = (GDALDataset *)GDALOpen(pszNewFileName, GA_ReadOnly);
	printf("Reopen Image Projection is:\n\t %s.\n", poReopenDS->GetProjectionRef());

	//close dataset
	delete poReopenDS;
}

int _tmain(int argc, _TCHAR* argv[])
{
	CreateNewImg();
	return 0;
}

comment:6 by warmerdam, 19 years ago

I have confirmed a bug reading HFA files with spherical ellipsoids. I have
corrected this in hfadataset.cpp now and I am now able to copy the source file
to a new .img file with no ill effects.

Please re-open if you find a problem persists.

Note: See TracTickets for help on using tickets.