C++ OGRCoordinateTransformation类代码示例

/** 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 );
}

// [[Rcpp::export]]
Rcpp::List CPL_transform(Rcpp::List sfc, Rcpp::CharacterVector proj4) {

	// import proj4string:
	OGRSpatialReference *dest = new OGRSpatialReference;
	handle_error(dest->importFromProj4((const char *) (proj4[0])));

	// transform geometries:
	std::vector<OGRGeometry *> g = ogr_from_sfc(sfc, NULL);
	if (g.size() == 0) {
		dest->Release(); // #nocov
		Rcpp::stop("CPL_transform: zero length geometry list"); // #nocov
	}
	OGRCoordinateTransformation *ct = 
		OGRCreateCoordinateTransformation(g[0]->getSpatialReference(), dest);
	if (ct == NULL) {
		dest->Release(); // #nocov
		Rcpp::stop("OGRCreateCoordinateTransformation() returned NULL: PROJ.4 available?"); // #nocov
	}
	for (size_t i = 0; i < g.size(); i++) {
		CPLPushErrorHandler(CPLQuietErrorHandler);
		OGRErr err = 0;
		if (! g[i]->IsEmpty())
			err = g[i]->transform(ct);
		CPLPopErrorHandler();
		if (err == 1 || err == 6) {
			OGRwkbGeometryType geomType = g[i]->getGeometryType();
			OGRGeometryFactory f;
			f.destroyGeometry(g[i]);
			g[i] = f.createGeometry(geomType);
		} else
			handle_error(err);
	}

	Rcpp::List ret = sfc_from_ogr(g, true); // destroys g;
	ct->DestroyCT(ct);
	dest->Release();
	return ret; 
}

/**
*@brief 点对点的转换
*@param poDataset [in] gdal数据集
*@param xyPeak [in] 左上顶点，和右下顶点的坐标{xl,yl,xr,yr}
*@param band [in] 波段数
*@param typeParam [in] 类型重载变量
*@return
*/
double *GdalProjection::PointToPoint(const char *fromWkt,const char *toWkt,double x,double y)
{
    OGRSpatialReference oSourceSRS,oTargetSRS;
    OGRCoordinateTransformation *poct;

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

    poct = OGRCreateCoordinateTransformation(&oSourceSRS,&oTargetSRS);
    double *xy = (double *)malloc(sizeof(double) * 2);
    if (poct == NULL || !poct->Transform(1,&x,&y))
    {
        xy[0] = 0.0;
        xy[1] = 0.0;
    }
    else
    {
        xy[0] = x;
        xy[1] = y;
    }
    delete poct;
    return xy;
}

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;
}

int main( int nArgc, char ** papszArgv )
{
    // register drivers
    GDALAllRegister();

    if( nArgc < 2 )
        return EXIT_FAILURE;

    double dfaCornersX[5] = {0};
    double dfaCornersY[5] = {0};
    CPLString sFileName;

    // parse input values
    for( int iArg = 1; iArg < nArgc; iArg++ )
    {
        if( EQUAL(papszArgv[iArg],"-nw"))
        {
            CHECK_HAS_ENOUGH_ADDITIONAL_ARGS(2);
            const char* pszCoord = papszArgv[++iArg];
            dfaCornersY[1] = CPLAtofM(pszCoord);
            pszCoord = papszArgv[++iArg];
            dfaCornersX[1] = CPLAtofM(pszCoord);
        }
        else if( EQUAL(papszArgv[iArg],"-ne"))
        {
            CHECK_HAS_ENOUGH_ADDITIONAL_ARGS(2);
            const char* pszCoord = papszArgv[++iArg];
            dfaCornersY[2] = CPLAtofM(pszCoord);
            pszCoord = papszArgv[++iArg];
            dfaCornersX[2] = CPLAtofM(pszCoord);
        }
        else if( EQUAL(papszArgv[iArg],"-se"))
        {
            CHECK_HAS_ENOUGH_ADDITIONAL_ARGS(2);
            const char* pszCoord = papszArgv[++iArg];
            dfaCornersY[3] = CPLAtofM(pszCoord);
            pszCoord = papszArgv[++iArg];
            dfaCornersX[3] = CPLAtofM(pszCoord);
        }
        else if( EQUAL(papszArgv[iArg],"-sw"))
        {
            CHECK_HAS_ENOUGH_ADDITIONAL_ARGS(2);
            const char* pszCoord = papszArgv[++iArg];
            dfaCornersY[4] = CPLAtofM(pszCoord);
            pszCoord = papszArgv[++iArg];
            dfaCornersX[4] = CPLAtofM(pszCoord);
        }
        else if( EQUAL(papszArgv[iArg],"-c"))
        {
            CHECK_HAS_ENOUGH_ADDITIONAL_ARGS(2);
            const char* pszCoord = papszArgv[++iArg];
            dfaCornersY[0] = CPLAtofM(pszCoord);
            pszCoord = papszArgv[++iArg];
            dfaCornersX[0] = CPLAtofM(pszCoord);
        }
        else if(sFileName.empty())
            sFileName = papszArgv[iArg];
    }

    OGRSpatialReference oOGRSpatialReference(SRS_WKT_WGS84);
    int nZoneNo = ceil( (180.0 + dfaCornersX[0]) / 6.0 );
    OGRSpatialReference oDstSpatialReference(SRS_WKT_WGS84);
    oDstSpatialReference.SetUTM(nZoneNo, dfaCornersY[0] > 0);

    // transform coordinates from WGS84 to UTM
    OGRCoordinateTransformation *poCT = OGRCreateCoordinateTransformation( &oOGRSpatialReference, &oDstSpatialReference);
    if(!poCT)
    {
        Usage("get coordinate transformation failed");
        return EXIT_FAILURE;
    }

    int nResult = poCT->Transform(5, dfaCornersX, dfaCornersY, NULL);
    if(!nResult)
    {
        Usage("transformation failed");
        return EXIT_FAILURE;
    }

    // open input dataset
    GDALDataset *poSrcDataset = (GDALDataset *) GDALOpen( sFileName, GA_ReadOnly ); // GA_Update
    char* pszSpaRefDef = NULL;
    if( oDstSpatialReference.exportToWkt(&pszSpaRefDef) != OGRERR_NONE)
    {
        CPLFree( pszSpaRefDef );
        GDALClose( (GDALDatasetH) poSrcDataset );
        return EXIT_FAILURE;
    }

    // search point along image
    // add GCP to opened raster
    OGRPoint ptCenter(dfaCornersX[0], dfaCornersY[0]);
    OGRPoint pt1(dfaCornersX[1], dfaCornersY[1]); // NW Cormer
    OGRPoint pt2(dfaCornersX[2], dfaCornersY[2]); // NE Corner
    OGRPoint pt3(dfaCornersX[3], dfaCornersY[3]); // SE Corner
    OGRPoint pt4(dfaCornersX[4], dfaCornersY[4]); // SW Corner
    int nGCPCount = 0;
    OGREnvelope DstEnv;
    GDAL_GCP *paGSPs = PrepareGCP(sFileName, &pt1, &pt2, &pt3, &pt4, &ptCenter, oDstSpatialReference, poSrcDataset->GetRasterXSize(), poSrcDataset->GetRasterYSize(), nGCPCount, DstEnv);

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

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
  }
}

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);
//.........这里部分代码省略.........

std::string FetchTimeZone( double dfX, double dfY, const char *pszWkt )
{
    CPLDebug( "WINDNINJA", "Fetching timezone for  %lf,%lf", dfX, dfY );
    if( pszWkt != NULL )
    {
        OGRSpatialReference oSourceSRS, oTargetSRS;
        OGRCoordinateTransformation *poCT;

        oSourceSRS.SetWellKnownGeogCS( "WGS84" );
        oTargetSRS.importFromWkt( (char**)&pszWkt );

        poCT = OGRCreateCoordinateTransformation( &oSourceSRS, &oTargetSRS );
        if( poCT == NULL )
        {
            CPLError( CE_Failure, CPLE_AppDefined,
                      "OGR coordinate transformation failed" );
            return std::string();
        }
        if( !poCT->Transform( 1, &dfX, &dfY ) )
        {
            CPLError( CE_Failure, CPLE_AppDefined,
                      "OGR coordinate transformation failed" );
            return std::string();
        }
        OGRCoordinateTransformation::DestroyCT( poCT );
    }
    OGRGeometryH hGeometry = OGR_G_CreateGeometry( wkbPoint );
    OGR_G_SetPoint_2D( hGeometry, 0, dfX, dfY );

    OGRDataSourceH hDS;
    OGRLayerH hLayer;
    OGRFeatureH hFeature;

    std::string oTzFile = FindDataPath( "tz_world.zip" );
    oTzFile = "/vsizip/" + oTzFile + "/world/tz_world.shp";

    hDS = OGROpen( oTzFile.c_str(), 0, NULL );
    if( hDS == NULL )
    {
        CPLError( CE_Failure, CPLE_AppDefined,
                  "Failed to open datasource: %s", oTzFile.c_str() );
        return std::string();
    }
    hLayer = OGR_DS_GetLayer( hDS, 0 );
    OGR_L_SetSpatialFilter( hLayer, hGeometry );
    OGR_L_ResetReading( hLayer );
    int nMaxTries = 5;
    int nTries = 0;
    OGRGeometryH hBufferGeometry;
    do
    {
        if( nTries == 0 )
        {
            hBufferGeometry = OGR_G_Clone( hGeometry );
        }
        else
        {
            hBufferGeometry = OGR_G_Buffer( hGeometry, 0.2 * nTries, 30 );
        }
        OGR_L_SetSpatialFilter( hLayer, hBufferGeometry );
        hFeature = OGR_L_GetNextFeature( hLayer );
        OGR_G_DestroyGeometry( hBufferGeometry );
        nTries++;
    }
    while( hFeature == NULL && nTries < nMaxTries );
    std::string oTimeZone;
    if( hFeature == NULL )
    {
        oTimeZone = std::string();
        CPLError( CE_Failure, CPLE_AppDefined,
                  "Failed to find timezone" );
    }
    else
    {
        oTimeZone = std::string( OGR_F_GetFieldAsString( hFeature, 0 ) );
    }
    OGR_F_Destroy( hFeature );
    OGR_G_DestroyGeometry( hGeometry );
    OGR_DS_Destroy( hDS );
    return oTimeZone;
}

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);
}