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

[lky]

Hello,everyone:

I'm a new learner of the gdal. During my study, I met some problems. I hope I can find help here.

I found these is one sentence "If a .prj files in old Arc/Info style or new ESRI OGC WKT style is present, it will be read and used to associate a projection with features." in \ogr\dev_shapefile.html.

And I found there are some *.prj¡¡files with the same name shapefile in directory C:\arcgis\arcexe83\Reference Systems(The directory which I install arcinfo).In the ArcMap, that the frame map file is *.mxt.

I try to use the OGRSFDriverRegistrar::Open function to open these *.prj files but I failed.

So I have the following questions.

1. Who can explain the use of these *.prj files in this directory? Which shapefiles do they support? Or they only support project mode for AcrMap?
2. What is the relation with the *.prj and *mxt files?
3. How can I read the coordinate system and projection from the *.prj files? And whether this prj file support the coordinate system and projection for the same name shapefile or for the shapefiles in the same directory?

Thanks a lot.

Kunyang Li
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://lists.osgeo.org/pipermail/gdal-dev/attachments/20051031/f40cdde4/attachment.html

[Frank Warmerdam]

Post by lky
I'm a new learner of the gdal. During my study, I met some problems. I hope
I can find help here.
I found these is one sentence "If a .prj files in old Arc/Info style or new
ESRI OGC WKT style is present, it will be read and used to associate a
projection with features." in \ogr\dev_shapefile.html.
And I found there are some *.prj?files with the same name shapefile in
directory C:\arcgis\arcexe83\Reference Systems(The directory which I install
arcinfo).In the ArcMap, that the frame map file is *.mxt.
I try to use the OGRSFDriverRegistrar::Open function to open these *.prj
files but I failed.
So I have the following questions.
1. Who can explain the use of these *.prj files in this directory? Which
shapefiles do they support? Or they only support project mode for AcrMap?

Kunyang Li,

You should open the shapefile normally with OGRSFDriverRegistrar::Open()
and then fetch the layers from the datasource. To find the coordinate
system of a layer you can call OGRLayer::GetSpatialRef() to get a handle
to the corresponding OGRSpatialReference object.

This method will read it from a correspondingly named .prj file if it
is present.

Post by lky
2. What is the relation with the *.prj and *mxt files?

I don't know anything about mxt files.

Post by lky
3. How can I read the coordinate system and projection from the *.prj files?
And whether this prj file support the coordinate system and projection for
the same name shapefile or for the shapefiles in the same directory?

OGR assumes a .prj file applies to the shapefile with the same
base filename. If there is no corresponding .prj file, then GetSpatialRef()
will return NULL.

Best regards,
--
---------------------------------------+--------------------------------------
I set the clouds in motion - turn up | Frank Warmerdam, ***@pobox.com
light and sound - activate the windows | http://pobox.com/~warmerdam
and watch the world go round - Rush | Geospatial Programmer for Rent

[lky]

Frank Warmerdam,

Thanks for your help.

I have understood how to get the coordinate system and projection from .prj file. But I have another question about this. Since I can get the coordinate system and projection from .prj file, I also can write these information to it. But I didn't find the support about this. Can you tell me how to write the files?

Best regards

Kunyang Li

----- Original Message -----
From: "Frank Warmerdam" <***@pobox.com>
To: "lky" <***@163.com>
Cc: "gdal ????" <gdal-***@lists.maptools.org>
Sent: Monday, October 31, 2005 11:03 PM
Subject: Re: [Gdal-dev] How can I get coordinate system and projection from shape files?

Post by Frank Warmerdam

Post by lky
I'm a new learner of the gdal. During my study, I met some problems. I hope
I can find help here.
I found these is one sentence "If a .prj files in old Arc/Info style or new
ESRI OGC WKT style is present, it will be read and used to associate a
projection with features." in \ogr\dev_shapefile.html.
And I found there are some *.prj?files with the same name shapefile in
directory C:\arcgis\arcexe83\Reference Systems(The directory which I install
arcinfo).In the ArcMap, that the frame map file is *.mxt.
I try to use the OGRSFDriverRegistrar::Open function to open these *.prj
files but I failed.
So I have the following questions.
1. Who can explain the use of these *.prj files in this directory? Which
shapefiles do they support? Or they only support project mode for AcrMap?

Kunyang Li,
You should open the shapefile normally with OGRSFDriverRegistrar::Open()
and then fetch the layers from the datasource. To find the coordinate
system of a layer you can call OGRLayer::GetSpatialRef() to get a handle
to the corresponding OGRSpatialReference object.
This method will read it from a correspondingly named .prj file if it
is present.

Post by lky
2. What is the relation with the *.prj and *mxt files?

I don't know anything about mxt files.

Post by lky
3. How can I read the coordinate system and projection from the *.prj files?
And whether this prj file support the coordinate system and projection for
the same name shapefile or for the shapefiles in the same directory?

OGR assumes a .prj file applies to the shapefile with the same
base filename. If there is no corresponding .prj file, then GetSpatialRef()
will return NULL.
Best regards,
--
---------------------------------------+--------------------------------------
light and sound - activate the windows | http://pobox.com/~warmerdam
and watch the world go round - Rush | Geospatial Programmer for Rent

[Frank Warmerdam]

Post by Frank Warmerdam
Frank Warmerdam,
Thanks for your help.
I have understood how to get the coordinate system and projection from .prj file. But I have another question about this. Since I can get the coordinate system and projection from .prj file, I also can write these information to it. But I didn't find the support about this. Can you tell me how to write the files?

Kunyang Li,

The only way to write the projection to a shapefile layer
through OGR is to pass it at the point the layer is created.
That is, pass the appropriate OGRSpatialReference into the
OGRDataSource::CreateLayer() call.

If you really need to generate a .prj file for a shapefile later
you can do it manually. To do that, setup an appropriate
OGRSpatialReference, then invoke morphToESRI() on it
to convert the format to match what ESRI software expects.
Then exportToWkt() and manually write the resulting string
to the .prj file yourself.

Best regards,
--
---------------------------------------+--------------------------------------
I set the clouds in motion - turn up | Frank Warmerdam, ***@pobox.com
light and sound - activate the windows | http://pobox.com/~warmerdam
and watch the world go round - Rush | Geospatial Programmer for Rent

[lky]

Frank Warmerdam:

Thanks for your bounteous help.

Using the methods you teld me, I have finished the program that I read the data from one shapefile ,and wrote it to the other one after I transformed their coordinate. But there still some problems in the program. I can use the "transform" method to transform the obejct of class OGRMultiPolygon. But I can't give the coordinate system and projection which I order to the OGRDataSource pointer.

Here, poSourceSRS is a pointer of the coordinate system and projection of the source data,and poTargetSRS is a pointer of the coordinate system and projection of the target data.I have tried two ways, but both failed.Firstly, I got the poSourceSRS by the OGRLayer::GetSpatialRef() method, but I can't change this pointer. Secondly I use the TargetSRS in the OGRDataSource::CreateLayer() method to specify the coordinate system and projection, I will get the result. But this is not useful in the shapefile. When I open the shapefile by ArcMap, the note pop up that is "Note: one or more layers is missing spatial reference information. Data from these layers cannot be projected.".

Perhaps I can't let you clear my meaning because my poor English. I repeat it again.

I transform one data geometry by changing its coordinate system and projection. I have transform the geometry, but I don't know how to write the new coordinate system and projection to the .prj file.

Best Regards

Kunyang Li
----- Original Message -----
From: "Frank Warmerdam" <***@pobox.com>
To: "lky" <***@163.com>
Cc: "gdal ????" <gdal-***@lists.maptools.org>
Sent: Tuesday, November 01, 2005 10:44 PM
Subject: Re: [Gdal-dev] How can I get coordinate system and projection from shape files?

Post by Frank Warmerdam
Frank Warmerdam,
Thanks for your help.
I have understood how to get the coordinate system and projection from .prj file. But I have another question about this. Since I can get the coordinate system and projection from .prj file, I also can write these information to it. But I didn't find the support about this. Can you tell me how to write the files?

Kunyang Li,

The only way to write the projection to a shapefile layer
through OGR is to pass it at the point the layer is created.
That is, pass the appropriate OGRSpatialReference into the
OGRDataSource::CreateLayer() call.

If you really need to generate a .prj file for a shapefile later
you can do it manually. To do that, setup an appropriate
OGRSpatialReference, then invoke morphToESRI() on it
to convert the format to match what ESRI software expects.
Then exportToWkt() and manually write the resulting string
to the .prj file yourself.

Best regards,
--
---------------------------------------+--------------------------------------
I set the clouds in motion - turn up | Frank Warmerdam, ***@pobox.com
light and sound - activate the windows | http://pobox.com/~warmerdam
and watch the world go round - Rush | Geospatial Programmer for Rent

[Frank Warmerdam]

Post by lky
Here, poSourceSRS is a pointer of the coordinate system and projection of the source data,and poTargetSRS is a pointer of the coordinate system and projection of the target data.I have tried two ways, but both failed.Firstly, I got the poSourceSRS by the OGRLayer::GetSpatialRef() method, but I can't change this pointer. Secondly I use the TargetSRS in the OGRDataSource::CreateLayer() method to specify the coordinate system and projection, I will get the result. But this is not useful in the shapefile. When I open the shapefile by ArcMap, the note pop up that is "Note: one or more layers is missing spatial reference information. Data from these layers cannot be projected.".

Kunyang Li,

OK, how did you prepare poTargetSRS? If you pass a pointer to a valid
coordinate system to CreateLayer() it should get written out as a .prj
file (after morphing to ESRI format). Does a .prj file get written out?
Perhaps it is just incompatible with the ArcMap software in someway
I will need to correct.

BTW, what you want to do should also be easily accomplished
at the commandline with "ogr2ogr".

Best regards,
--
---------------------------------------+--------------------------------------
I set the clouds in motion - turn up | Frank Warmerdam, ***@pobox.com
light and sound - activate the windows | http://pobox.com/~warmerdam
and watch the world go round - Rush | Geospatial Programmer for Rent

[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

[Frank Warmerdam]

Post by lky
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.

Kunyang Li,

You just transformed the point location, you didn't do anything to
change the coordinate system of the layer as a whole. You also
didn't do anything to actually write the results back out to the original
layer (if that was your intention). I'm not sure that in-place update is
even supported by the Shapefile driver.

Post by lky
How can I change the spatial reference?

o Create a new layer for the output.
o Make sure you pass the altered OGRSpatialReference to the
CreateLayer() call.
o You might as well just call transform() on the geometry
instead of extracting the point and transforming that.

Make sure you follow the translation example on the OGR page.
It is important to create create all the fields on the new layer
with CreateLayer(), and to make sure the feature get converted
from the old layer definition to the new layer definition (even
though they seem the same).

Best regards,
--
---------------------------------------+--------------------------------------
I set the clouds in motion - turn up | Frank Warmerdam, ***@pobox.com
light and sound - activate the windows | http://pobox.com/~warmerdam
and watch the world go round - Rush | Geospatial Programmer for Rent