[Gdal-dev] How can I get coordinate system and projection from shape files?

[lky]

Frank Warmerdam,

Thanks for your help.

I met another question about coordinate system and projection.

I read the data from the shapefile by OGRSFDriverRegistrar::Open() method , then I transform the data by OGRCoordinateTransformation::Transform() method. But the coordinate system and projection of the OGRDataSource is not changed at all.

How can I change the spatial reference?

The following is my program.
/-----------------------------------------------------------------------------------------------------------------------------------------------------
 OGRDataSource*    poDS = OGRSFDriverRegistrar::Open(sourcefile,FALSE);                            //Open Source Data
 for (i=0;i<poDS->GetLayerCount();i++)
 {
  OGRLayer*                    poLayer = poDS->GetLayer(i);                                                                    //Open Layer
  OGRSpatialReference*    pSourceSRS = poLayer->GetSpatialRef();                                                   //Open Spatial Reference
  OGRSpatialReference*    pTargetSRS = pSourceSRS->Clone();                                                         //Define target  Spatial Reference
  pTargetSRS->SetProjCS("TM"); 
  pTargetSRS->SetTM(0,0,1,0,0);
  OGRCoordinateTransformation*    poCT = OGRCreateCoordinateTransformation(pSourceSRS,pTargetSRS);    //Create Coordiante Transform
  poLayer->ResetReading();
 OGRFeature      *poFeature;
  while ((poFeature = poLayer->GetNextFeature()) != NULL)                        //Open Feature
  {
   poGeometry = poFeature->GetGeometryRef();                                            //Open Geometry
   if (poGeometry !=NULL)   
   {
    if (wkbFlatten(poGeometry->getGeometryType()) == wkbPoint)                //point
    {
     OGRPoint *poPoint = (OGRPoint *) poGeometry;                                    //Get coordiante
     x = poPoint->getX();
     y = poPoint->getY();
     if (poCT == NULL || !poCT->Transform(1,&x,&y))                                //transform x,y
     {
      err = 3;   
      return NULL;   
    }
     poPoint->setX(x);                                                                                    //Set coordinate
     poPoint->setY(y);
    }//if
   OGRFeature::DestroyFeature(poFeature);                                                //deal feature
  }//while -- feature loop
 }//for -- layer loop
 return poDS;                                                                                                //return data source
------------------------------------------------------------------------------------------------------------------------------------------------/
Best regards,

Kunyang Li