Changeset 13780

Show
Ignore:
Timestamp:
02/14/08 10:25:46 (5 months ago)
Author:
warmerdam
Message:

add support for projected NRL products (#2225)

Files:

Legend:

Unmodified
Added
Removed
Modified
Copied
Moved

  • branches/1.5/gdal/frmts/hdf4/hdf4imagedataset.cpp

    r13227 r13780  
    971971/*                       CaptureNRLGeoTransform()                       */ 
    972972/*                                                                      */ 
    973 /*      Capture geotransform and coordinate system from NRL (Navel      */ 
     973/*      Capture geotransform and coordinate system from NRL (Naval      */ 
    974974/*      Research Laboratory, Stennis Space Center) metadata.            */ 
    975975/************************************************************************/ 
    …… 
    10291029 
    10301030{ 
     1031/* -------------------------------------------------------------------- */ 
     1032/*      Collect the four corners.                                       */ 
     1033/* -------------------------------------------------------------------- */ 
    10311034    double adfXY[8]; 
    10321035    static const char *apszItems[] = { 
    10331036        "mapUpperLeft", "mapUpperRight", "mapLowerLeft", "mapLowerRight" }; 
    10341037    int iCorner; 
     1038    int bLLPossible = TRUE; 
    10351039 
    10361040    for( iCorner = 0; iCorner < 4; iCorner++ ) 
    …… 
    10501054        adfXY[iCorner*2+1] = CPLAtof( papszTokens[0] ); 
    10511055 
     1056        if( adfXY[iCorner*2+0] < -360 || adfXY[iCorner*2+0] > 360  
     1057            || adfXY[iCorner*2+1] < -90 || adfXY[iCorner*2+1] > 90 ) 
     1058            bLLPossible = FALSE; 
     1059 
    10521060        CSLDestroy( papszTokens ); 
    10531061    } 
    10541062 
    1055     if( adfXY[0*2+0] == adfXY[2*2+0] && adfXY[0*2+1] == adfXY[1*2+1] ) 
     1063/* -------------------------------------------------------------------- */ 
     1064/*      Does this look like nice clean "northup" lat/long data?         */ 
     1065/* -------------------------------------------------------------------- */ 
     1066    if( adfXY[0*2+0] == adfXY[2*2+0] && adfXY[0*2+1] == adfXY[1*2+1]  
     1067        && bLLPossible ) 
    10561068    { 
    10571069        bHasGeoTransform = TRUE; 
    …… 
    10661078        CPLFree( pszProjection ); 
    10671079        oSRS.exportToWkt( &pszProjection ); 
     1080    } 
     1081 
     1082/* -------------------------------------------------------------------- */ 
     1083/*      Can we find the USGS Projection Parameters?                     */ 
     1084/* -------------------------------------------------------------------- */ 
     1085    int  bGotGCTPProjection = FALSE; 
     1086    int  iSDSIndex = FAIL, iSDS = FAIL; 
     1087    const char *mapProjection = CSLFetchNameValue( papszGlobalMetadata,  
     1088                                                   "mapProjection" ); 
     1089     
     1090    if( mapProjection ) 
     1091        iSDSIndex = SDnametoindex( hSD, mapProjection ); 
     1092 
     1093    if( iSDSIndex != FAIL ) 
     1094        iSDS = SDselect( hSD, iSDSIndex ); 
     1095        
     1096    if( iSDS != FAIL ) 
     1097    { 
     1098        char        szName[HDF4_SDS_MAXNAMELEN]; 
     1099        int32       iRank, iNumType, nAttrs; 
     1100        int32       aiDimSizes[MAX_VAR_DIMS]; 
     1101 
     1102        double adfGCTP[29]; 
     1103        int32 aiStart[MAX_NC_DIMS], aiEdges[MAX_NC_DIMS]; 
     1104 
     1105        aiStart[0] = 0; 
     1106        aiEdges[0] = 29; 
     1107 
     1108        if( SDgetinfo( iSDS, szName, &iRank, aiDimSizes, &iNumType,  
     1109                       &nAttrs) == 0  
     1110            && iNumType == DFNT_FLOAT64  
     1111            && iRank == 1 
     1112            && aiDimSizes[0] >= 29  
     1113            && SDreaddata( iSDS, aiStart, NULL, aiEdges, adfGCTP ) == 0  
     1114            && oSRS.importFromUSGS( (long) adfGCTP[1], (long) adfGCTP[2],  
     1115                                    adfGCTP+4,  
     1116                                    (long) adfGCTP[3] ) == OGRERR_NONE ) 
     1117        { 
     1118            CPLDebug( "HDF4Image", "GCTP Parms = %g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g", 
     1119                      adfGCTP[0],  
     1120                      adfGCTP[1],  
     1121                      adfGCTP[2],  
     1122                      adfGCTP[3],  
     1123                      adfGCTP[4],  
     1124                      adfGCTP[5],  
     1125                      adfGCTP[6],  
     1126                      adfGCTP[7],  
     1127                      adfGCTP[8],  
     1128                      adfGCTP[9],  
     1129                      adfGCTP[10],  
     1130                      adfGCTP[11],  
     1131                      adfGCTP[12],  
     1132                      adfGCTP[13],  
     1133                      adfGCTP[14],  
     1134                      adfGCTP[15],  
     1135                      adfGCTP[16],  
     1136                      adfGCTP[17],  
     1137                      adfGCTP[18],  
     1138                      adfGCTP[19],  
     1139                      adfGCTP[20],  
     1140                      adfGCTP[21],  
     1141                      adfGCTP[22],  
     1142                      adfGCTP[23],  
     1143                      adfGCTP[24],  
     1144                      adfGCTP[25],  
     1145                      adfGCTP[26],  
     1146                      adfGCTP[27],  
     1147                      adfGCTP[28] ); 
     1148             
     1149            CPLFree( pszProjection ); 
     1150            oSRS.exportToWkt( &pszProjection ); 
     1151            bGotGCTPProjection = TRUE; 
     1152        } 
     1153    } 
     1154 
     1155/* -------------------------------------------------------------------- */ 
     1156/*      If we derived a GCTP based projection, then we need to          */ 
     1157/*      transform the lat/long corners into this projection and use     */ 
     1158/*      them to establish the geotransform.                             */ 
     1159/* -------------------------------------------------------------------- */ 
     1160    if( bLLPossible && bGotGCTPProjection ) 
     1161    { 
     1162        double dfULX, dfULY, dfLRX, dfLRY; 
     1163        OGRSpatialReference oWGS84; 
     1164 
     1165        oWGS84.SetWellKnownGeogCS( "WGS84" ); 
     1166 
     1167        OGRCoordinateTransformation *poCT =  
     1168            OGRCreateCoordinateTransformation( &oWGS84, &oSRS ); 
     1169 
     1170        dfULX = adfXY[0*2+0]; 
     1171        dfULY = adfXY[0*2+1]; 
     1172         
     1173        dfLRX = adfXY[3*2+0]; 
     1174        dfLRY = adfXY[3*2+1]; 
     1175         
     1176        if( poCT->Transform( 1, &dfULX, &dfULY )  
     1177            && poCT->Transform( 1, &dfLRX, &dfLRY ) ) 
     1178        { 
     1179            bHasGeoTransform = TRUE; 
     1180            adfGeoTransform[0] = dfULX; 
     1181            adfGeoTransform[1] = (dfLRX - dfULX) / nRasterXSize; 
     1182            adfGeoTransform[2] = 0.0; 
     1183            adfGeoTransform[3] = dfULY; 
     1184            adfGeoTransform[4] = 0.0; 
     1185            adfGeoTransform[5] = (dfLRY - dfULY) / nRasterYSize; 
     1186        } 
     1187 
     1188        delete poCT; 
    10681189    } 
    10691190}