Help with projections and transformations

Topics: Algorithms, SharpMap Project, SharpMap v0.9 / v1.x
Mar 6, 2010 at 7:42 AM

Dear John/Fobermaier,

For now, I have decided to not work on SharpMap v2.0 because there are a lot of gaps here and there. For example the MapViewControl throws a NotImplementedException for the MouseWheel event. There are a lot of such minor issues and hence we are kind of wary to jump into the v2.0 bandwagon at the moment.

Now, coming to the point, I am able to load some polygons from the database. I am using GeometryFeatureProvider to do this. The SRID of the layer is set to 4326 (WGS 84). I have also been successful in creating another layer and setting it's DataSource to ShapeFile.

But the problem is the two layers are way off each other. To give you the perspective, the ShapeFile contains polygons pertaining to Land Use patterns as of 2003 for an entire district. The polygons I create using the data in the database has data for current land usage pattern. Now I would like to overlay the current land usage pattern over 2003's land usage pattern.

But the problem is, the polygons from the shape file show up on top-right corner and database polygons show in bottom-left corner of the map. They were supposed be on top of each other.

The SRID of the polygons in the database are always constant (WGS 84). But the shape files we use might be different. I am not even able to figure out how do I find out the projection used in a shape file.

I know this is something to do with the Projections and Transformations. But any reading material, sample source code you could point me to will be of great help.

Thanks,
Raghu

Coordinator
Mar 6, 2010 at 9:37 AM
Edited Mar 6, 2010 at 9:40 AM

Hi Raghu, first you need to decide on the spatial reference system you want for the map. Lets assume WGS84.

We can look the wkt up on http://spatialreference.org which gives us:

 

GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],
AUTHORITY["EPSG","6326"]],
PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.01745329251994328,AUTHORITY["EPSG","9122"]],
AUTHORITY["EPSG","4326"]]
(note I have added random line breaks so that it displays ok on codeplex)

 

Because the data in the database is also in WGS84 we do not need to transform the database layers at all so the data can be added without change.

For the Shapefile layers we will need to transform them into the same spatial reference system we do this by creating an instance of an ICoordinateTransformation.

See http://sharpmap.codeplex.com/wikipage?title=Apply%20on%20the%20fly%20transformation%20to%20a%20layer&referringTitle=How%20to...

Basically we reference ProjNet from the External References folder.

Then we create an ICoordinateSystem by parsing the wkt for WGS84 (see link above).

Shapefile provider already has a CoordinateSystem property which is fed by the .prj file associated with the .shp file - we do not need to do anything special here.

We then create the ICoordinateTransform which converts geometries from the the source projection (i.e the one from the shapefile) to the target (i.e the one from the WGS84 wkt) also detailed in the link above.

Then we set the CoordinateTransformation property on the Layer that the Shapefile feeds.

That should fix up your projections.

The main issue that you may come across is if there is no .prj file accompanying the .shp file. In this case the CoordinateSystem property on the Shapefile provider will be null and you will have to find out what spatial reference system the data actually belongs in, locate and parse the wkt for the srs into a suitable ICoordinateSystem before creating the transform.

hth jd

Mar 6, 2010 at 11:30 AM

Hi John,

Thanks for the detailed explanation.

You're right! The SHP files I have in hand right now do not have a PRJ file! I'll try to get them from the source.

Anyways, this should keep me rolling :)

Thank You
Raghu