Changeset 12293

Show
Ignore:
Timestamp:
10/02/07 01:41:35 (10 months ago)
Author:
warmerdam
Message:

backported 1.5 ENVI driver

Files:

Legend:

Unmodified
Added
Removed
Modified
Copied
Moved

  • sandbox/warmerdam/1.4-esri/gdal/frmts/raw/envidataset.cpp

    r10646 r12293  
    201201    } 
    202202 
    203     return 0; 
     203    return nESRIZone; // perhaps it *is* the USGS zone? 
    204204} 
    205205 
    …… 
    228228    int         ReadHeader( FILE * ); 
    229229    int         ProcessMapinfo( const char * ); 
     230    void        SetENVIDatum( OGRSpatialReference *, const char * ); 
     231    void        SetENVIEllipse( OGRSpatialReference *, char ** ); 
     232    void        WriteProjectionInfo(); 
    230233     
    231234    char        **SplitList( const char * ); 
    …… 
    409412/*      Write the rest of header.                                       */ 
    410413/* -------------------------------------------------------------------- */ 
    411     if ( pszProjection && !EQUAL(pszProjection, "") ) 
    412     { 
    413         const char      *pszHemisphere; 
    414         double          dfPixelY; 
    415         int             bNorth; 
    416         int             iUTMZone; 
    417         OGRSpatialReference oSRS; 
    418  
    419         char    *pszProj = pszProjection; 
    420  
    421         oSRS.importFromWkt( &pszProj ); 
    422         iUTMZone = oSRS.GetUTMZone( &bNorth ); 
    423         if ( iUTMZone ) 
    424         { 
    425             if ( bNorth ) 
    426             { 
    427                 pszHemisphere = "North"; 
    428                 dfPixelY = -adfGeoTransform[5]; 
    429             } 
    430             else 
    431             { 
    432                 pszHemisphere = "South"; 
    433                 dfPixelY = adfGeoTransform[5]; 
    434             } 
    435             VSIFPrintf( fp, "map info = {UTM, 1, 1, %f, %f, %f, %f, %d, %s}\n", 
    436                         adfGeoTransform[0], adfGeoTransform[3], adfGeoTransform[1], 
    437                         dfPixelY, iUTMZone, pszHemisphere); 
    438         } 
    439     } else { 
    440         // Suppose we are in North hemisphere. 
    441         double dfPixelY = -adfGeoTransform[5]; 
    442         const char* pszHemisphere = "North"; 
    443         VSIFPrintf( fp, "map info = {Unknown, 1, 1, %f, %f, %f, %f, %d, %s}\n", 
    444                     adfGeoTransform[0], adfGeoTransform[3], adfGeoTransform[1], 
    445                     dfPixelY, 0, pszHemisphere); 
    446     } 
    447      
    448  
     414    WriteProjectionInfo(); 
    449415    VSIFPrintf( fp, "band names = {\n" ); 
    450416    for ( int i = 1; i <= nBands; i++ ) 
    …… 
    462428 
    463429/************************************************************************/ 
     430/*                           GetEPSGGeogCS()                            */ 
     431/*                                                                      */ 
     432/*      Try to establish what the EPSG code for this coordinate         */ 
     433/*      systems GEOGCS might be.  Returns -1 if no reasonable guess     */ 
     434/*      can be made.                                                    */ 
     435/*                                                                      */ 
     436/*      TODO: We really need to do some name lookups.                   */ 
     437/************************************************************************/ 
     438 
     439static int ENVIGetEPSGGeogCS( OGRSpatialReference *poThis ) 
     440 
     441{ 
     442    const char *pszAuthName = poThis->GetAuthorityName( "GEOGCS" ); 
     443 
     444/* -------------------------------------------------------------------- */ 
     445/*      Do we already have it?                                          */ 
     446/* -------------------------------------------------------------------- */ 
     447    if( pszAuthName != NULL && EQUAL(pszAuthName,"epsg") ) 
     448        return atoi(poThis->GetAuthorityCode( "GEOGCS" )); 
     449 
     450/* -------------------------------------------------------------------- */ 
     451/*      Get the datum and geogcs names.                                 */ 
     452/* -------------------------------------------------------------------- */ 
     453    const char *pszGEOGCS = poThis->GetAttrValue( "GEOGCS" ); 
     454    const char *pszDatum = poThis->GetAttrValue( "DATUM" ); 
     455 
     456    // We can only operate on coordinate systems with a geogcs. 
     457    if( pszGEOGCS == NULL || pszDatum == NULL ) 
     458        return -1; 
     459 
     460/* -------------------------------------------------------------------- */ 
     461/*      Is this a "well known" geographic coordinate system?            */ 
     462/* -------------------------------------------------------------------- */ 
     463    int bWGS, bNAD; 
     464 
     465    bWGS = strstr(pszGEOGCS,"WGS") != NULL 
     466        || strstr(pszDatum, "WGS") 
     467        || strstr(pszGEOGCS,"World Geodetic System") 
     468        || strstr(pszGEOGCS,"World_Geodetic_System") 
     469        || strstr(pszDatum, "World Geodetic System") 
     470        || strstr(pszDatum, "World_Geodetic_System");  
     471 
     472    bNAD = strstr(pszGEOGCS,"NAD") != NULL 
     473        || strstr(pszDatum, "NAD") 
     474        || strstr(pszGEOGCS,"North American") 
     475        || strstr(pszGEOGCS,"North_American") 
     476        || strstr(pszDatum, "North American") 
     477        || strstr(pszDatum, "North_American");  
     478 
     479    if( bWGS && (strstr(pszGEOGCS,"84") || strstr(pszDatum,"84")) ) 
     480        return 4326; 
     481 
     482    if( bWGS && (strstr(pszGEOGCS,"72") || strstr(pszDatum,"72")) ) 
     483        return 4322; 
     484 
     485    if( bNAD && (strstr(pszGEOGCS,"83") || strstr(pszDatum,"83")) ) 
     486        return 4269; 
     487 
     488    if( bNAD && (strstr(pszGEOGCS,"27") || strstr(pszDatum,"27")) ) 
     489        return 4267; 
     490 
     491/* -------------------------------------------------------------------- */ 
     492/*      If we know the datum, associate the most likely GCS with        */ 
     493/*      it.                                                             */ 
     494/* -------------------------------------------------------------------- */ 
     495    pszAuthName = poThis->GetAuthorityName( "GEOGCS|DATUM" ); 
     496 
     497    if( pszAuthName != NULL  
     498        && EQUAL(pszAuthName,"epsg")  
     499        && poThis->GetPrimeMeridian() == 0.0 ) 
     500    { 
     501        int nDatum = atoi(poThis->GetAuthorityCode("GEOGCS|DATUM")); 
     502         
     503        if( nDatum >= 6000 && nDatum <= 6999 ) 
     504            return nDatum - 2000; 
     505    } 
     506 
     507    return -1; 
     508} 
     509 
     510/************************************************************************/ 
     511/*                        WriteProjectionInfo()                         */ 
     512/************************************************************************/ 
     513 
     514void ENVIDataset::WriteProjectionInfo() 
     515 
     516{ 
     517/* -------------------------------------------------------------------- */ 
     518/*      Format the location (geotransform) portion of the map info      */ 
     519/*      line.                                                           */ 
     520/* -------------------------------------------------------------------- */ 
     521    CPLString   osLocation; 
     522 
     523    osLocation.Printf( "1, 1, %.15g, %.15g, %.15g, %.15g",  
     524                       adfGeoTransform[0], adfGeoTransform[3],  
     525                       adfGeoTransform[1], fabs(adfGeoTransform[5]) ); 
     526                        
     527/* -------------------------------------------------------------------- */ 
     528/*      Minimal case - write out simple geotransform if we have a       */ 
     529/*      non-default geotransform.                                       */ 
     530/* -------------------------------------------------------------------- */ 
     531    if( pszProjection == NULL || strlen(pszProjection) == 0 ) 
     532    { 
     533        if( adfGeoTransform[0] != 0.0 || adfGeoTransform[1] != 1.0 
     534            || adfGeoTransform[2] != 0.0 || adfGeoTransform[3] != 0.0 
     535            || adfGeoTransform[4] != 0.0 || adfGeoTransform[5] != 1.0 ) 
     536        { 
     537            const char* pszHemisphere = "North"; 
     538            VSIFPrintf( fp, "map info = {Unknown, %s, %d, %s}\n", 
     539                        osLocation.c_str(), 0, pszHemisphere); 
     540        } 
     541        return; 
     542    } 
     543 
     544/* -------------------------------------------------------------------- */ 
     545/*      Ingest WKT.                                                     */ 
     546/* -------------------------------------------------------------------- */ 
     547    OGRSpatialReference oSRS; 
     548     
     549    char        *pszProj = pszProjection; 
     550     
     551    if( oSRS.importFromWkt( &pszProj ) != OGRERR_NONE ) 
     552        return; 
     553 
     554/* -------------------------------------------------------------------- */ 
     555/*      Try to translate the datum and get major/minor ellipsoid        */ 
     556/*      values.                                                         */ 
     557/* -------------------------------------------------------------------- */ 
     558    int nEPSG_GCS = ENVIGetEPSGGeogCS( &oSRS ); 
     559    CPLString osDatum, osCommaDatum; 
     560    double dfA, dfB; 
     561 
     562    if( nEPSG_GCS == 4326 ) 
     563        osDatum = "WGS-84"; 
     564    else if( nEPSG_GCS == 4322 ) 
     565        osDatum = "WGS-72"; 
     566    else if( nEPSG_GCS == 4269 ) 
     567        osDatum = "North America 1983"; 
     568    else if( nEPSG_GCS == 4267 ) 
     569        osDatum = "North America 1927"; 
     570    else if( nEPSG_GCS == 4230 ) 
     571        osDatum = "European 1950"; 
     572    else if( nEPSG_GCS == 4277 ) 
     573        osDatum = "Ordnance Survey of Great Britain '36"; 
     574    else if( nEPSG_GCS == 4291 ) 
     575        osDatum = "SAD-69/Brazil"; 
     576    else if( nEPSG_GCS == 4283 ) 
     577        osDatum = "Geocentric Datum of Australia 1994"; 
     578    else if( nEPSG_GCS == 4275 ) 
     579        osDatum = "Nouvelle Triangulation Francaise IGN"; 
     580 
     581    if( osDatum != "" ) 
     582        osCommaDatum.Printf( ",%s", osDatum.c_str() ); 
     583 
     584    dfA = oSRS.GetSemiMajor(); 
     585    dfB = oSRS.GetSemiMinor(); 
     586 
     587/* -------------------------------------------------------------------- */ 
     588/*      Do we have unusual linear units?                                */ 
     589/* -------------------------------------------------------------------- */ 
     590    CPLString osOptionalUnits; 
     591    if( fabs(oSRS.GetLinearUnits()-0.3048) < 0.0001 ) 
     592        osOptionalUnits = ", units=Feet"; 
     593 
     594/* -------------------------------------------------------------------- */ 
     595/*      Handle UTM case.                                                */ 
     596/* -------------------------------------------------------------------- */ 
     597    const char  *pszHemisphere; 
     598    const char  *pszProjName = oSRS.GetAttrValue("PROJECTION"); 
     599    int         bNorth; 
     600    int         iUTMZone; 
     601 
     602    iUTMZone = oSRS.GetUTMZone( &bNorth ); 
     603    if ( iUTMZone ) 
     604    { 
     605        if ( bNorth ) 
     606            pszHemisphere = "North"; 
     607        else 
     608            pszHemisphere = "South"; 
     609 
     610        VSIFPrintf( fp, "map info = {UTM, %s, %d, %s%s%s}\n", 
     611                    osLocation.c_str(), iUTMZone, pszHemisphere, 
     612                    osCommaDatum.c_str(), osOptionalUnits.c_str() ); 
     613    } 
     614    else if( oSRS.IsGeographic() ) 
     615    { 
     616        VSIFPrintf( fp, "map info = {Geographic Lat/Lon, %s%s}\n", 
     617                    osLocation.c_str(), osCommaDatum.c_str()); 
     618    } 
     619    else if( pszProjName == NULL ) 
     620    { 
     621        // what to do?  
     622    } 
     623    else if( EQUAL(pszProjName,SRS_PT_NEW_ZEALAND_MAP_GRID) ) 
     624    { 
     625        VSIFPrintf( fp, "map info = {New Zealand Map Grid, %s%s%s}\n", 
     626                    osLocation.c_str(),  
     627                    osCommaDatum.c_str(), osOptionalUnits.c_str() ); 
     628 
     629        VSIFPrintf( fp, "projection info = {39, %.16g, %.16g, %.16g, %.16g, %.16g, %.16g%s, New Zealand Map Grid}\n", 
     630                    dfA, dfB,  
     631                    oSRS.GetNormProjParm(SRS_PP_LATITUDE_OF_ORIGIN,0.0), 
     632                    oSRS.GetNormProjParm(SRS_PP_CENTRAL_MERIDIAN,0.0), 
     633                    oSRS.GetNormProjParm(SRS_PP_FALSE_EASTING,0.0), 
     634                    oSRS.GetNormProjParm(SRS_PP_FALSE_NORTHING,0.0), 
     635                    osCommaDatum.c_str() ); 
     636    } 
     637    else if( EQUAL(pszProjName,SRS_PT_TRANSVERSE_MERCATOR) ) 
     638    { 
     639        VSIFPrintf( fp, "map info = {Transverse Mercator, %s%s%s}\n", 
     640                    osLocation.c_str(),  
     641                    osCommaDatum.c_str(), osOptionalUnits.c_str() ); 
     642 
     643        VSIFPrintf( fp, "projection info = {3, %.16g, %.16g, %.16g, %.16g, %.16g, %.16g, %.16g%s, Transverse Mercator}\n", 
     644                    dfA, dfB,  
     645                    oSRS.GetNormProjParm(SRS_PP_LATITUDE_OF_ORIGIN,0.0), 
     646                    oSRS.GetNormProjParm(SRS_PP_CENTRAL_MERIDIAN,0.0), 
     647                    oSRS.GetNormProjParm(SRS_PP_FALSE_EASTING,0.0), 
     648                    oSRS.GetNormProjParm(SRS_PP_FALSE_NORTHING,0.0), 
     649                    oSRS.GetNormProjParm(SRS_PP_SCALE_FACTOR,1.0), 
     650                    osCommaDatum.c_str() ); 
     651    } 
     652    else if( EQUAL(pszProjName,SRS_PT_LAMBERT_CONFORMAL_CONIC_2SP) 
     653         || EQUAL(pszProjName,SRS_PT_LAMBERT_CONFORMAL_CONIC_2SP_BELGIUM) ) 
     654    { 
     655        VSIFPrintf( fp, "map info = {Lambert Conformal Conic, %s%s%s}\n", 
     656                    osLocation.c_str(),  
     657                    osCommaDatum.c_str(), osOptionalUnits.c_str() ); 
     658 
     659        VSIFPrintf( fp, "projection info = {4, %.16g, %.16g, %.16g, %.16g, %.16g, %.16g, %.16g, %.16g%s, Lambert Conformal Conic}\n", 
     660                    dfA, dfB,  
     661                    oSRS.GetNormProjParm(SRS_PP_LATITUDE_OF_ORIGIN,0.0), 
     662                    oSRS.GetNormProjParm(SRS_PP_CENTRAL_MERIDIAN,0.0), 
     663                    oSRS.GetNormProjParm(SRS_PP_FALSE_EASTING,0.0), 
     664                    oSRS.GetNormProjParm(SRS_PP_FALSE_NORTHING,0.0), 
     665                    oSRS.GetNormProjParm(SRS_PP_STANDARD_PARALLEL_1,0.0), 
     666                    oSRS.GetNormProjParm(SRS_PP_STANDARD_PARALLEL_2,0.0), 
     667                    osCommaDatum.c_str() ); 
     668    } 
     669    else if( EQUAL(pszProjName, 
     670                   SRS_PT_HOTINE_OBLIQUE_MERCATOR_TWO_POINT_NATURAL_ORIGIN) ) 
     671    { 
     672        VSIFPrintf( fp, "map info = {Hotine Oblique Mercator A, %s%s%s}\n", 
     673                    osLocation.c_str(),  
     674                    osCommaDatum.c_str(), osOptionalUnits.c_str() ); 
     675 
     676        VSIFPrintf( fp, "projection info = {5, %.16g, %.16g, %.16g, %.16g, %.16g, %.16g, %.16g, %.16g, %.16g, %.16g%s, Hotine Oblique Mercator A}\n", 
     677                    dfA, dfB,  
     678                    oSRS.GetNormProjParm(SRS_PP_LATITUDE_OF_ORIGIN,0.0), 
     679                    oSRS.GetNormProjParm(SRS_PP_LATITUDE_OF_POINT_1,0.0), 
     680                    oSRS.GetNormProjParm(SRS_PP_LONGITUDE_OF_POINT_1,0.0), 
     681                    oSRS.GetNormProjParm(SRS_PP_LATITUDE_OF_POINT_2,0.0), 
     682                    oSRS.GetNormProjParm(SRS_PP_LONGITUDE_OF_POINT_2,0.0), 
     683                    oSRS.GetNormProjParm(SRS_PP_FALSE_EASTING,0.0), 
     684                    oSRS.GetNormProjParm(SRS_PP_FALSE_NORTHING,0.0), 
     685                    oSRS.GetNormProjParm(SRS_PP_SCALE_FACTOR,1.0), 
     686                    osCommaDatum.c_str() ); 
     687    } 
     688    else if( EQUAL(pszProjName,SRS_PT_HOTINE_OBLIQUE_MERCATOR) ) 
     689    { 
     690        VSIFPrintf( fp, "map info = {Hotine Oblique Mercator B, %s%s%s}\n", 
     691                    osLocation.c_str(),  
     692                    osCommaDatum.c_str(), osOptionalUnits.c_str() ); 
     693 
     694        VSIFPrintf( fp, "projection info = {6, %.16g, %.16g, %.16g, %.16g, %.16g, %.16g, %.16g, %.16g%s, Hotine Oblique Mercator B}\n", 
     695                    dfA, dfB,  
     696                    oSRS.GetNormProjParm(SRS_PP_LATITUDE_OF_ORIGIN,0.0), 
     697                    oSRS.GetNormProjParm(SRS_PP_CENTRAL_MERIDIAN,0.0), 
     698                    oSRS.GetNormProjParm(SRS_PP_AZIMUTH,0.0), 
     699                    oSRS.GetNormProjParm(SRS_PP_FALSE_EASTING,0.0), 
     700                    oSRS.GetNormProjParm(SRS_PP_FALSE_NORTHING,0.0), 
     701                    oSRS.GetNormProjParm(SRS_PP_SCALE_FACTOR,1.0), 
     702                    osCommaDatum.c_str() ); 
     703    } 
     704    else if( EQUAL(pszProjName,SRS_PT_STEREOGRAPHIC)  
     705             || EQUAL(pszProjName,SRS_PT_OBLIQUE_STEREOGRAPHIC) ) 
     706    { 
     707        VSIFPrintf( fp, "map info = {Stereographic (ellipsoid), %s%s%s}\n", 
     708                    osLocation.c_str(),  
     709                    osCommaDatum.c_str(), osOptionalUnits.c_str() ); 
     710 
     711        VSIFPrintf( fp, "projection info = {7, %.16g, %.16g, %.16g, %.16g, %.16g, %.16g, %.16g, %.16g%s, Stereographic (ellipsoid)}\n", 
     712                    dfA, dfB,  
     713                    oSRS.GetNormProjParm(SRS_PP_LATITUDE_OF_ORIGIN,0.0), 
     714                    oSRS.GetNormProjParm(SRS_PP_CENTRAL_MERIDIAN,0.0), 
     715                    oSRS.GetNormProjParm(SRS_PP_FALSE_EASTING,0.0), 
     716                    oSRS.GetNormProjParm(SRS_PP_FALSE_NORTHING,0.0), 
     717                    oSRS.GetNormProjParm(SRS_PP_SCALE_FACTOR,1.0), 
     718                    osCommaDatum.c_str() ); 
     719    } 
     720    else if( EQUAL(pszProjName,SRS_PT_ALBERS_CONIC_EQUAL_AREA) ) 
     721    { 
     722        VSIFPrintf( fp, "map info = {Albers Conical Equal Area, %s%s%s}\n", 
     723                    osLocation.c_str(),  
     724                    osCommaDatum.c_str(), osOptionalUnits.c_str() ); 
     725 
     726        VSIFPrintf( fp, "projection info = {9, %.16g, %.16g, %.16g, %.16g, %.16g, %.16g, %.16g, %.16g%s, Albers Conical Equal Area}\n", 
     727                    dfA, dfB,  
     728                    oSRS.GetNormProjParm(SRS_PP_LATITUDE_OF_ORIGIN,0.0), 
     729                    oSRS.GetNormProjParm(SRS_PP_CENTRAL_MERIDIAN,0.0), 
     730                    oSRS.GetNormProjParm(SRS_PP_FALSE_EASTING,0.0), 
     731                    oSRS.GetNormProjParm(SRS_PP_FALSE_NORTHING,0.0), 
     732                    oSRS.GetNormProjParm(SRS_PP_STANDARD_PARALLEL_1,0.0), 
     733                    oSRS.GetNormProjParm(SRS_PP_STANDARD_PARALLEL_2,0.0), 
     734                    osCommaDatum.c_str() ); 
     735    } 
     736    else if( EQUAL(pszProjName,SRS_PT_POLYCONIC) ) 
     737    { 
     738        VSIFPrintf( fp, "map info = {Polyconic, %s%s%s}\n", 
     739                    osLocation.c_str(),  
     740                    osCommaDatum.c_str(), osOptionalUnits.c_str() ); 
     741 
     742        VSIFPrintf( fp, "projection info = {10, %.16g, %.16g, %.16g, %.16g, %.16g, %.16g%s, Polyconic}\n", 
     743                    dfA, dfB,  
     744                    oSRS.GetNormProjParm(SRS_PP_LATITUDE_OF_ORIGIN,0.0), 
     745                    oSRS.GetNormProjParm(SRS_PP_CENTRAL_MERIDIAN,0.0), 
     746                    oSRS.GetNormProjParm(SRS_PP_FALSE_EASTING,0.0), 
     747                    oSRS.GetNormProjParm(SRS_PP_FALSE_NORTHING,0.0), 
     748                    osCommaDatum.c_str() ); 
     749    } 
     750    else if( EQUAL(pszProjName,SRS_PT_LAMBERT_AZIMUTHAL_EQUAL_AREA) ) 
     751    { 
     752        VSIFPrintf( fp, "map info = {Lambert Azimuthal Equal Area, %s%s%s}\n", 
     753                    osLocation.c_str(),  
     754                    osCommaDatum.c_str(), osOptionalUnits.c_str() ); 
     755 
     756        VSIFPrintf( fp, "projection info = {11, %.16g, %.16g, %.16g, %.16g, %.16g, %.16g%s, Lambert Azimuthal Equal Area}\n", 
     757                    dfA, dfB,  
     758                    oSRS.GetNormProjParm(SRS_PP_LATITUDE_OF_ORIGIN,0.0), 
     759                    oSRS.GetNormProjParm(SRS_PP_CENTRAL_MERIDIAN,0.0), 
     760                    oSRS.GetNormProjParm(SRS_PP_FALSE_EASTING,0.0), 
     761                    oSRS.GetNormProjParm(SRS_PP_FALSE_NORTHING,0.0), 
     762                    osCommaDatum.c_str() ); 
     763    } 
     764    else if( EQUAL(pszProjName,SRS_PT_AZIMUTHAL_EQUIDISTANT) ) 
     765    { 
     766        VSIFPrintf( fp, "map info = {Azimuthal Equadistant, %s%s%s}\n", 
     767                    osLocation.c_str(),  
     768                    osCommaDatum.c_str(), osOptionalUnits.c_str() ); 
     769 
     770        VSIFPrintf( fp, "projection info = {12, %.16g, %.16g, %.16g, %.16g, %.16g, %.16g%s, Azimuthal Equadistant}\n", 
     771                    dfA, dfB,  
     772                    oSRS.GetNormProjParm(SRS_PP_LATITUDE_OF_ORIGIN,0.0), 
     773                    oSRS.GetNormProjParm(SRS_PP_CENTRAL_MERIDIAN,0.0), 
     774                    oSRS.GetNormProjParm(SRS_PP_FALSE_EASTING,0.0), 
     775                    oSRS.GetNormProjParm(SRS_PP_FALSE_NORTHING,0.0), 
     776                    osCommaDatum.c_str() ); 
     777    } 
     778    else if( EQUAL(pszProjName,SRS_PT_POLAR_STEREOGRAPHIC) ) 
     779    { 
     780        VSIFPrintf( fp, "map info = {Polar Stereographic, %s%s%s}\n", 
     781                    osLocation.c_str(),  
     782                    osCommaDatum.c_str(), osOptionalUnits.c_str() ); 
     783 
     784        VSIFPrintf( fp, "projection info = {31, %.16g, %.16g, %.16g, %.16g, %.16g, %.16g%s, Polar Stereographic}\n", 
     785                    dfA, dfB,  
     786                    oSRS.GetNormProjParm(SRS_PP_LATITUDE_OF_ORIGIN,90.0), 
     787                    oSRS.GetNormProjParm(SRS_PP_CENTRAL_MERIDIAN,0.0), 
     788                    oSRS.GetNormProjParm(SRS_PP_FALSE_EASTING,0.0), 
     789                    oSRS.GetNormProjParm(SRS_PP_FALSE_NORTHING,0.0), 
     790                    osCommaDatum.c_str() ); 
     791    } 
     792    else 
     793    { 
     794        VSIFPrintf( fp, "map info = {%s, %s}\n", 
     795                    pszProjName, osLocation.c_str()); 
     796    } 
     797} 
     798 
     799 
     800/************************************************************************/ 
    464801/*                          GetProjectionRef()                          */ 
    465802/************************************************************************/ 
    …… 
    568905 
    569906/************************************************************************/ 
     907/*                            SetENVIDatum()                            */ 
     908/************************************************************************/ 
     909 
     910void ENVIDataset::SetENVIDatum( OGRSpatialReference *poSRS,  
     911                                const char *pszENVIDatumName ) 
     912 
     913{ 
     914    // datums 
     915    if( EQUAL(pszENVIDatumName, "WGS-84") ) 
     916        poSRS->SetWellKnownGeogCS( "WGS84" ); 
     917    else if( EQUAL(pszENVIDatumName, "WGS-72") ) 
     918        poSRS->SetWellKnownGeogCS( "WGS72" ); 
     919    else if( EQUAL(pszENVIDatumName, "North America 1983") ) 
     920        poSRS->SetWellKnownGeogCS( "NAD83" ); 
     921    else if( EQUAL(pszENVIDatumName, "North America 1927")  
     922             || strstr(pszENVIDatumName,"NAD27")  
     923             || strstr(pszENVIDatumName,"NAD-27") ) 
     924        poSRS->SetWellKnownGeogCS( "NAD27" ); 
     925    else if( EQUALN(pszENVIDatumName, "European 1950",13) ) 
     926        poSRS->SetWellKnownGeogCS( "EPSG:4230" ); 
     927    else if( EQUAL(pszENVIDatumName, "Ordnance Survey of Great Britain '36") ) 
     928        poSRS->SetWellKnownGeogCS( "EPSG:4277" ); 
     929    else if( EQUAL(pszENVIDatumName, "SAD-69/Brazil") ) 
     930        poSRS->SetWellKnownGeogCS( "EPSG:4291" ); 
     931    else if( EQUAL(pszENVIDatumName, "Geocentric Datum of Australia 1994") ) 
     932        poSRS->SetWellKnownGeogCS( "EPSG:4283" ); 
     933    else if( EQUAL(pszENVIDatumName, "Nouvelle Triangulation Francaise IGN") ) 
     934        poSRS->SetWellKnownGeogCS( "EPSG:4275" ); 
     935 
     936    // Ellipsoids 
     937    else if( EQUAL(pszENVIDatumName, "GRS 80") ) 
     938        poSRS->SetWellKnownGeogCS( "NAD83" ); 
     939    else if( EQUAL(pszENVIDatumName, "Airy") ) 
     940        poSRS->SetWellKnownGeogCS( "EPSG:4001" ); 
     941    else if( EQUAL(pszENVIDatumName, "Australian National") ) 
     942        poSRS->SetWellKnownGeogCS( "EPSG:4003" ); 
     943    else if( EQUAL(pszENVIDatumName, "Bessel 1841") ) 
     944        poSRS->SetWellKnownGeogCS( "EPSG:4004" ); 
     945    else if( EQUAL(pszENVIDatumName, "Clark 1866") ) 
     946        poSRS->SetWellKnownGeogCS( "EPSG:4008" ); 
     947    else  
     948    { 
     949        CPLError( CE_Warning, CPLE_AppDefined, 
     950                  "Unrecognised datum '%s', defaulting to WGS84." ); 
     951        poSRS->SetWellKnownGeogCS( "WGS84" ); 
     952    } 
     953} 
     954 
     955/************************************************************************/ 
     956/*                           SetENVIEllipse()                           */ 
     957/************************************************************************/ 
     958 
     959void ENVIDataset::SetENVIEllipse( OGRSpatialReference *poSRS,  
     960                                  char **papszPI_EI ) 
     961 
     962{ 
     963    double dfA = CPLAtofM(papszPI_EI[0]); 
     964    double dfB = CPLAtofM(papszPI_EI[1]); 
     965    double dfInvF; 
     966 
     967    if( fabs(dfA-dfB) < 0.1 ) 
     968        dfInvF = 0.0; // sphere 
     969    else 
     970        dfInvF = dfA / (dfA - dfB); 
     971     
     972     
     973    poSRS->SetGeogCS( "Ellipse Based", "Ellipse Based", "Unnamed",  
     974                      dfA, dfInvF ); 
     975} 
     976 
     977/************************************************************************/ 
    570978/*                           ProcessMapinfo()                           */ 
    571979/*                                                                      */ 
    …… 
    590998    } 
    591999 
     1000/* -------------------------------------------------------------------- */ 
     1001/*      Check if we have projection info, and if so parse it.           */ 
     1002/* -------------------------------------------------------------------- */ 
     1003    char        **papszPI = NULL; 
     1004    int         nPICount = 0; 
     1005    if( CSLFetchNameValue( papszHeader, "projection_info" ) != NULL ) 
     1006    { 
     1007        papszPI = SplitList(  
     1008            CSLFetchNameValue( papszHeader, "projection_info" ) ); 
     1009        nPICount = CSLCount(papszPI); 
     1010    } 
     1011 
     1012/* -------------------------------------------------------------------- */ 
     1013/*      Capture geotransform.                                           */ 
     1014/* -------------------------------------------------------------------- */ 
    5921015    adfGeoTransform[1] = atof(papszFields[5]);      // Pixel width 
    5931016    adfGeoTransform[5] = -atof(papszFields[6]);     // Pixel height 
    …… 
    5991022    adfGeoTransform[4] = 0.0; 
    6001023 
     1024/* -------------------------------------------------------------------- */ 
     1025/*      Capture projection.                                             */ 
     1026/* -------------------------------------------------------------------- */ 
    6011027    if( EQUALN(papszFields[0],"UTM",3) && nCount >= 9 ) 
    6021028    { 
    6031029        oSRS.SetUTM( atoi(papszFields[7]),  
    6041030                     !EQUAL(papszFields[8],"South") ); 
    605         oSRS.SetWellKnownGeogCS( "WGS84" ); 
     1031        if( nCount >= 10 && strstr(papszFields[9],"=") == NULL ) 
     1032            SetENVIDatum( &oSRS, papszFields[9] ); 
     1033        else 
     1034            oSRS.SetWellKnownGeogCS( "WGS84" ); 
    6061035    } 
    6071036    else if( EQUALN(papszFields[0],"State Plane (NAD 27)",19) 
    608              && nCount >= 8 ) 
     1037             && nCount >= 7 ) 
    6091038    { 
    6101039        oSRS.SetStatePlane( ESRIToUSGSZone(atoi(papszFields[7])), FALSE ); 
    6111040    } 
    6121041    else if( EQUALN(papszFields[0],"State Plane (NAD 83)",19) 
    613              && nCount >= 8 ) 
     1042             && nCount >= 7 ) 
    6141043    { 
    6151044        oSRS.SetStatePlane( ESRIToUSGSZone(atoi(papszFields[7])), TRUE ); 
    …… 
    6181047             && nCount >= 8 ) 
    6191048    { 
    620         oSRS.SetWellKnownGeogCS( "WGS84" ); 
    621     } 
    622  
     1049        if( nCount >= 8 && strstr(papszFields[7],"=") == NULL ) 
     1050            SetENVIDatum( &oSRS, papszFields[7] ); 
     1051        else 
     1052            oSRS.SetWellKnownGeogCS( "WGS84" ); 
     1053    } 
     1054    else if( nPICount > 8 && atoi(papszPI[0]) == 3 ) // TM 
     1055    { 
     1056        oSRS.SetTM( CPLAtofM(papszPI[3]), CPLAtofM(papszPI[4]), 
     1057                    CPLAtofM(papszPI[7]), 
     1058                    CPLAtofM(papszPI[5]), CPLAtofM(papszPI[6]) ); 
     1059    } 
     1060    else if( nPICount > 8 && atoi(papszPI[0]) == 4 ) // Lambert Conformal Conic 
     1061    { 
     1062        oSRS.SetLCC( CPLAtofM(papszPI[7]), CPLAtofM(papszPI[8]), 
     1063                     CPLAtofM(papszPI[3]), CPLAtofM(papszPI[4]), 
     1064                     CPLAtofM(papszPI[5]), CPLAtofM(papszPI[6]) ); 
     1065    } 
     1066    else if( nPICount > 10 && atoi(papszPI[0]) == 5 ) // Oblique Merc (2 point) 
     1067    { 
     1068        oSRS.SetHOM2PNO( CPLAtofM(papszPI[3]),  
     1069                         CPLAtofM(papszPI[4]), CPLAtofM(papszPI[5]),  
     1070                         CPLAtofM(papszPI[6]), CPLAtofM(papszPI[7]),  
     1071                         CPLAtofM(papszPI[10]),  
     1072                         CPLAtofM(papszPI[8]), CPLAtofM(papszPI[9]) ); 
     1073    } 
     1074    else if( nPICount > 8 && atoi(papszPI[0]) == 6 ) // Oblique Merc  
     1075    { 
     1076        oSRS.SetHOM(CPLAtofM(papszPI[3]), CPLAtofM(papszPI[4]),  
     1077                    CPLAtofM(papszPI[5]), 0.0, 
     1078                    CPLAtofM(papszPI[8]),  
     1079                    CPLAtofM(papszPI[6]), CPLAtofM(papszPI[7]) ); 
     1080    } 
     1081    else if( nPICount > 8 && atoi(papszPI[0]) == 7 ) // Stereographic 
     1082    { 
     1083        oSRS.SetStereographic( CPLAtofM(papszPI[3]), CPLAtofM(papszPI[4]), 
     1084                               CPLAtofM(papszPI[7]), 
     1085                               CPLAtofM(papszPI[5]), CPLAtofM(papszPI[6]) ); 
     1086    } 
     1087    else if( nPICount > 8 && atoi(papszPI[0]) == 9 ) // Albers Equal Area 
     1088    { 
     1089        oSRS.SetACEA( CPLAtofM(papszPI[7]), CPLAtofM(papszPI[8]), 
     1090                      CPLAtofM(papszPI[3]), CPLAtofM(papszPI[4]), 
     1091                      CPLAtofM(papszPI[5]), CPLAtofM(papszPI[6]) ); 
     1092    } 
     1093    else if( nPICount > 6 && atoi(papszPI[0]) == 10 ) // Polyconic 
     1094    { 
     1095        oSRS.SetPolyconic(CPLAtofM(papszPI[3]), CPLAtofM(papszPI[4]),  
     1096                          CPLAtofM(papszPI[5]), CPLAtofM(papszPI[6]) ); 
     1097    } 
     1098    else if( nPICount > 6 && atoi(papszPI[0]) == 11 ) // LAEA 
     1099    { 
     1100        oSRS.SetLAEA(CPLAtofM(papszPI[3]), CPLAtofM(papszPI[4]),  
     1101                     CPLAtofM(papszPI[5]), CPLAtofM(papszPI[6]) ); 
     1102    } 
     1103    else if( nPICount > 6 && atoi(papszPI[0]) == 12 ) // Azimuthal Equid. 
     1104    { 
     1105        oSRS.SetAE(CPLAtofM(papszPI[3]), CPLAtofM(papszPI[4]),  
     1106                   CPLAtofM(papszPI[5]), CPLAtofM(papszPI[6]) ); 
     1107    } 
     1108    else if( nPICount > 6 && atoi(papszPI[0]) == 31 ) // Polar Stereographic 
     1109    { 
     1110        oSRS.SetPS(CPLAtofM(papszPI[3]), CPLAtofM(papszPI[4]),  
     1111                   1.0, 
     1112                   CPLAtofM(papszPI[5]), CPLAtofM(papszPI[6]) ); 
     1113    } 
     1114 
     1115    // Still lots more that could be added for someone with the patience. 
     1116 
     1117/* -------------------------------------------------------------------- */ 
     1118/*      fallback to localcs if we don't recognise things.               */ 
     1119/* -------------------------------------------------------------------- */ 
    6231120    if( oSRS.GetRoot() == NULL ) 
    6241121        oSRS.SetLocalCS( papszFields[0] ); 
    6251122 
     1123/* -------------------------------------------------------------------- */ 
     1124/*      Try to set datum from projection info line if we have a         */ 
     1125/*      projected coordinate system without a GEOGCS.                   */ 
     1126/* -------------------------------------------------------------------- */ 
     1127    if( oSRS.IsProjected() && oSRS.GetAttrNode("GEOGCS") == NULL  
     1128        && nPICount > 3 ) 
     1129    { 
     1130        // Do we have a datum on the projection info line? 
     1131        int iDatum = nPICount-1; 
     1132 
     1133        // Ignore units= items. 
     1134        if( strstr(papszPI[iDatum],"=") != NULL ) 
     1135            iDatum--; 
     1136 
     1137        // Skip past the name.  
     1138        iDatum--; 
     1139     
     1140        CPLString osDatumName = papszPI[iDatum]; 
     1141        if( osDatumName.find_first_of("abcdefghijklmnopqrstuvwxyz" 
     1142                                      "ABCDEFGHIJKLMNOPQRSTUVWXYZ")  
     1143                                      != std::string::npos ) 
     1144        { 
     1145            SetENVIDatum( &oSRS, osDatumName ); 
     1146        } 
     1147        else 
     1148        { 
     1149            SetENVIEllipse( &oSRS, papszPI + 1 ); 
     1150        } 
     1151    } 
     1152 
     1153 
     1154/* -------------------------------------------------------------------- */ 
     1155/*      Try to process specialized units.                               */ 
     1156/* -------------------------------------------------------------------- */ 
    6261157    if( EQUAL(papszFields[nCount-1],"units=Feet") ) 
    6271158    { 
    …… 
    6401171    } 
    6411172 
     1173/* -------------------------------------------------------------------- */ 
     1174/*      Turn back into WKT.                                             */ 
     1175/* -------------------------------------------------------------------- */ 
    6421176    if( oSRS.GetRoot() != NULL ) 
    6431177    { 
    …… 
    6521186 
    6531187    CSLDestroy( papszFields ); 
     1188    CSLDestroy( papszPI ); 
    6541189    return TRUE; 
    6551190} 
    …… 
    7561291/* -------------------------------------------------------------------- */ 
    7571292    const char  *pszMode; 
    758     const char *pszHdrFilename; 
     1293    CPLString   osHdrFilename; 
    7591294    FILE        *fpHeader; 
    7601295 
    …… 
    7641299        pszMode = "r"; 
    7651300 
    766     pszHdrFilename = CPLResetExtension( poOpenInfo->pszFilename, "hdr" ); 
    767     fpHeader = VSIFOpen( pszHdrFilename, pszMode ); 
     1301    osHdrFilename = CPLResetExtension( poOpenInfo->pszFilename, "hdr" ); 
     1302    fpHeader = VSIFOpen( osHdrFilename, pszMode ); 
    7681303 
    7691304#ifndef WIN32 
    7701305    if( fpHeader == NULL ) 
    7711306    { 
    772         pszHdrFilename = CPLResetExtension( poOpenInfo->pszFilename, "HDR" ); 
    773         fpHeader = VSIFOpen( pszHdrFilename, pszMode ); 
     1307        osHdrFilename = CPLResetExtension( poOpenInfo->pszFilename, "HDR" ); 
     1308        fpHeader = VSIFOpen( osHdrFilename, pszMode ); 
    7741309    } 
    7751310#endif 
    7761311    if( fpHeader == NULL ) 
    7771312    { 
    778         pszHdrFilename = CPLFormFilename( NULL, poOpenInfo->pszFilename,  
     1313        osHdrFilename = CPLFormFilename( NULL, poOpenInfo->pszFilename,  
    7791314                                          "hdr" ); 
    780         fpHeader = VSIFOpen( pszHdrFilename, pszMode ); 
     1315        fpHeader = VSIFOpen( osHdrFilename, pszMode ); 
    7811316    } 
    7821317#ifndef WIN32 
    7831318    if( fpHeader == NULL ) 
    7841319    { 
    785         pszHdrFilename = CPLFormFilename( NULL, poOpenInfo->pszFilename,  
     1320        osHdrFilename = CPLFormFilename( NULL, poOpenInfo->pszFilename,  
    7861321                                          "HDR" ); 
    787         fpHeader = VSIFOpen( pszHdrFilename, pszMode ); 
     1322        fpHeader = VSIFOpen( osHdrFilename, pszMode ); 
    7881323    } 
    7891324#endif 
    …… 
    8141349 
    8151350    poDS = new ENVIDataset(); 
    816     poDS->pszHDRFilename = pszHdrFilename; 
     1351    poDS->pszHDRFilename = CPLStrdup(osHdrFilename); 
    8171352    poDS->fp = fpHeader; 
    8181353 
    …… 
    8311366    if( EQUAL(CPLGetExtension(poOpenInfo->pszFilename), "hdr") ) 
    8321367    { 
     1368        delete poDS; 
    8331369        CPLError( CE_Failure, CPLE_AppDefined,  
    8341370                  "The selected file is an ENVI header file, but to\n"