00001 #include "nlapsfrmt.h" 00002 #include <proj.h> 00003 #include <fstream> 00004 #include <iostream> 00005 #include "TraceLog.h" 00006 00007 const std::string NlapsFrmt::GCTPProjMap[] = { 00008 "latlong", 00009 "utm", 00010 "UNKNOWN STATE PLANE", 00011 "aea", 00012 "lcc", 00013 "merc", 00014 "ups", 00015 "poly", 00016 "eqdc", 00017 "tmerc", 00018 "stere", 00019 "laea", 00020 "aeqd", 00021 "gnom", 00022 "ortho", 00023 "nsper", 00024 "sinu", 00025 "UNKNOWN EQUIRECTANGULAR", 00026 "mill", 00027 "vandg", 00028 "omerc", 00029 "robin", 00030 "omerc", 00031 "robin", 00032 "lsat", 00033 "alsk", 00034 "goode", 00035 "hammer", 00036 "wag4", 00037 "wag7", 00038 "oea", 00039 "UNKNOWN INTEGERIZED SINUSOIDAL GRID" 00040 }; 00041 00042 const char* NlapsFrmt::GCTPParmMap[][15] = { 00043 { "", "", "", "", "", "", "", "", "", "", "", "", "", "", ""}, 00044 { "a=", "b=", "", "", "", "", "", "", "", "", "", "", "", "", ""}, 00045 { "", "", "", "", "", "", "", "", "", "", "", "", "", "", ""}, 00046 { "a=", "b=", "lat_1=", "lat_2=", "lon_0=", "lat_0=", "x_0=", "y_0=", "", "", "", "", "", "", ""}, 00047 { "a=", "b=", "lat_1=", "lat_2=", "lon_0=", "lat_0=", "x_0=", "y_0=", "", "", "", "", "", "", ""}, 00048 { "a=", "b=", "", "", "lon_0=", "k_0=", "x_0=", "y_0=", "", "", "", "", "", "", ""}, 00049 { "a=", "b=", "", "", "lon_0=", "k_0=", "x_0=", "y_0=", "", "", "", "", "", "", ""}, 00050 { "a=", "b=", "", "", "lon_0=", "lat_0=", "x_0=", "y_0=", "", "", "", "", "", "", ""}, 00051 { "a=", "b=", "lat_1=", "lat_2=", "lon_0=", "lat_0=", "x_0=", "y_0=", "", "", "", "", "", "", ""}, 00052 { "a=", "", "k=", "", "lon_0=", "lat_0=", "x_0=", "y_0=", "", "", "", "", "", "", ""}, 00053 { "a=1", "b=", "", "", "lon_0=", "lat_0=", "x_0=", "y_0=", "", "", "", "", "", "", ""}, 00054 { "a=1", "", "", "", "lon_0=", "lat_0=", "x_0=", "y_0=", "", "", "", "", "", "", ""}, 00055 { "a=1", "", "", "", "lon_0=", "lat_0=", "x_0=", "y_0=", "", "", "", "", "", "", ""}, 00056 { "a=1", "", "", "", "lon_0=", "lat_0=", "x_0=", "y_0=", "", "", "", "", "", "", ""}, 00057 { "a=1", "", "", "", "lon_0=", "lat_0=", "x_0=", "y_0=", "", "", "", "", "", "", ""}, 00058 { "a=1", "", "", "", "lon_0=", "lat_0=", "x_0=", "y_0=", "", "", "", "", "", "", ""}, 00059 { "a=1", "", "", "", "lon_0=", "", "x_0=", "y_0=", "", "", "", "", "", "", ""}, 00060 { "a=1", "", "", "", "lon_0=", "lat_ts=","x_0=", "y_0=", "", "", "", "", "", "", ""}, 00061 { "a=1", "", "", "", "lon_0=", "lat_0=", "x_0=", "y_0=", "", "", "", "", "", "", ""}, 00062 { "a=1", "", "", "", "lon_0=", "", "x_0=", "y_0=", "", "", "", "", "", "", ""}, 00063 { "a=", "b=", "k_0=", "", "", "lat_0=", "x_0=", "y_0=", "", "", "", "", "", "", ""}, 00064 { "a=1", "", "", "", "lon_0=", "", "x_0=", "y_0=", "", "", "", "", "", "", ""}, 00065 { "a=", "b=", "", "", "", "", "", "", "", "", "", "", "", "", ""}, 00066 { "a=", "b=", "", "", "", "", "x_0=", "y_0=", "", "", "", "", "", "", ""}, 00067 { "a=1", "", "", "", "", "", "", "", "", "", "", "", "", "", ""}, 00068 { "a=1", "", "", "", "lon_0=", "", "x_0=", "y_0=", "", "", "", "", "", "", ""}, 00069 { "a=1", "", "", "", "", "", "", "", "", "", "", "", "", "", ""}, 00070 { "a=1", "", "", "", "lon_0=", "", "x_0=", "y_0=", "", "", "", "", "", "", ""}, 00071 { "a=1", "", "", "", "lon_0=", "", "x_0=", "y_0=", "", "", "", "", "", "", ""}, 00072 { "a=1", "", "", "", "lon_0=", "", "x_0=", "y_0=", "", "", "", "", "", "", ""}, 00073 { "a=1", "", "", "", "", "", "", "", "", "", "", "", "", "", ""} 00074 }; 00075 00076 bool NlapsFrmt::open( std::string fname ) { 00077 mFName = fname; 00078 if ( !readFile() ) 00079 return false; 00080 00081 fillHeader(); 00082 fillBandInf(); 00083 fillImageInf(); 00084 return true; 00085 } 00086 00087 void* NlapsFrmt::getBand( int iband, 00088 int x1, 00089 int y1, 00090 int w1, 00091 int h1, 00092 int w2, 00093 int h2, 00094 ScaleType st ) const { 00095 std::string bfname; 00096 int ifs = -1; 00097 int pl; 00098 int scale; 00099 00100 int w = mBandSet[iband].width; 00101 // int h = mBandSet[iband].height; 00102 int dsz = dtypeSize( mBandSet[iband].type ); 00103 00104 unsigned char *data; 00105 unsigned char *line; 00106 00107 if ( w2 <= 0 || h2 <= 0 ) return NULL; 00108 00109 scale = (w1 / w2) + 1; 00110 00111 ifs = mFName.find_last_of( "\\" ); 00112 if ( ifs == -1 ) 00113 ifs = mFName.find_last_of( "/" ); 00114 if ( ifs == -1 ) return NULL; 00115 00116 bfname = mFName.substr( 0, ifs + 1 ); 00117 00118 std::ifstream in( ( bfname + mBFNameSet[iband] ).c_str(), 00119 std::istream::binary | 00120 std::istream::in ); 00121 if ( !in ) 00122 return NULL; 00123 00124 00125 data = new unsigned char[ ((w1 * h1) / scale) * dsz ]; 00126 00127 pl = ( y1 * w * dsz ) + ( x1 * dsz ); 00128 int sw = w1 / scale; 00129 int sh = h1 / scale; 00130 // int off = ( w * dsz ) * scale; 00131 line = new unsigned char[ w1 * dsz ]; 00132 00133 for ( int i=0; i < sh; i++ ) { 00134 00135 in.seekg( pl, std::istream::beg ); 00136 in.read( (char*)line, w1 * dsz ); 00137 00138 for ( int x=0; x < sw; x++ ) 00139 data[ (i * sw) + x ] = line[ x * scale ]; 00140 00141 pl += ( w * dsz ) * scale; 00142 } 00143 in.close(); 00144 00145 scaleData( (char**)&data, sw, sh, w2, h2, mBandSet[iband].type, st ); 00146 return(void*)data; 00147 00148 } 00149 00150 void NlapsFrmt::getLL( int iband, 00151 double &i, 00152 double &j ) { 00153 xyToLL( iband, i, j ); 00154 } 00155 00156 void NlapsFrmt::getXY( int iband, 00157 double &i, 00158 double &j ) { 00159 llToXY( iband, i, j ); 00160 } 00161 00162 bool NlapsFrmt::readFile( ) { 00163 // int fsz = 0; 00164 mText = ""; 00165 char ch; 00166 std::string fl; 00167 00168 std::ifstream in( mFName.c_str() ); 00169 if ( !in ) return false; 00170 00171 for ( int i=0; i < 20 && !in.eof(); i++ ) { 00172 in.get( ch ); 00173 fl += ch; 00174 } 00175 00176 if ( (int)fl.find( "NDF_REVISION" ) == -1 ) 00177 return false; 00178 00179 while ( !in.get( ch ).eof() ) 00180 mText += ch; 00181 00182 return true; 00183 } 00184 00185 std::vector<std::string> 00186 NlapsFrmt::getValue( std::string id ) { 00187 std::vector<std::string> ret; 00188 std::string field, tid, tval; 00189 int line = 0; 00190 bool foundOne = false; 00191 00192 do { 00193 if ( !getToken( mText, "\n", line++, field ) ) 00194 break; 00195 00196 field = trimStr( field ); 00197 getToken( field, "=", 0, tid ); 00198 tid = trimStr( tid ); 00199 00200 if ( tid == id ) { 00201 getToken( field, "=", 1, tval ); 00202 tval = trimStr( tval ); 00203 foundOne = true; 00204 } 00205 } while ( !foundOne ); 00206 00207 if ( foundOne ) { 00208 int i=0; 00209 std::string val; 00210 while ( getToken( tval, ",", i++, val ) ) 00211 ret.push_back( val ); 00212 } 00213 00214 return ret; 00215 } 00216 00217 void NlapsFrmt::fillHeader( ) { 00218 std::vector<std::string> valset; 00219 00220 // DATA_FILE_INTERLEAVE 00221 valset = getValue( "DATA_FILE_INTERLEAVING" ); 00222 if ( valset[0] == "BSQ" ) 00223 mHead.bsq = true; 00224 00225 valset = getValue( "USGS_PROJECTION_NUMBER" ); 00226 mHead.gctpID = atoi( valset[0].c_str() ); 00227 00228 valset = getValue( "USGS_MAP_ZONE" ); 00229 if ( valset.size() > 0 ) 00230 mHead.utmz = atoi( valset[0].c_str() ); 00231 00232 valset = getValue( "HORIZONTAL_DATUM" ); 00233 mHead.datum = valset[0]; 00234 00235 valset = getValue( "UPPER_LEFT_CORNER" ); 00236 mHead.ul_lon = parseDeg( valset[0] ); 00237 mHead.ul_lat = parseDeg( valset[1] ); 00238 } 00239 00240 void NlapsFrmt::fillBandInf( ) { 00241 BandInf binf; 00242 std::vector<std::string> val; 00243 char id[255]; 00244 00245 binf.proj = GCTPProjMap[mHead.gctpID]; 00246 binf.UTMZ = mHead.utmz; 00247 binf.rot = 0.0; 00248 binf.ulx = mHead.ul_lon; 00249 binf.uly = mHead.ul_lat; 00250 00251 val = getValue( "PIXEL_FORMAT" ); 00252 if ( val[0] == "BYTE" ) 00253 binf.type = UINT8; 00254 else if ( val[0] == "2BYTEINT" ) 00255 binf.type = INT16; 00256 else if ( val[0] == "4BYTEINT" ) 00257 binf.type = INT32; 00258 else if ( val[0] == "REAL" ) 00259 binf.type = INT32; 00260 else if ( val[0] == "DOUBLE" ) 00261 binf.type = FLOAT32; 00262 else 00263 binf.type = NOT_DEFINED; 00264 00265 val = getValue( "USGS_PROJECTION_PARAMETERS" ); 00266 binf.parms = getProjParms( val, mHead.gctpID ); 00267 00268 00269 val = getValue( "PIXEL_SPACING" ); 00270 binf.res = atoi( val[0].c_str() ); 00271 00272 val = getValue( "LINES_PER_DATA_FILE" ); 00273 binf.height = atoi( val[0].c_str() ); 00274 00275 val = getValue( "PIXELS_PER_LINE" ); 00276 binf.width = atoi( val[0].c_str() ); 00277 00278 binf.ndims = 2; 00279 00280 int cnt = 1; 00281 while ( true ) { 00282 BandInf cbinf = binf; 00283 sprintf( id, "BAND%d_NAME", cnt ); 00284 val = getValue( id ); 00285 if ( val.size() < 1 ) break; // THIS SHOULD STOP THE LOOP 00286 cbinf.name = val[0]; 00287 00288 sprintf( id, "BAND%d_FILENAME", cnt ); 00289 val = getValue( id ); 00290 mBFNameSet.push_back( val[0] ); 00291 00292 mBandSet.push_back( cbinf ); 00293 00294 cnt++; 00295 } 00296 00297 getULCorners(); 00298 } 00299 00300 void NlapsFrmt::fillImageInf( ) { 00301 double x, y; 00302 std::vector<std::string> val; 00303 00304 x = mBandSet[0].width / 2; 00305 y = mBandSet[0].height / 2; 00306 xyToLL( 0, x, y ); 00307 00308 mImageInf.mCenterLat = y; 00309 mImageInf.mCenterLon = x; 00310 00311 val = getValue( "MAP_PROJECTION_NAME" ); 00312 if ( (int)val.size() > 0 ) 00313 mImageInf.mProjection = val[0]; 00314 00315 val = getValue( "ACQUISITION_DATE/TIME" ); 00316 if ( (int)val.size() > 0 && mImageInf.mAquiDateTime.size() == 0) 00317 mImageInf.mAquiDateTime = formatDateTime( val[0] ); 00318 } 00319 00320 std::string NlapsFrmt::formatDateTime( std::string str ) { 00321 // 1999-11-04T17:49:07Z 00322 00323 std::string ret; 00324 00325 ret = str.substr( 0, 4 ) + 00326 str.substr( 5, 2 ) + 00327 str.substr( 8, 2 ) + " " + 00328 str.substr( 11, 8 ) + "-0000"; 00329 00330 return ret; 00331 } 00332 00333 double NlapsFrmt::parseDeg( std::string sdeg ) { 00334 double deg = 0.0; 00335 double min = 0.0; 00336 double sec = 0.0; 00337 char dir; 00338 double val; 00339 00340 deg = atoi( sdeg.substr( 0, 3 ).c_str() ); 00341 min = atoi( sdeg.substr( 3, 2 ).c_str() ); 00342 sec = atof( sdeg.substr( 5, 7 ).c_str() ); 00343 dir = sdeg[12]; 00344 00345 val = deg + (min/60.0) + (sec/60.0/60.0); 00346 if ( dir == 'S' || dir == 'W' ) 00347 val = 0 - val; 00348 00349 return val; 00350 } 00351 00352 bool NlapsFrmt::getToken( std::string str, 00353 const std::string delim, 00354 int n, 00355 std::string &tok ) { 00356 int l_ind=0, r_ind=0; 00357 int count=-1; 00358 while ( ( r_ind = str.find(delim, l_ind) ) != -1 && ++count < n ) 00359 l_ind = r_ind + 1; 00360 00361 if ( count == n || (n == count+1 && r_ind == -1) ) { 00362 if ( r_ind == -1 ) 00363 r_ind = str.length(); 00364 tok = str.substr(l_ind, r_ind - l_ind); 00365 if ( tok == "" ) 00366 return false; 00367 return true; 00368 } 00369 00370 tok = ""; 00371 return false; 00372 } 00373 00374 std::string NlapsFrmt::trimStr( std::string str ) { 00375 int a=0; 00376 int strsz = str.length(); 00377 while ( a < strsz ) { 00378 char ch = str[a]; 00379 if ( ch != '\r' && ch != '\n' && ch != ' ' && ch != '\t' && ch != ';' ) 00380 break; 00381 a++; 00382 } 00383 00384 if ( a >= strsz ) 00385 return ""; 00386 00387 int b = strsz; 00388 while ( b > a ) { 00389 char ch = str[b-1]; 00390 if ( ch != '\r' && ch != '\n' && ch != ' ' && ch != '\t' && ch != ';' ) 00391 break; 00392 b--; 00393 } 00394 00395 if ( b <= a ) 00396 return ""; 00397 00398 return str.substr( a, b-a ); 00399 } 00400 00401 std::vector<std::string> 00402 NlapsFrmt::getProjParms( std::vector<std::string> &raw, int gctpid ) { 00403 std::vector<std::string> ret; 00404 char buf[100]; 00405 00406 if ( (int)raw.size() != 15 ) 00407 return ret; 00408 00409 for ( int i=0; i < 15; i++ ) { 00410 std::string pjparm = GCTPParmMap[gctpid][i]; 00411 if ( pjparm != "" ) { 00412 if ( pjparm[pjparm.length()-1] == '=' ) { 00413 double val = atof( raw[i].c_str() ); 00414 if ( (int)pjparm.find( "lat" ) != -1 || (int)pjparm.find( "lon" ) != -1 ) { 00415 int deg = (int)(val / 1000000); 00416 double min = (val - (deg*1000000)) / 1000; 00417 double sec = val - (deg*1000000) - (min*1000); 00418 00419 min /= 60; 00420 sec /= 360.0; 00421 00422 val = deg + min + sec; 00423 } 00424 00425 sprintf( buf, "%s%.5f", pjparm.c_str(), val ); 00426 ret.push_back( std::string( buf ) ); 00427 } else { 00428 ret.push_back( pjparm ); 00429 } 00430 } 00431 } 00432 00433 return ret; 00434 } 00435 00436 00437 00438 00439 00440
Home |
Search |
Disclaimers & Privacy |
Contact Us GLIMSView Maintainer: dsoltesz@usgs.gov |