GDALDriver::CreateCopy() fails for "USGSDEM"

[harishkumargvs]

Trying to create a DEM from in-memory data but, GDALDriver::CreateCopy() fails by returning NULL. Here's the code:

USES_CONVERSION;

GDALDataset *poDST = NULL;
	OGRSpatialReference oSRS;
	char *prj4 = NULL;
	GDALAllRegister();

	// Get the GDAL driver
	GDALDriver *poDriver = NULL;
	poDriver = GetGDALDriverManager()->GetDriverByName("USGSDEM");
	if (!poDriver)	goto cleanup;

	if (!poDriver->GetMetadataItem(GDAL_DCAP_CREATECOPY) ||
		poDriver->GetMetadataItem(GDAL_DCAP_CREATECOPY)[0] != 'Y') goto cleanup;

	// GDAL georeferencing
	double trf[6] = { 423002.41999999975,
		4.9999999999999973,
		0.00000000000000000,
		1544418.0500000000,
		0.00000000000000000,
		-4.9999999999999973 };

	// Create an in memory dataset of the copy
	char inmem[256];
	long nCols(2194), nRows(2188);
	float *pfData = (float*)malloc(nCols * nRows * sizeof(float));
	if (!pfData) goto cleanup;

	for (long idx = 0; idx < nCols * nRows; idx++)
		pfData[idx] = (float) 600;

	sprintf_s(inmem, 256,
		"MEM:::DATAPOINTER=0x%p,"
		"PIXELS=%d,"
		"LINES=%d,"
		"BANDS=1,"
		"DATATYPE=Float32,"
		"PIXELOFFSET=4",
		pfData, nCols, nRows);

	GDALDataset *poSRS = NULL;
	poSRS = (GDALDataset*)GDALOpen(inmem, GA_ReadOnly);
	if (!poSRS) goto cleanup;

	if (poSRS->SetMetadataItem(GDALMD_AREA_OR_POINT, GDALMD_AOP_AREA)) goto cleanup;

	if (poSRS->SetGeoTransform(trf)) goto cleanup;

	GDALRasterBand* poBand = poSRS->GetRasterBand(1);
	if (!poBand) goto cleanup;

	if (poBand->SetNoDataValue(FLT_MAX)) goto cleanup;

	const wchar_t* path = L"Test.dem";

	char szCSInfo[] = "PROJCS[\"NAD_1983_StatePlane_Alabama_East_FIPS_0101_Feet\",GEOGCS[\"GCS_North_American_1983\",DATUM[\"D_North_American_1983\",SPHEROID[\"GRS_1980\",6378137,298.257222101]],PRIMEM[\"Greenwich\",0],UNIT[\"Degree\",0.017453292519943295]],PROJECTION[\"Transverse_Mercator\"],PARAMETER[\"latitude_of_origin\",30.5],PARAMETER[\"central_meridian\",-85.833333333333329],PARAMETER[\"scale_factor\",0.99995999999999996],PARAMETER[\"false_easting\",656166.66666666674],PARAMETER[\"false_northing\",0],UNIT[\"Foot_US\",0.30480060960121924]]";
	char *sz = szCSInfo;

	if (oSRS.importFromWkt(&sz)) goto cleanup;

	if (oSRS.morphFromESRI()) goto cleanup;

	if (oSRS.exportToWkt(&prj4) || !prj4) goto cleanup;

	if (poSRS->SetProjection(prj4)) goto cleanup;

	// Write the data to file
	poDST = poDriver->CreateCopy(W2A(path), poSRS, FALSE, NULL, NULL, NULL);
	if (!poDST){
		std::cout << "poDriver->CreateCopy Failed!" << std::endl;
		goto cleanup;
	}

cleanup:
	if (pfData)	{
		free(pfData);
		pfData = NULL;
	}
	if (prj4)
		OGRFree(prj4);

	if (poDST)
		GDALClose(poDST);

	if (poSRS)
		GDALClose(poSRS);

Note: GDALDriver::CreateCopy() used to work fine with gdal17.dll versioned 1.7.2.0 but with gdal110.dll versioned 1.10.1.0 it fails...

[Even Rouault]

I'm surprised this worked with previous GDAL releases. I believe that the WKT you provide is incorrect. It has an empty string for the PROJECTION[""] parameter. Hence the warping done internaly in the USGSDEM driver fails to work not managing to transform this WKT to a PROJ.4 string. Perhaps there was some change in the ESRI morphing code, that used to work with this strange WKT before...

Try directly with the OGC WKT at ​http://spatialreference.org/ref/esri/102654/ogcwkt/. I actually tried it but it gave some other error during reprojection which makes me think it doesn't expect that kind of source SRS.

[harishkumargvs]

I have modified the test-driver code slightly to replicate the test data I am using excepting the Z-Values, I have hard-coded it as 600. As you can see that having appropriate PROJECTION string in place, the latest GDAL Driver's CreateCopy() is now able to create the output but I get the following error along with a warning message: ERROR 1: Too many points (441 out of 441) failed to transform, unable to compute output bounds. Warning 1: Unable to compute source region for output window 0,0,2194,2188, skipping. As a result, the output DEM is not usable i.e., it simply crashes our elevation data viewer application. The same code, if I link it with the older version of gdal_i.lib to let our application use gdal17.dll at run-time, I get the same error message but no warning this time. The output DEM is not usable even in this case.

Now, I tried with an empty PROJECTION string and ran it with gdal17.dll, I get the following error message: ERROR 6: No translation for to PROJ.4 format is known. Interestingly, this time the output DEM is viewable with our elevation data viewer application.

I would like to know the reason as to why the latest GDAL fails to create correct output even after providing it with the correct WKTEXT?

[Kyle Shannon]

harishkumargvs,

Is GDAL_DATA properly set? Some of those errors are frequently reported if the support files cannot be found.

[harishkumargvs]

Kyle, I understand that you are asking whether we are setting GDAL_DATA System Environmental Variable or not. If so, we are not doing that. Now, as you said that that these errors might result due to not setting it, I have tried setting GDAL_DATA with the path that contains miscellaneous CSV files, namely - coordinate_axis, ellipsoid, gcs, gcs.override, gdal_datum, pcs, pcs.override, prime_meridian, stateplane, unit_of_measure along with three wkt extension files - cubewerx_extra, epsg, esri_extra. Even then my test-driver application throws the same errors. I believe that these are the support files that you were talking about. If you think that there are some other support files that are missing in my aforementioned list which might be causing these errors, please let me know. If not what else could be the reason for these errors. Did you tried executing the test-code on your machine? As I have already mentioned, without setting GDAL_DATA, and without providing PROJECTION string gdal17.dll was able to create the output that can be viewed in our elevation data viewer of course, I get an error - ERROR 6: No translation for to PROJ.4 format.

[harishkumargvs]

Replying to harishkumargvs: I would request someone to please comment on what is the reason for the errors in my previous comment. Thanks in advance...

Kyle, I understand that you are asking whether we are setting GDAL_DATA System Environmental Variable or not. If so, we are not doing that. Now, as you said that that these errors might result due to not setting it, I have tried setting GDAL_DATA with the path that contains miscellaneous CSV files, namely - coordinate_axis, ellipsoid, gcs, gcs.override, gdal_datum, pcs, pcs.override, prime_meridian, stateplane, unit_of_measure along with three wkt extension files - cubewerx_extra, epsg, esri_extra. Even then my test-driver application throws the same errors. I believe that these are the support files that you were talking about. If you think that there are some other support files that are missing in my aforementioned list which might be causing these errors, please let me know. If not what else could be the reason for these errors. Did you tried executing the test-code on your machine? As I have already mentioned, without setting GDAL_DATA, and without providing PROJECTION string gdal17.dll was able to create the output that can be viewed in our elevation data viewer of course, I get an error - ERROR 6: No translation for to PROJ.4 format.

[Kyle Shannon]

harishkumargvs, I get the same error with GDAL_DATA set properly. Unfortuneately, I don't have the time to investigate further.

[harishkumargvs]

Kyle, thank you for confirming that you too get the same error. May I know that you are the only person responsible for this module? I am not sure if changing the Priority of this ticket would help me have someone to investigate this issue but I am changing it as "highest" priority as this issue is a release blocker for us...

[Even Rouault]

harishkumargvs,

I've tested your code with GDAL 1.7 branch, and I also get the "ERROR 1: Too many points (441 out of 441) failed to transform,". I can't understand how it can work for you

The help page of the USGSDEM driver mentions "Input data must already be sampled in a geographic or UTM coordinate system.", and indeed when checking the source code of the driver, I can see it is not ready to receive an input dataset with other kind of coordinate system. So the error you get is logical.

You have 2 options :

• either provide a geotransform and SRS that are consistent with a geographic or UTM coordinate system
• or fix the driver to accept other coordinate systems. Regarding the later, I could be contracted to do such a fix. If you are interested, email me at contact at spatialys.com.

[Jukka Rahkonen]

Do you all agree that with the current title this ticket can be closed as invalid? A new ticket could be created as "Enhance the coordinate system support of USGSDEM driver".

[Jukka Rahkonen]

Closing, create a new ticket "Enhance the coordinate system support of USGSDEM driver" if needed.