Why are the output point shapefile PRJ files empty?

Question

I have tried everything I could find (using C++) to create a new .prj file for an output point shapefile using OGR. While I get no compilation or execution errors and ArcMap accepts the file as a point file, ArcMap still claims that my output shapefile has no spatial reference:

Data Type:  Shapefile Feature Class 
Shapefile:  G:\GIS\Working\FIrani\Ubuntu\0_From_Ubuntu\Seeds001.shp
Geometry Type:  Point
Coordinates have Z values:  No 
Coordinates have measures:  No 

Coordinate System:  <Undefined>

How do I get OGR to actually apply my ACEA specified coordinate information to its output shapefile? Do I actually have to delete the current .prj and create and write out a new .prj file myself?

My current (kitchen sink) approach is as follows:

    /*
    Define projection for output seed file.
    --------------------------------------- */
    char *pszWkt = NULL;
    oSRS.SetProjection("ACEA");
    oSRS.SetProjCS("ACEA");
    oSRS.SetWellKnownGeogCS("ACEA");
    oSRS.SetACEA(29.5, 45.5, 23.0, -96.0, 0.0, 0.0); 
    oSRS.morphToESRI();
    oSRS.importFromWkt( &pszWkt );
  pioSeedLayer = pioSeedDS->CreateLayer( pszSeedsfn, &oSRS, wkbPoint, NULL ); 
  // Ouput projection via oSRS.
  if( pioSeedLayer == NULL )
  {
    printf( "Layer creation for seeds failed.\n" );
    exit( 1 );
  }

Answer 1

I recommend you look through http://www.gdal.org/ogr/osr_tutorial.html which discusses the how of your problem.

The other issue you have is that ACEA is not enough information to determine which projection your data is in. http://spatialreference.org/ref/?search=albers+equal will help you narrow down the projection you want (and give you the WKT string you need to create the object).

I think you just need to do something like:

char *pszWkt = "PROJCS["Albers Conical Equal Area",GEOGCS["GCS_North_American_1927",DATUM["D_North_American_1927",SPHEROID["Clarke_1866",6378206.4,294.9786982]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Albers"],PARAMETER["False_Easting",0.0],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",-105.6855555555556],PARAMETER["Standard_Parallel_1",40.38611111111113],PARAMETER["Standard_Parallel_2",47.27722222222222],PARAMETER["Latitude_Of_Origin",37.0],UNIT["Meter",1.0]]";
oSRS.importFromWkt( &pszWkt );


pioSeedLayer = pioSeedDS->CreateLayer( pszSeedsfn, &oSRS, wkbPoint, NULL ); 
  // Ouput projection via oSRS.
  if( pioSeedLayer == NULL )
  {
    printf( "Layer creation for seeds failed.\n" );
    exit( 1 );
  }

Answer 2

Thanks for your Suggestion.

I think I found a nicer way to transfer GDAL projection Information to an OGR Spatial Reference. I ported a GDAL example code snippet written in C to C++ as below. It would be cleaner if I didn't have to call GetProjectionRef AND GDALGetProjectionRef, but the crude code snippet I wrote below works for me and I am moving on for now. Perhaps this will help someone. [I tried using block Quotes ">" below but they didn't work. I hope this is legible] :

OGRSpatialReference oSRS = NULL;

GDALDataset  *inputDS;

const char  *pszInput_filename;

/*

Open Input Probability Surface.

------------------------------- */

inputDS = (GDALDataset * ) GDALOpen( pszInput_filename, GA_ReadOnly );

if ( inputDS != NULL )

{

  /*

  Capture input Projection.

 ------------------------- */


if( inputDS->GetProjectionRef()  != NULL )
 {


pszWkt = GDALGetProjectionRef(inputDS);

if (pszWkt && pszWkt[0])

{

  oSRS = OGRSpatialReference(pszWkt);  // Pass GDAL Projection to OGR.
}

else

{

  printf( "Unable to get ProjectionReference for %s\n", pszInput_filename );

  exit(0);

   }    

}

else

{

printf( "Unable to get ProjectionReference for %s\n", pszInput_filename );

exit(0);

}