00001 00002 00003 00004 00005 00006 00007 00008 00009 00010 00011 00012 00013 00014 00015 00016 00017 00018 00019 00020 00021 00022 00023 00024 00025 00026 00027 00028 00029 00030 00031 00032 00032 00032 00032 00033 00034 00035 00036 00037 00038 00039 00040 00041 00042 00043 00044 00045 00046 00047 00048 00049 00050 00051 00052 00053 00054 00055
00056
00057 #include "gdal_priv.h"
00058 #include "cpl_conv.h"
00059 #include "cpl_string.h"
00060 #include "ogr_spatialref.h"
00061
00062 typedef struct { double u, v; } UV;
00063
00064 #define PJ void
00065
00066 static PJ *(*pfn_pj_init)(int, char**) = NULL;
00067 static UV (*pfn_pj_fwd)(UV, PJ *) = NULL;
00068 static UV (*pfn_pj_inv)(UV, PJ *) = NULL;
00069 static void (*pfn_pj_free)(PJ *) = NULL;
00070
00071 #define RAD_TO_DEG 57.29577951308232
00072 #define DEG_TO_RAD .0174532925199432958
00073
00074 #ifdef WIN32
00075 # define LIBNAME "proj.dll"
00076 #else
00077 # define LIBNAME "libproj.so"
00078 #endif
00079
00080
00081
00082
00083
00084 static int LoadProjLibrary()
00085
00086 {
00087 static int bTriedToLoad = FALSE;
00088
00089 if( bTriedToLoad )
00090 return( pfn_pj_init != NULL );
00091
00092 bTriedToLoad = TRUE;
00093
00094 CPLPushErrorHandler( CPLQuietErrorHandler );
00095 pfn_pj_init = (PJ *(*)(int, char**)) CPLGetSymbol( LIBNAME,
00096 "pj_init" );
00097 CPLPopErrorHandler();
00098
00099 if( pfn_pj_init == NULL )
00100 return( FALSE );
00101
00102 pfn_pj_fwd = (UV (*)(UV,PJ*)) CPLGetSymbol( LIBNAME, "pj_fwd" );
00103 pfn_pj_inv = (UV (*)(UV,PJ*)) CPLGetSymbol( LIBNAME, "pj_inv" );
00104 pfn_pj_free = (void (*)(PJ*)) CPLGetSymbol( LIBNAME, "pj_free" );
00105
00106 return( TRUE );
00107 }
00108
00109
00110
00111
00112
00113 GDALProjDef::GDALProjDef( const char * pszProjectionIn )
00114
00115 {
00116 if( pszProjectionIn == NULL )
00117 pszProjection = CPLStrdup( "" );
00118 else
00119 pszProjection = CPLStrdup( pszProjectionIn );
00120
00121 psPJ = NULL;
00122
00123 SetProjectionString( pszProjectionIn );
00124 }
00125
00126
00127
00128
00129
00130 GDALProjDefH GDALCreateProjDef( const char * pszProjection )
00131
00132 {
00133 return( (GDALProjDefH) (new GDALProjDef( pszProjection )) );
00134 }
00135
00136
00137
00138
00139
00140 GDALProjDef::~GDALProjDef()
00141
00142 {
00143 CPLFree( pszProjection );
00144 if( psPJ != NULL )
00145 pfn_pj_free( psPJ );
00146 }
00147
00148
00149
00150
00151
00152 void GDALDestroyProjDef( GDALProjDefH hProjDef )
00153
00154 {
00155 delete (GDALProjDef *) hProjDef;
00156 }
00157
00158
00159
00160
00161
00162 CPLErr GDALProjDef::SetProjectionString( const char * pszProjectionIn )
00163
00164 {
00165 char **args;
00166
00167 if( psPJ != NULL && pfn_pj_free != NULL )
00168 {
00169 pfn_pj_free( psPJ );
00170 psPJ = NULL;
00171 }
00172
00173 if( pszProjection != NULL )
00174 CPLFree( pszProjection );
00175
00176 pszProjection = CPLStrdup( pszProjectionIn );
00177
00178 if( !LoadProjLibrary() )
00179 {
00180 return CE_Failure;
00181 }
00182
00183
00184
00185
00186
00187 char *pszProj4Projection = NULL;
00188
00189 if( EQUALN(pszProjection,"PROJCS",6) || EQUALN(pszProjection,"GEOGCS",6) )
00190 {
00191 OGRSpatialReference oSRS;
00192 char *pszProjRef = pszProjection;
00193
00194 if( oSRS.importFromWkt( &pszProjRef ) != OGRERR_NONE
00195 || oSRS.exportToProj4( &pszProj4Projection ) != OGRERR_NONE )
00196 return CE_Failure;
00197 }
00198 else
00199 {
00200 pszProj4Projection = CPLStrdup( pszProjection );
00201 }
00202
00203
00204
00205
00206 args = CSLTokenizeStringComplex( pszProj4Projection, " +", TRUE, FALSE );
00207
00208 psPJ = pfn_pj_init( CSLCount(args), args );
00209
00210
00211
00212
00213 CSLDestroy( args );
00214 CPLFree( pszProj4Projection );
00215
00216 if( psPJ == NULL )
00217 return CE_Failure;
00218 else
00219 return CE_None;
00220 }
00221
00222
00223
00224
00225
00226 CPLErr GDALProjDef::ToLongLat( double * padfX, double * padfY )
00227
00228 {
00229 UV uv;
00230
00231 CPLAssert( padfX != NULL && padfY != NULL );
00232
00233 if( strstr(pszProjection,"+proj=longlat") != NULL
00234 || strstr(pszProjection,"+proj=latlong") != NULL
00235 || EQUALN(pszProjection,"GEOGCS",6) )
00236 return CE_None;
00237
00238 if( psPJ == NULL )
00239 return CE_Failure;
00240
00241 uv.u = *padfX;
00242 uv.v = *padfY;
00243
00244 uv = pfn_pj_inv( uv, psPJ );
00245
00246 *padfX = uv.u * RAD_TO_DEG;
00247 *padfY = uv.v * RAD_TO_DEG;
00248
00249 return( CE_None );
00250 }
00251
00252
00253
00254
00255
00256 CPLErr GDALReprojectToLongLat( GDALProjDefH hProjDef,
00257 double * padfX, double * padfY )
00258
00259 {
00260 return( ((GDALProjDef *) hProjDef)->ToLongLat(padfX, padfY) );
00261 }
00262
00263
00264
00265
00266
00267 CPLErr GDALProjDef::FromLongLat( double * padfX, double * padfY )
00268
00269 {
00270 UV uv;
00271
00272 CPLAssert( padfX != NULL && padfY != NULL );
00273
00274 if( strstr(pszProjection,"+proj=longlat") != NULL
00275 || strstr(pszProjection,"+proj=latlong") != NULL
00276 || EQUALN(pszProjection,"GEOGCS",6) )
00277 return CE_None;
00278
00279 if( psPJ == NULL )
00280 return CE_Failure;
00281
00282 uv.u = *padfX * DEG_TO_RAD;
00283 uv.v = *padfY * DEG_TO_RAD;
00284
00285 uv = pfn_pj_fwd( uv, psPJ );
00286
00287 *padfX = uv.u;
00288 *padfY = uv.v;
00289
00290 return( CE_None );
00291 }
00292
00293
00294
00295
00296
00297 CPLErr GDALReprojectFromLongLat( GDALProjDefH hProjDef,
00298 double * padfX, double * padfY )
00299
00300 {
00301 return( ((GDALProjDef *) hProjDef)->FromLongLat(padfX, padfY) );
00302 }
00303
00304
00305
00306
00307
00308
00309
00310
00311 const char *GDALDecToDMS( double dfAngle, const char * pszAxis,
00312 int nPrecision )
00313
00314 {
00315 int nDegrees, nMinutes;
00316 double dfSeconds;
00317 char szFormat[30];
00318 static char szBuffer[50];
00319 const char *pszHemisphere;
00320
00321
00322 nDegrees = (int) ABS(dfAngle);
00323 nMinutes = (int) ((ABS(dfAngle) - nDegrees) * 60);
00324 dfSeconds = (ABS(dfAngle) * 3600 - nDegrees*3600 - nMinutes*60);
00325
00326 if( EQUAL(pszAxis,"Long") && dfAngle < 0.0 )
00327 pszHemisphere = "W";
00328 else if( EQUAL(pszAxis,"Long") )
00329 pszHemisphere = "E";
00330 else if( dfAngle < 0.0 )
00331 pszHemisphere = "S";
00332 else
00333 pszHemisphere = "N";
00334
00335 sprintf( szFormat, "%%3dd%%2d\'%%.%df\"%s", nPrecision, pszHemisphere );
00336 sprintf( szBuffer, szFormat, nDegrees, nMinutes, dfSeconds );
00337
00338 return( szBuffer );
00339 }
00340