C++ OGRCoordinateTransformation::Transform方法代码示例

/** Fetch the longitude/latitude bounds of an image
 * @param poDS a pointer to a valid GDAL Dataset
 * @param boundsLonLat a pointer to a double[4] n, e, s, w order
 * @return true on valid population of boundsLonLat in n, e, s, w order
 */
bool GDALGetBounds( GDALDataset *poDS, double *boundsLonLat )
{
    char* pszPrj;
    double adfGeoTransform[6];
    int xSize, ySize;
    double nLat, eLon, sLat, wLon;

    OGRSpatialReference oSourceSRS, oTargetSRS;
    OGRCoordinateTransformation *poCT;

    if( poDS == NULL )
	return false;

    xSize = poDS->GetRasterXSize ();
    ySize = poDS->GetRasterYSize();

    poDS->GetGeoTransform( adfGeoTransform );

    if( poDS->GetProjectionRef() == NULL )
	return false;
    else
	pszPrj = (char*)poDS->GetProjectionRef();

    oSourceSRS.importFromWkt( &pszPrj );
    oTargetSRS.SetWellKnownGeogCS( "WGS84" );

    poCT = OGRCreateCoordinateTransformation( &oSourceSRS, &oTargetSRS );
    if( poCT == NULL )
	return false;

    nLat = adfGeoTransform[3] + adfGeoTransform[4] * 0 + adfGeoTransform[5] * 0;
    eLon = adfGeoTransform[0] + adfGeoTransform[1] * xSize + adfGeoTransform[2] * 0;
    sLat = adfGeoTransform[3] + adfGeoTransform[4] * 0 + adfGeoTransform[5] * ySize;
    wLon = adfGeoTransform[0] + adfGeoTransform[1] * 0 + adfGeoTransform[2] * 0;

    if( !poCT->Transform( 1, &eLon, &nLat ) ) {
	OGRCoordinateTransformation::DestroyCT( poCT );
	return false;
    }

    boundsLonLat[0] = nLat;
    boundsLonLat[1] = eLon;

    if( !poCT->Transform( 1, &wLon, &sLat ) ) {
	OGRCoordinateTransformation::DestroyCT( poCT );
	return false;
    }

    boundsLonLat[2] = sLat;
    boundsLonLat[3] = wLon;

    OGRCoordinateTransformation::DestroyCT( poCT );
    return true;
}

CPLErr HDF5ImageDataset::CreateODIMH5Projection()
{
    const char* const pszProj4String = GetMetadataItem("where_projdef");
    const char* const pszLL_lon = GetMetadataItem("where_LL_lon");
    const char* const pszLL_lat = GetMetadataItem("where_LL_lat");
    const char* const pszUR_lon = GetMetadataItem("where_UR_lon");
    const char* const pszUR_lat = GetMetadataItem("where_UR_lat");
    if( pszProj4String == NULL ||
            pszLL_lon == NULL || pszLL_lat == NULL ||
            pszUR_lon == NULL || pszUR_lat == NULL )
        return CE_Failure;

    if( oSRS.importFromProj4( pszProj4String ) != OGRERR_NONE )
        return CE_Failure;

    OGRSpatialReference oSRSWGS84;
    oSRSWGS84.SetWellKnownGeogCS( "WGS84" );

    OGRCoordinateTransformation* poCT =
        OGRCreateCoordinateTransformation( &oSRSWGS84, &oSRS );
    if( poCT == NULL )
        return CE_Failure;

    /* Reproject corners from long,lat WGS84 to the target SRS */
    double dfLLX = CPLAtof(pszLL_lon);
    double dfLLY = CPLAtof(pszLL_lat);
    double dfURX = CPLAtof(pszUR_lon);
    double dfURY = CPLAtof(pszUR_lat);
    if( !poCT->Transform(1, &dfLLX, &dfLLY) ||
            !poCT->Transform(1, &dfURX, &dfURY) )
    {
        delete poCT;
        return CE_Failure;
    }
    delete poCT;

    /* Compute the geotransform now */
    const double dfPixelX = (dfURX - dfLLX) / nRasterXSize;
    const double dfPixelY = (dfURY - dfLLY) / nRasterYSize;

    bHasGeoTransform = true;
    adfGeoTransform[0] = dfLLX;
    adfGeoTransform[1] = dfPixelX;
    adfGeoTransform[2] = 0;
    adfGeoTransform[3] = dfURY;
    adfGeoTransform[4] = 0;
    adfGeoTransform[5] = -dfPixelY;

    CPLFree( pszProjection );
    oSRS.exportToWkt( &pszProjection );

    return CE_None;
}

bool GDALPointToLatLon( double &x, double &y, GDALDataset *poSrcDS,
				const char *datum )
{
    char* pszPrj = NULL;

    OGRSpatialReference oSourceSRS, oTargetSRS;
    OGRCoordinateTransformation *poCT;

    if( poSrcDS == NULL )
	return false;

    if( poSrcDS->GetProjectionRef() == NULL )
	return false;
    else
	pszPrj = (char*)poSrcDS->GetProjectionRef();

    oSourceSRS.importFromWkt( &pszPrj );
    oTargetSRS.SetWellKnownGeogCS( datum );

    poCT = OGRCreateCoordinateTransformation( &oSourceSRS, &oTargetSRS );

    if( poCT == NULL )
	return false;

    if( !poCT->Transform( 1, &x, &y ) ) {
	OGRCoordinateTransformation::DestroyCT( poCT );
	return false;
    }
    OGRCoordinateTransformation::DestroyCT( poCT );
    return true;
}

std::pair<double,double> Terrain::getGeoTransformFromDEMs(Terrain *large, Terrain *small)
{
    QString proj = large->dataset->GetProjectionRef();

    OGRSpatialReference sr1;
    auto t = proj.toLatin1();
    char *test = t.data();
    sr1.importFromWkt(&test);

    proj = small->dataset->GetProjectionRef();

    t = proj.toLatin1();
    OGRSpatialReference sr2;
    test = t.data();
    sr2.importFromWkt(&test);

    OGRCoordinateTransformation* poTransform = OGRCreateCoordinateTransformation( &sr2, &sr1 );
    large->geot.resize(6);
    small->geot.resize(6);
    large->dataset->GetGeoTransform(large->geot.data());
    small->dataset->GetGeoTransform(small->geot.data());
    //
    double x = small->geot[0];
    double y = small->geot[3];
    // DCEWsqrext.tif upperleft hand corner is convert to tl2p5_dem.tif coordinate system.
    poTransform->Transform (1, &x, &y);

    return std::pair<double,double>(x,y);
}

/**
*@brief 数组对数组的转换
*@param fromWkt [in] 原点的wkt字符串
*@param toWkt [in] 目标点的wkt字符串
*@param x [in] 经度数组
*@param y [in] 纬度数组
*@param nCount [in] 数组大小
*@return 转换后的经纬度坐标
*/
vector<double*> GdalProjection::ArrayToArray(const char *fromWkt,const char *toWkt,double *x,double *y,int nCount)
{
    vector<double*> ve;
    OGRSpatialReference oSourceSRS,oTargetSRS;
    OGRCoordinateTransformation *poct;

    oSourceSRS.importFromWkt((char **) &fromWkt);
    oTargetSRS.importFromWkt((char **) &toWkt);

    poct = OGRCreateCoordinateTransformation(&oSourceSRS,&oTargetSRS);

    if (poct == NULL || !poct->Transform(nCount,x,y))
    {
        x = NULL;
        y = NULL;
        ve.push_back(x);
        ve.push_back(y);
    }
    else
    {
        ve.push_back(x);
        ve.push_back(y);
    }
    return ve;
}

/** Fetch the center of a domain.
 * Fetch the center of a domain from any valid GDAL dataset
 * @param poDS a pointer to a valid GDAL Dataset
 * @param centerLonLat a pointer to a double size of double * 2
 * @return true on valid population of the double*
 */
bool GDALGetCenter( GDALDataset *poDS, double *centerLonLat )
{
    char* pszPrj;
    double adfGeoTransform[6];
    int xSize, ySize;
    double xCenter, yCenter;
    double lon, lat;

    OGRSpatialReference oSourceSRS, oTargetSRS;
    OGRCoordinateTransformation *poCT;

    if( poDS == NULL )
	return false;

    xSize = poDS->GetRasterXSize( );
    ySize = poDS->GetRasterYSize( );

    if( poDS->GetGeoTransform( adfGeoTransform ) != CE_None )
	return false;

    if( poDS->GetProjectionRef(  ) == NULL )
	return false;
    else
	pszPrj = (char*)poDS->GetProjectionRef();

    oSourceSRS.importFromWkt( &pszPrj );
    oTargetSRS.SetWellKnownGeogCS( "WGS84" );

    poCT = OGRCreateCoordinateTransformation( &oSourceSRS, &oTargetSRS );
    if( poCT == NULL )
	return false;

    xCenter = xSize / 2;
    yCenter = ySize / 2;

    lon = adfGeoTransform[0] + adfGeoTransform[1] * xCenter
	+ adfGeoTransform[2] * yCenter;

    lat = adfGeoTransform[3] + adfGeoTransform[4] * xCenter
	+ adfGeoTransform[5] * yCenter;

    if( !poCT->Transform( 1, &lon, &lat ) ) {
	OGRCoordinateTransformation::DestroyCT( poCT );
	return false;
    }

    centerLonLat[0] = lon;
    centerLonLat[1] = lat;

    OGRCoordinateTransformation::DestroyCT( poCT );
    return true;
}

void Convert_Geographic_To_UTM( const int&          number_coordinates,
                                const double*       latitudes_degrees,
                                const double*       longitudes_degrees,
                                const std::string&  input_datum,
                                const int&          grid_zone,
                                const bool&         is_northern,
                                double*             eastings_meters,
                                double*             northings_meters,
                                const std::string&  output_datum )
{
    // Create the Spatial Reference Objects
    OGRSpatialReference sourceSRS, targetSRS;

    sourceSRS.SetWellKnownGeogCS(  input_datum.c_str() );
    targetSRS.SetWellKnownGeogCS( output_datum.c_str() );

    // Configure the Projected Coordinate Components
    targetSRS.SetUTM( grid_zone,
                      is_northern );

    // Build the Transform Engine
    OGRCoordinateTransformation* transform;
    transform = OGRCreateCoordinateTransformation( &sourceSRS,
                                                   &targetSRS );
    double* elevations_meters = NULL;

    memcpy( eastings_meters,  longitudes_degrees, sizeof(double)*number_coordinates );
    memcpy( northings_meters, latitudes_degrees,  sizeof(double)*number_coordinates );

    if( !transform->Transform( number_coordinates, 
                               eastings_meters,
                               northings_meters,
                               NULL ))
    {
        throw std::runtime_error("Transformation Failed.");
    }

    
    // Destroy the Transform
    OCTDestroyCoordinateTransformation( transform );
}

bool OGRPointToLatLon(double &x, double &y, OGRDataSourceH hDS,
                      const char *datum) {
  char *pszPrj = NULL;

  OGRSpatialReference *poSrcSRS;
  OGRSpatialReference oSourceSRS, oTargetSRS;
  OGRCoordinateTransformation *poCT;

  if (hDS == NULL) {
    return false;
  }

  OGRLayer *poLayer;

  poLayer = (OGRLayer *)OGR_DS_GetLayer(hDS, 0);
  poLayer->ResetReading();

  poSrcSRS = poLayer->GetSpatialRef();
  if (poSrcSRS == NULL) {
    return false;
  }

  oTargetSRS.SetWellKnownGeogCS(datum);

  poCT = OGRCreateCoordinateTransformation(poSrcSRS, &oTargetSRS);

  if (poCT == NULL) {
    return false;
  }

  if (!poCT->Transform(1, &x, &y)) {
    OGRCoordinateTransformation::DestroyCT(poCT);
    return false;
  }
  OGRCoordinateTransformation::DestroyCT(poCT);
  return true;
}

OGREnvelope wxGISRasterRGBRenderer::TransformEnvelope(OGREnvelope* pEnvelope, OGRSpatialReference* pSrsSpatialReference, OGRSpatialReference* pDstSpatialReference)
{
	//get new envelope - it may rotate
	OGRCoordinateTransformation *poCT = OGRCreateCoordinateTransformation( pSrsSpatialReference, pDstSpatialReference);
	double pointsx[4];
	double pointsy[4];
	pointsx[0] = pEnvelope->MaxX;
	pointsy[0] = pEnvelope->MaxY;
	pointsx[1] = pEnvelope->MinX;
	pointsy[1] = pEnvelope->MinY;
	pointsx[2] = pEnvelope->MaxX;
	pointsy[2] = pEnvelope->MinY;
	pointsx[3] = pEnvelope->MinX;
	pointsy[3] = pEnvelope->MaxY;
	//get real envelope
	poCT->Transform(4, pointsx, pointsy);
	OCTDestroyCoordinateTransformation(poCT);
	OGREnvelope out;
	out.MinX = MIN(pointsx[0], MIN(pointsx[1], MIN(pointsx[2], pointsx[3])));
	out.MaxX = MAX(pointsx[0], MAX(pointsx[1], MAX(pointsx[2], pointsx[3])));
	out.MinY = MIN(pointsy[0], MIN(pointsy[1], MIN(pointsy[2], pointsy[3])));
	out.MaxY = MAX(pointsy[0], MAX(pointsy[1], MAX(pointsy[2], pointsy[3])));
	return out;
}

void GRIBDataset::SetGribMetaData(grib_MetaData* meta)
{
    nRasterXSize = meta->gds.Nx;
    nRasterYSize = meta->gds.Ny;

/* -------------------------------------------------------------------- */
/*      Image projection.                                               */
/* -------------------------------------------------------------------- */
    OGRSpatialReference oSRS;

    switch(meta->gds.projType)
    {
      case GS3_LATLON:
      case GS3_GAUSSIAN_LATLON:
          // No projection, only latlon system (geographic)
          break;
      case GS3_MERCATOR:
        oSRS.SetMercator(meta->gds.meshLat, meta->gds.orientLon,
                         1.0, 0.0, 0.0);
        break;
      case GS3_POLAR:
        oSRS.SetPS(meta->gds.meshLat, meta->gds.orientLon,
                   meta->gds.scaleLat1,
                   0.0, 0.0);
        break;
      case GS3_LAMBERT:
        oSRS.SetLCC(meta->gds.scaleLat1, meta->gds.scaleLat2,
                    0.0, meta->gds.orientLon,
                    0.0, 0.0); // set projection
        break;
			

      case GS3_ORTHOGRAPHIC:

        //oSRS.SetOrthographic(0.0, meta->gds.orientLon,
        //											meta->gds.lon2, meta->gds.lat2);
        //oSRS.SetGEOS(meta->gds.orientLon, meta->gds.stretchFactor, meta->gds.lon2, meta->gds.lat2);
        oSRS.SetGEOS(  0, 35785831, 0, 0 ); // hardcoded for now, I don't know yet how to parse the meta->gds section
        break;
      case GS3_EQUATOR_EQUIDIST:
        break;
      case GS3_AZIMUTH_RANGE:
        break;
    }

/* -------------------------------------------------------------------- */
/*      Earth model                                                     */
/* -------------------------------------------------------------------- */
    double a = meta->gds.majEarth * 1000.0; // in meters
    double b = meta->gds.minEarth * 1000.0;
    if( a == 0 && b == 0 )
    {
        a = 6377563.396;
        b = 6356256.910;
    }

    if (meta->gds.f_sphere)
    {
        oSRS.SetGeogCS( "Coordinate System imported from GRIB file",
                        NULL,
                        "Sphere",
                        a, 0.0 );
    }
    else
    {
        double fInv = a/(a-b);
        oSRS.SetGeogCS( "Coordinate System imported from GRIB file",
                        NULL,
                        "Spheroid imported from GRIB file",
                        a, fInv );
    }

    OGRSpatialReference oLL; // construct the "geographic" part of oSRS
    oLL.CopyGeogCSFrom( &oSRS );

    double rMinX;
    double rMaxY;
    double rPixelSizeX;
    double rPixelSizeY;
    if (meta->gds.projType == GS3_ORTHOGRAPHIC)
    {
        //rMinX = -meta->gds.Dx * (meta->gds.Nx / 2); // This is what should work, but it doesn't .. Dx seems to have an inverse relation with pixel size
        //rMaxY = meta->gds.Dy * (meta->gds.Ny / 2);
        const double geosExtentInMeters = 11137496.552; // hardcoded for now, assumption: GEOS projection, full disc (like MSG)
        rMinX = -(geosExtentInMeters / 2);
        rMaxY = geosExtentInMeters / 2;
        rPixelSizeX = geosExtentInMeters / meta->gds.Nx;
        rPixelSizeY = geosExtentInMeters / meta->gds.Ny;
    }
    else if( oSRS.IsProjected() )
    {
        rMinX = meta->gds.lon1; // longitude in degrees, to be transformed to meters (or degrees in case of latlon)
        rMaxY = meta->gds.lat1; // latitude in degrees, to be transformed to meters 
        OGRCoordinateTransformation *poTransformLLtoSRS = OGRCreateCoordinateTransformation( &(oLL), &(oSRS) );
        if ((poTransformLLtoSRS != NULL) && poTransformLLtoSRS->Transform( 1, &rMinX, &rMaxY )) // transform it to meters
        {
            if (meta->gds.scan == GRIB2BIT_2) // Y is minY, GDAL wants maxY
                rMaxY += (meta->gds.Ny - 1) * meta->gds.Dy; // -1 because we GDAL needs the coordinates of the centre of the pixel
            rPixelSizeX = meta->gds.Dx;
            rPixelSizeY = meta->gds.Dy;
//.........这里部分代码省略.........

//.........这里部分代码省略.........

        else if( EQUALN(pszPR,"PR=TRANSVERSE MERCATOR", 22)
                 && osPP.size() > 0 )
        {
            
            osUnderlyingSRS.Printf( 
                "PROJCS[\"unnamed\",%s,PROJECTION[\"Transverse_Mercator\"],PARAMETER[\"latitude_of_origin\",0],PARAMETER[\"central_meridian\",%s],PARAMETER[\"scale_factor\",1],PARAMETER[\"false_easting\",0],PARAMETER[\"false_northing\",0]]",
                pszGEOGCS, osPP.c_str() );
        }

        else if( EQUALN(pszPR,"PR=UNIVERSAL TRANSVERSE MERCATOR", 32)
                 && osPP.size() > 0 )
        {
            // This is not *really* UTM unless the central meridian 
            // matches a zone which it does not in some (most?) maps. 
            osUnderlyingSRS.Printf( 
                "PROJCS[\"unnamed\",%s,PROJECTION[\"Transverse_Mercator\"],PARAMETER[\"latitude_of_origin\",0],PARAMETER[\"central_meridian\",%s],PARAMETER[\"scale_factor\",0.9996],PARAMETER[\"false_easting\",500000],PARAMETER[\"false_northing\",0]]", 
                pszGEOGCS, osPP.c_str() );
        }

        else if( EQUALN(pszPR,"PR=POLYCONIC", 12) && osPP.size() > 0 )
        {
            osUnderlyingSRS.Printf( 
                "PROJCS[\"unnamed\",%s,PROJECTION[\"Polyconic\"],PARAMETER[\"latitude_of_origin\",0],PARAMETER[\"central_meridian\",%s],PARAMETER[\"false_easting\",0],PARAMETER[\"false_northing\",0]]", 
                pszGEOGCS, osPP.c_str() );
        }
        
        else if( EQUALN(pszPR,"PR=LAMBERT CONFORMAL CONIC", 26) 
                 && osPP.size() > 0 && pszKNQ != NULL )
        {
            CPLString osP2, osP3;
        
            // Capture the KNQ/P2 string.
            pszValue = strstr(pszKNQ,"P2=");
            if( pszValue )
                pszEnd = strstr(pszValue,",");
            if( pszValue && pszEnd )
                osP2.assign(pszValue+3,pszEnd-pszValue-3);
            
            // Capture the KNQ/P3 string.
            pszValue = strstr(pszKNQ,"P3=");
            if( pszValue )
                pszEnd = strstr(pszValue,",");
            if( pszValue )
            {
                if( pszEnd )
                    osP3.assign(pszValue+3,pszEnd-pszValue-3);
                else
                    osP3.assign(pszValue+3);
            }

            if( osP2.size() > 0 && osP3.size() > 0 )
                osUnderlyingSRS.Printf( 
                    "PROJCS[\"unnamed\",%s,PROJECTION[\"Lambert_Conformal_Conic_2SP\"],PARAMETER[\"standard_parallel_1\",%s],PARAMETER[\"standard_parallel_2\",%s],PARAMETER[\"latitude_of_origin\",0.0],PARAMETER[\"central_meridian\",%s],PARAMETER[\"false_easting\",0.0],PARAMETER[\"false_northing\",0.0]]",
                    pszGEOGCS, osP2.c_str(), osP3.c_str(), osPP.c_str() );

        }
    }

/* -------------------------------------------------------------------- */
/*      If we got an alternate underlying coordinate system, try        */
/*      converting the GCPs to that coordinate system.                  */
/* -------------------------------------------------------------------- */
    if( osUnderlyingSRS.length() > 0 )
    {
        OGRSpatialReference oGeog_SRS, oProjected_SRS;
        OGRCoordinateTransformation *poCT;
        
        oProjected_SRS.SetFromUserInput( osUnderlyingSRS );
        oGeog_SRS.CopyGeogCSFrom( &oProjected_SRS );
        
        poCT = OGRCreateCoordinateTransformation( &oGeog_SRS, 
                                                  &oProjected_SRS );
        if( poCT != NULL )
        {
            for( i = 0; i < nGCPCount; i++ )
            {
                poCT->Transform( 1, 
                                 &(pasGCPList[i].dfGCPX), 
                                 &(pasGCPList[i].dfGCPY), 
                                 &(pasGCPList[i].dfGCPZ) );
            }

            osGCPProjection = osUnderlyingSRS;

            delete poCT;
        }
        else
            CPLErrorReset();
    }

/* -------------------------------------------------------------------- */
/*      Attempt to prepare a geotransform from the GCPs.                */
/* -------------------------------------------------------------------- */
    if( GDALGCPsToGeoTransform( nGCPCount, pasGCPList, adfGeoTransform, 
                                FALSE ) )
    {
        bGeoTransformSet = TRUE;
    }
}

int main( int nArgc, char ** papszArgv )

{
    OGRSpatialReference oSRS;
    int i;
    int bReportXML = FALSE;

    /* -------------------------------------------------------------------- */
    /*      Processing command line arguments.                              */
    /* -------------------------------------------------------------------- */
    nArgc = OGRGeneralCmdLineProcessor( nArgc, &papszArgv, 0 );

    if( nArgc < 2 )
        Usage();

    for( i = 1; i < nArgc; i++ )
    {
        if( EQUAL(papszArgv[i],"-xml") )
            bReportXML = TRUE;

        else if( EQUAL(papszArgv[i],"-t") && i < nArgc - 4 )
        {
            OGRSpatialReference oSourceSRS, oTargetSRS;
            OGRCoordinateTransformation *poCT;
            double                      x, y, z_orig, z;
            int                         nArgsUsed = 4;

            if( oSourceSRS.SetFromUserInput(papszArgv[i+1]) != OGRERR_NONE )
            {
                CPLError( CE_Failure, CPLE_AppDefined,
                          "SetFromUserInput(%s) failed.",
                          papszArgv[i+1] );
                continue;
            }
            if( oTargetSRS.SetFromUserInput(papszArgv[i+2]) != OGRERR_NONE )
            {
                CPLError( CE_Failure, CPLE_AppDefined,
                          "SetFromUserInput(%s) failed.",
                          papszArgv[i+2] );
                continue;
            }

            poCT = OGRCreateCoordinateTransformation( &oSourceSRS,
                    &oTargetSRS );
            x = CPLAtof( papszArgv[i+3] );
            y = CPLAtof( papszArgv[i+4] );
            if( i < nArgc - 5
                    && (CPLAtof(papszArgv[i+5]) > 0.0 || papszArgv[i+5][0] == '0') )
            {
                z_orig = z = CPLAtof(papszArgv[i+5]);
                nArgsUsed++;
            }
            else
                z_orig = z = 0;

            if( poCT == NULL || !poCT->Transform( 1, &x, &y, &z ) )
                printf( "Transformation failed.\n" );
            else
                printf( "(%f,%f,%f) -> (%f,%f,%f)\n",
                        CPLAtof( papszArgv[i+3] ),
                        CPLAtof( papszArgv[i+4] ),
                        z_orig,
                        x, y, z );

            i += nArgsUsed;
        }
        else
        {
            /* coverity[tainted_data] */
            if( oSRS.SetFromUserInput(papszArgv[i]) != OGRERR_NONE )
                CPLError( CE_Failure, CPLE_AppDefined,
                          "Error occurred translating %s.\n",
                          papszArgv[i] );
            else
            {
                char  *pszWKT = NULL;

                if( oSRS.Validate() != OGRERR_NONE )
                    printf( "Validate Fails.\n" );
                else
                    printf( "Validate Succeeds.\n" );

                oSRS.exportToPrettyWkt( &pszWKT, FALSE );
                printf( "WKT[%s] =\n%s\n",
                        papszArgv[i], pszWKT );
                CPLFree( pszWKT );

                printf( "\n" );

                oSRS.exportToPrettyWkt( &pszWKT, TRUE );
                printf( "Simplified WKT[%s] =\n%s\n",
                        papszArgv[i], pszWKT );
                CPLFree( pszWKT );

                printf( "\n" );

                OGRSpatialReference *poSRS2;

                poSRS2 = oSRS.Clone();
                poSRS2->StripCTParms();
//.........这里部分代码省略.........

bool shape::createMesh(OGRSpatialReference* sr2,glm::vec2 origin,glm::vec2 scale,terrain& t)
{
	// Now to pair the points
	vector<Vertex> vertexs = {};
	vector<int> indicies = {};
    rend.init();
    GLuint vertex_shader2;
    GLuint frag_shader2;
    vertex_shader2 = glCreateShader(GL_VERTEX_SHADER);

    frag_shader2 = glCreateShader(GL_FRAGMENT_SHADER);
    cout << "LINE RENDERER CREATOR" << endl;
    if(!rend.addShader(vertex_shader2,"../shaders/linevert.shader"))
    {
      cout << "Failed to add shader" << endl;
      return false;
    }
    if(!rend.addShader(frag_shader2,"../shaders/linefrag.shader"))
    {
        cout << "Failed to add shader" << endl;
        return false;
    }

    if(!rend.compile())
    {
        cout << "Failed to compile shaders" << endl;
        return false;
    }
    if(!rend.link())
    {
        cout << "Failed to link shaders" << endl;
        return false;
    }
	OGRCoordinateTransformation* transform;
    // Create a coordinate transform
    if(sr != NULL)
    {
    	transform = OGRCreateCoordinateTransformation( sr, sr2 );
    }
    else
    {
    	return false;
    }
    /*Vertex tempo = {{0,0,0},{1,0,0},{1,1}};
    Vertex tempo2 = {{0,100,0},{1,0,0},{1,1}};
    Vertex tempo3 = {{100,0,100},{1,0,0},{1,1}};
    Vertex tempo4 = {{100,100,100},{1,0,0},{1,1}};
    vertexs.push_back(tempo);
    vertexs.push_back(tempo2);
    vertexs.push_back(tempo3);
    vertexs.push_back(tempo4);
    indicies.push_back(0);
    indicies.push_back(1);
    indicies.push_back(3);*/
    //indicies.push_back(2);
    //indicies.push_back(3);
    //indicies.push_back(0);
	// triangulate and create lines between vertices
    bool first = true;
	for(int i = 0; i < points.size(); i++)
	{
        //break;
        first = true;
		for(int j = 0; j < points[i].size(); j++)
		{
            if(j == points[i].size() - 1)
            {
                continue;
            }
			double x = points[i][j].x;
            double y = points[i][j].y;
            //points[points.size()-1].push_back(glm::vec2(x,y));
            transform->Transform (1, &x, &y);
            glm::vec2 orig = glm::vec2(x,y);
            x = abs(x -origin.x) *scale.x;
            y = abs(y-origin.y)*scale.y;

			Vertex temp = {{y,t.interpolateHeight(orig),x},{1,0,0},{1,1}};
            //cout << "Vertex: " << temp.position.x << " " << temp.position.y << " " << temp.position.z << endl;
            if(first)
            {
			 vertexs.push_back(temp);
			 temp.position.y += 2;
			 vertexs.push_back(temp);
             first = false;
            }
			if(j+1 < points[i].size())
			{
				x = points[i][j+1].x;
            	y = points[i][j+1].y;
                transform->Transform (1, &x, &y);
                orig = glm::vec2(x,y);
                x = abs(x -origin.x) *scale.x;
                y = abs(y-origin.y)*scale.y;

				Vertex temp2 = {{y,t.interpolateHeight(orig),x},{1,0,0},{1,1}};
                //cout << "Vertex: " << temp2.position.x << " " << temp2.position.y << " " << temp2.position.z << endl;
				vertexs.push_back(temp2);
				temp2.position.y += 2;
				vertexs.push_back(temp2);
//.........这里部分代码省略.........

void wxGISMapView::AddLayer(wxGISLayer* pLayer)
{
	IDisplayTransformation* pDisplayTransformation = pGISScreenDisplay->GetDisplayTransformation();
	if(pDisplayTransformation->GetSpatialReference() == NULL)
	{
		OGRSpatialReference* pSpaRef = pLayer->GetSpatialReference();
		if(pSpaRef)
			pDisplayTransformation->SetSpatialReference(pSpaRef);
		else
		{
			OGREnvelope* pEnv = pLayer->GetEnvelope();
			if(pEnv)
			{
				if(pEnv->MaxX <= 180 && pEnv->MaxY <= 90 && pEnv->MinX >= -180 && pEnv->MinY >= -90)
				{
					OGRSpatialReference* pSpaRef = new OGRSpatialReference();
					pSpaRef->SetWellKnownGeogCS("WGS84");
					pDisplayTransformation->SetSpatialReference(pSpaRef);
				}
			}
		}
	}
	
	OGREnvelope* pEnv = pLayer->GetEnvelope();
	if(pEnv == NULL)
		return;
	OGRSpatialReference* pSpaRef = pLayer->GetSpatialReference();

	if(pSpaRef && pDisplayTransformation->GetSpatialReference())
	{
		if(!pDisplayTransformation->GetSpatialReference()->IsSame(pSpaRef))
		{
			OGRCoordinateTransformation *poCT = OGRCreateCoordinateTransformation( pSpaRef, pDisplayTransformation->GetSpatialReference() );
			poCT->Transform(1, &pEnv->MaxX, &pEnv->MaxY);
			poCT->Transform(1, &pEnv->MinX, &pEnv->MinY);
		}
	}

	if(!pDisplayTransformation->IsBoundsSet())
		pDisplayTransformation->SetBounds(*pEnv);
	else
	{
		OGREnvelope Bounds = pDisplayTransformation->GetBounds();
		Bounds.Merge(*pEnv);
		pDisplayTransformation->SetBounds(Bounds);
	}

	//caching
	if(pLayer->GetCached())
	{
		pLayer->SetCacheID(pGISScreenDisplay->AddCache());
	}
	else
	{
		if(m_Layers.size() > 0 && m_Layers[m_Layers.size() - 1]->GetCached())
			pLayer->SetCacheID(pGISScreenDisplay->AddCache());
		else
			pLayer->SetCacheID(pGISScreenDisplay->GetLastCacheID());
	}
	wxGISMap::AddLayer(pLayer);
}

GDALTiler::GDALTiler(GDALDataset *poDataset, const Grid &grid, const TilerOptions &options):
  mGrid(grid),
  poDataset(poDataset),
  options(options)
{

  // Transformed bounds can give slightly different results on different threads unless mutexed
  static std::mutex mutex;
  std::lock_guard<std::mutex> lock(mutex);

  // if the dataset is set we need to initialise the tile bounds and raster
  // resolution from it.
  if (poDataset != NULL) {

    // Get the bounds of the dataset
    double adfGeoTransform[6];
    CRSBounds bounds;

    if (poDataset->GetGeoTransform(adfGeoTransform) == CE_None) {
      bounds = CRSBounds(adfGeoTransform[0],
                         adfGeoTransform[3] + (poDataset->GetRasterYSize() * adfGeoTransform[5]),
                         adfGeoTransform[0] + (poDataset->GetRasterXSize() * adfGeoTransform[1]),
                         adfGeoTransform[3]);
    } else {
      throw CTBException("Could not get transformation information from source dataset");
    }

    // Find out whether the dataset SRS matches that of the grid
    const char *srcWKT = poDataset->GetProjectionRef();
    if (!strlen(srcWKT))
      throw CTBException("The source dataset does not have a spatial reference system assigned");

    OGRSpatialReference srcSRS = OGRSpatialReference(srcWKT);
    OGRSpatialReference gridSRS = mGrid.getSRS();

    if (!srcSRS.IsSame(&gridSRS)) { // it doesn't match
      // Check the srs is valid
      switch(srcSRS.Validate()) {
      case OGRERR_NONE:
        break;
      case OGRERR_CORRUPT_DATA:
        throw CTBException("The source spatial reference system appears to be corrupted");
        break;
      case OGRERR_UNSUPPORTED_SRS:
        throw CTBException("The source spatial reference system is not supported");
        break;
      default:
        throw CTBException("There is an unhandled return value from `srcSRS.Validate()`");
      }

      // We need to transform the bounds to the grid SRS
      double x[4] = { bounds.getMinX(), bounds.getMaxX(), bounds.getMaxX(), bounds.getMinX() };
      double y[4] = { bounds.getMinY(), bounds.getMinY(), bounds.getMaxY(), bounds.getMaxY() };

      OGRCoordinateTransformation *transformer = OGRCreateCoordinateTransformation(&srcSRS, &gridSRS);
      if (transformer == NULL) {
        throw CTBException("The source dataset to tile grid coordinate transformation could not be created");
      } else if (transformer->Transform(4, x, y) != true) {
        delete transformer;
        throw CTBException("Could not transform dataset bounds to tile spatial reference system");
      }
      delete transformer;

      // Get the min and max values of the transformed coordinates
      double minX = std::min(std::min(x[0], x[1]), std::min(x[2], x[3])),
        maxX = std::max(std::max(x[0], x[1]), std::max(x[2], x[3])),
        minY = std::min(std::min(y[0], y[1]), std::min(y[2], y[3])),
        maxY = std::max(std::max(y[0], y[1]), std::max(y[2], y[3]));

      mBounds = CRSBounds(minX, minY, maxX, maxY); // set the bounds
      mResolution = mBounds.getWidth() / poDataset->GetRasterXSize(); // set the resolution

      // cache the SRS string for use in reprojections later
      char *srsWKT = NULL;
      if (gridSRS.exportToWkt(&srsWKT) != OGRERR_NONE) {
        CPLFree(srsWKT);
        throw CTBException("Could not create grid WKT string");
      }
      crsWKT = srsWKT;
      CPLFree(srsWKT);
      srsWKT = NULL;

    } else {                    // no reprojection is necessary
      mBounds = bounds;         // use the existing dataset bounds
      mResolution = std::abs(adfGeoTransform[1]); // use the existing dataset resolution
    }

    poDataset->Reference();     // increase the refcount of the dataset
  }
}