Reprojecting ShapeFile Data to sync up with OSM-Tiles

Topics: Algorithms
Jan 18, 2011 at 12:36 PM

Hello everyone!

I implemented BruTile to add a Backgroundmap to my programm.

The Requesting and displaying of the Map works just fine but the the Lines I put out from a Shapefile do not sync. (Vienna in front of the affrican coast for example)

I know it is a projection problem and I need to project my ShapeFile into SphericalMercator to have it line up with the map.

The Problem I face with this is that I have no Idea about how to do this.

Any help would be highly appreciated.

Kind Regards

Marc Reinmayr

Coordinator
Jan 19, 2011 at 7:01 AM

Hello Marc,

How reprojection with SharpMap and ProjNet works is described in the Documentation | HowTo section (here). Since reprojection is an expensive task I'd suggest you convert the shapefiles to some other projection using ogr2ogr from FWTools.

If that is not an option, you could compile sharpmap yourself and setting configuration to use DotSpatialProjection instead of ProjNet. You'll find examples in the WinFormSamples project how that works. If you do that, you're required to use .Net 4.0

Hth FObermaier

 

Jan 21, 2011 at 10:00 AM

Thank you very much your help.

I implemented the repojection so far but it seems that my parameters I use for the target CS aren not those for OSM.

I knew that the WKT I used is from Google but I thought it might be the same.

Is there any source of information what WKT is used by OSM so I can reproject my files correctly?

Kind regards

Marc

Marked as answer by FObermaier on 10/23/2013 at 10:33 AM
Developer
Jan 21, 2011 at 10:21 AM
Edited Jan 21, 2011 at 10:21 AM

OSM and GOOGLE uses the same mercator projection named EPSG:3857 (A.K.A. as 900913: http://trac.osgeo.org/openlayers/wiki/SphericalMercator)

http://spatialreference.org/ref/sr-org/6864/

Jan 23, 2011 at 2:26 PM

Hello again!

Thank you for the detailed reply I run some tests and saw that it does not really matter wheter or not I reproject my data, even if I have the .prj File for a Shapefile

 

My Code looks like this:

 

 private bool reprojectLayers()
        {
            bool noProjectionDetected = false;
            //Check if there is any layer without a CS
            for(int i =0; i < mbMap.Map.Layers.Count;i++)
            {
                if( ((mbMap.Map.Layers[i] as VectorLayer).DataSource as ShapeFile).CoordinateSystem == null )
                    noProjectionDetected = true;  
            }
            //Ask User if he wants to continue
            if(noProjectionDetected)
            {
                DialogResult dr = MessageBox.Show("If you continue to reproject, the output might not be usable! Do you want to continue?", "Shapefiles without Projection Data Found", MessageBoxButtons.YesNo);
                if (dr.CompareTo(DialogResult.Yes)==0)
                {
                    for(int i = 0; i < mbMap.Map.Layers.Count;i++)
                    {
                        IGeographicCoordinateSystem sourceCS = (((mbMap.Map.Layers[i] as VectorLayer).DataSource as ShapeFile).CoordinateSystem) as IGeographicCoordinateSystem;
                        if(sourceCS != null)
                            (mbMap.Map.Layers[i] as VectorLayer).CoordinateTransformation = CoordinateSystemManager.ReProjection(sourceCS,CoordinateSystemManager.c_EPSG900913(sourceCS));
                    }
                    return true;
                    MessageBox.Show("Reprojection has been succesfully performed","Reprojection finished!");
                }
            }
            else { MessageBox.Show("Reprojection canceled"); return false; } return false;
           
        }

 

CoordinateSystemManager Class:

//Returns a Projected CS for the target CS for: GoogleMaps and OpenStreetMap
        public static IProjectedCoordinateSystem c_EPSG900913(IGeographicCoordinateSystem sourceCS) // Google Mercator/OSM Mercator
        {
            CoordinateSystemFactory csFac = new CoordinateSystemFactory();
            List<ProjectionParameter> parameters = new List<ProjectionParameter>();
            parameters.Add(new ProjectionParameter("semi_major", 6378137.0));
            parameters.Add(new ProjectionParameter("semi_minor", 6378137.0));
            parameters.Add(new ProjectionParameter("latitude_of_origin", 0.0));
            parameters.Add(new ProjectionParameter("central_meridian", 0.0));
            parameters.Add(new ProjectionParameter("scale_factor", 1.0));
            parameters.Add(new ProjectionParameter("false_easting", 0.0));
            parameters.Add(new ProjectionParameter("false_northing", 0.0));
            IProjection projection = csFac.CreateProjection("OSM Mercator", "mercator_1sp", parameters);

            IProjectedCoordinateSystem epsg900913 = csFac.CreateProjectedCoordinateSystem(
                "OSM Mercator", sourceCS, projection, LinearUnit.Metre,
                new AxisInfo("West", AxisOrientationEnum.West), new AxisInfo("South", AxisOrientationEnum.South));
            return epsg900913;
        }

        //Returns the Transformation Object to apply for the layer
        public static ICoordinateTransformation ReProjection(ICoordinateSystem source, ICoordinateSystem target)
        {
            CoordinateTransformationFactory ctFac = new CoordinateTransformationFactory();
            return ctFac.CreateFromCoordinateSystems(source, target);
        }

 

I have no real experience with CoordinateSystems nor Projection, so I dont know y exactly it is not working. My Idea was simply to get the .prj File with the CS and let it reproject from there to a static epsg900913 CS

I hope someone can help me with this. Thanks in advance!

 

Coordinator
Jan 24, 2011 at 8:41 AM

Hello Marc,

the only time you do not need to reproject your (shapefile) vector data is when your vector data is already in the same coordinatesystem as the osm-tiles, that is epsg:900913 as you defined it.

I'd go about this as follows:

  • Check if shapefile datasource has CoordinateSystem defined. If not, ask user for CoordinateSystem as wkt / or Authority:Code string (EPSG:3857)
    • If wkt then parse that to CoordinateSystem using ProjNet converter
    • If Autority:Code than look that up wkt and parse that.
  • Check if CoordinateSystem is equal to epsg900013 (equalparams)
  • If not setup CoordinateTransformation for this layer using CoordinateTransformation factory.

Keep in mind that CoordinateTransformation is an expensive task, that is done every time you render some portion of the map. Therefore I suggest you do that beforehand and try to avoid it.

Hth FObermaier

Jan 24, 2011 at 10:13 AM

I get the idea now, but it still does not work.

With the above said code it doesnt seem to matter if I reproject.

If I take away the .prj File from the Shapefile I run my tests with the output is in the

middle of the black sea, the same when I add the .prj File and retrieve it.

It is still on the same spot middle of Black Sea

But thank you anyway, I will keep in mind to add a WKT input and a auto-check if the CS is already the same.

But from the files I got for testing, I know it is not the same CS so it at least put it somewhere else.

Thanks a lot!

Marc Reinmayr

Coordinator
Jan 24, 2011 at 10:50 AM

If you can, raise an issue and post a sample shapefile.

cheers FObermaier

Jan 24, 2011 at 12:44 PM

I would like to but the files we are given are confidential

and I cannot give them away. But I am able to upload the content of the .prj File

if that would be any help

Coordinator
Jan 24, 2011 at 1:10 PM

It would need one random geometry and the contents of the prj file. The attributes (*.dbf) are not required.

It might help if you state where that geometry should be.

 

 

Jan 24, 2011 at 1:35 PM

http://www.easy-share.com/1913657677/test.rar

These are the .shp and the .prj file if I am correct it should be in Germany when placed correctly.

Coordinator
Jan 25, 2011 at 9:59 AM
Edited Jun 22, 2011 at 12:15 PM

Maybe the problem is that the coordinate system of the shapefile is typeof(IProjectedCoordinateSystem).

This code reprojects the layer, but it needs fine tuning:

 

    var ctf = new CoordinateTransformationFactory();
    var epsg3857 = new CoordinateSystemFactory().CreateFromWkt(
        "PROJCS[\"Popular Visualisation CRS / Mercator\", "+
                 "GEOGCS[\"Popular Visualisation CRS\", "+
                         "DATUM[\"Popular Visualisation Datum\", "+
                                 "SPHEROID[\"Popular Visualisation Sphere\", 6378137, 0, AUTHORITY[\"EPSG\",\"7059\"]], "+
                                 "TOWGS84[0, 0, 0, 0, 0, 0, 0], "+
                                 "AUTHORITY[\"EPSG\",\"6055\"]], "+
                         "PRIMEM[\"Greenwich\", 0, AUTHORITY[\"EPSG\", \"8901\"]], "+
                         "UNIT[\"degree\", 0.0174532925199433, AUTHORITY[\"EPSG\", \"9102\"]], "+
                         "AXIS[\"E\", EAST], AXIS[\"N\", NORTH], "+
                         "AUTHORITY[\"EPSG\",\"4055\"]], "+
                 "PROJECTION[\"Mercator\"], "+
                              "PARAMETER[\"False_Easting\", 0], "+
                              "PARAMETER[\"False_Northing\", 0], "+
                              "PARAMETER[\"Central_Meridian\", 0], "+
                              "PARAMETER[\"Latitude_of_origin\", 0], "+
                              "UNIT[\"metre\", 1, AUTHORITY[\"EPSG\", \"9001\"]], "+
                              "AXIS[\"East\", EAST], AXIS[\"North\", NORTH], "+
                              "AUTHORITY[\"EPSG\",\"3857\"]]");

    Map map = new Map();
    map.Layers.Add(new TileLayer(new OsmTileSource(), "OSM"));
    var sf = new ShapeFile("GeoData\\World\\gasleitung.shp", false);
    sf.Open();
    var vl = new VectorLayer("Gasleitung", sf) {Style = {Line = new Pen(Brushes.Black, 2)}};
    vl.CoordinateTransformation = ctf.CreateFromCoordinateSystems(
        sf.CoordinateSystem, epsg3857);
    map.Layers.Add(vl);
    map.ZoomToBox(vl.Envelope);

Jan 25, 2011 at 12:43 PM

Thank you so much for your help!

I finally understood were my mistake was and changed the code accordingly and now It works!

With kind regards

Marc Reinmayr

Editor
Oct 23, 2013 at 1:40 PM
I just wanted to record a few of my most recent findings here for anyone who is looking at a similar problem to me.
I have an web based tool that uses openlayers to allow a user to digitise points over an EPSG:900913 projected base map.

I found that when I overlaid those points on an osm basemap in sharpmap, using the above projection transformation they were out by about 20km. after investigating I found this post http://alastaira.wordpress.com/2011/01/23/the-google-maps-bing-maps-spherical-mercator-projection/ which pointed me towards the solution.

So by specifying the semi-minor axis as suggested, it all overlaid perfectly. So the code now looks like this for me:

private ICoordinateSystem GetGoogleProjection()
    {
    var epsg3857 = new CoordinateSystemFactory().CreateFromWkt(
    "PROJCS[\"Popular Visualisation CRS / Mercator\"," +
    "GEOGCS[\"Popular Visualisation CRS\"," +
    "DATUM[\"WGS84\"," +
    "SPHEROID[\"WGS84\", 6378137.0, 298.257223563, AUTHORITY[\"EPSG\",\"7059\"]]," +
    "AUTHORITY[\"EPSG\",\"6055\"]]," +
    "PRIMEM[\"Greenwich\", 0, AUTHORITY[\"EPSG\", \"8901\"]]," +
    "UNIT[\"degree\", 0.0174532925199433, AUTHORITY[\"EPSG\", \"9102\"]]," +
    "AXIS[\"E\", EAST], AXIS[\"N\", NORTH], AUTHORITY[\"EPSG\",\"4055\"]]," +
    "PROJECTION[\"Mercator\"]," +
    "PARAMETER[\"semi_minor\",6378137]," +
    "PARAMETER[\"False_Easting\", 0]," +
    "PARAMETER[\"False_Northing\", 0]," +
    "PARAMETER[\"Central_Meridian\", 0]," +
    "PARAMETER[\"Latitude_of_origin\", 0]," +
    "UNIT[\"metre\", 1, AUTHORITY[\"EPSG\", \"9001\"]]," +
    "AXIS[\"East\", EAST], AXIS[\"North\", NORTH]," +
    "AUTHORITY[\"EPSG\",\"3785\"]]");

        return epsg3857;
    }
Developer
Oct 23, 2013 at 2:48 PM
Edited Oct 23, 2013 at 2:54 PM
here the code used in WMS.Demo to make programmatically an EPSG:3857:

this code works well using Proj.Net, when using WKT an exception is thrown.
Coordinator
Oct 23, 2013 at 3:59 PM
What is wrong with
ProjNet.CoordinateSystems.Projected.WebMercator
Developer
Oct 24, 2013 at 7:03 AM
Edited Oct 24, 2013 at 7:05 AM
FObermaier wrote:
What is wrong with
ProjNet.CoordinateSystems.Projected.WebMercator
Actually, nothing apart that I see this class for the first time :)
I fix immediately the SharpMap.Demo project to use this stuff.

EDIT: changeset 104305
Editor
Oct 25, 2013 at 9:52 AM
I was following the code in the how to and wasn't aware of the enumerate coordinate systems. that might make things easier!
Nov 22, 2014 at 11:59 AM
displaying of the Map works just fine when shape file is fetched from harddisk using the code mentioned below.


But now I want to fetch shapefile from postgis (geom column).
Not able to figure out how to modify following statement.
vl.CoordinateTransformation = ctf.CreateFromCoordinateSystems( sf.CoordinateSystem, epsg3857);

is it possible to fetch only single column of a record using postgis provider?
And is it possible to convert geometry column fetched from postgis to shapefile?



var ctf = new CoordinateTransformationFactory();
var epsg3857 = new CoordinateSystemFactory().CreateFromWkt(
    "PROJCS[\"Popular Visualisation CRS / Mercator\", "+
             "GEOGCS[\"Popular Visualisation CRS\", "+
                     "DATUM[\"Popular Visualisation Datum\", "+
                             "SPHEROID[\"Popular Visualisation Sphere\", 6378137, 0, AUTHORITY[\"EPSG\",\"7059\"]], "+
                             "TOWGS84[0, 0, 0, 0, 0, 0, 0], "+
                             "AUTHORITY[\"EPSG\",\"6055\"]], "+
                     "PRIMEM[\"Greenwich\", 0, AUTHORITY[\"EPSG\", \"8901\"]], "+
                     "UNIT[\"degree\", 0.0174532925199433, AUTHORITY[\"EPSG\", \"9102\"]], "+
                     "AXIS[\"E\", EAST], AXIS[\"N\", NORTH], "+
                     "AUTHORITY[\"EPSG\",\"4055\"]], "+
             "PROJECTION[\"Mercator\"], "+
                          "PARAMETER[\"False_Easting\", 0], "+
                          "PARAMETER[\"False_Northing\", 0], "+
                          "PARAMETER[\"Central_Meridian\", 0], "+
                          "PARAMETER[\"Latitude_of_origin\", 0], "+
                          "UNIT[\"metre\", 1, AUTHORITY[\"EPSG\", \"9001\"]], "+
                          "AXIS[\"East\", EAST], AXIS[\"North\", NORTH], "+
                          "AUTHORITY[\"EPSG\",\"3857\"]]");

Map map = new Map();
map.Layers.Add(new TileLayer(new OsmTileSource(), "OSM"));
var sf = new ShapeFile("GeoData\\World\\gasleitung.shp", false);
sf.Open();
var vl = new VectorLayer("Gasleitung", sf) {Style = {Line = new Pen(Brushes.Black, 2)}};
vl.CoordinateTransformation = ctf.CreateFromCoordinateSystems(
    sf.CoordinateSystem, epsg3857);
map.Layers.Add(vl);
map.ZoomToBox(vl.Envelope);

Any help would be highly appreciated.
Coordinator
Nov 22, 2014 at 7:46 PM
The geometry_columns table/view together with the spatial_ref_sys table provide the information you are seeking.
The geometry_columns provide the SRID for the table, the st_text column of the spatial_ref_sys table provide the wkt.

If I were you, I'd create a view that transforms your data to the desired projection and use that as source for the postgis provider.