Osm layer and SHP file coordinate conversion problem

Topics: General Topics, SharpMap v0.9 / v1.x
Feb 23, 2011 at 10:14 AM

Hi all, in my application I need to have a background OpenStreetMap (or google maps) layer and another layer from a SHP file over it.

After taking some good advices from other posts, I have a problem with the coordinate system conversion.

The result is that I have a good background layer and shp layer as well, but there is a distance of some kilometers between them. I had such kind of problems in previous GPS applications and those were related to decimal/sexagesimal representation in the WGS84 format.

A snapshot of this problem can be seen at: https://docs.google.com/document/pub?id=1Ke9e_AO1ZKfx5LcLYy3XQOp_NkFadAgMe2iPS3-ODc8

Has anyone had the same problem?

Here is my code:

using System;
using System.Collections.Generic;
using System.ComponentModel;
using System.Data;
using System.Drawing;
using System.Linq;
using System.Text;
using System.Windows.Forms;

using System.Configuration;
using Fep.Gui.Tests;
using SharpMap;
using SharpMap.Layers;
using SharpMap.Data.Providers;
using SharpMap.Geometries;
using System.IO;



using BruTile;
using BruTile.Web;
using BruTile.FileSystem;
using BruTile.PreDefined;
using BruTile.Cache;



using ProjNet;
using ProjNet.CoordinateSystems;
using ProjNet.Converters;
using ProjNet.CoordinateSystems.Transformations;


using System.Net;

//using SharpMap.Data.Providers
namespace Fep.Gui.Tests.Maps
{
    public partial class TestMaps : Form
    {
        public TestMaps()
        {
            InitializeComponent();
            InitializeMap();
        }

        #region Event handlers
        
        private void rbnQuitHandler(object sender, EventArgs e)
        {
            if (MessageBox.Show("Quit?", "Quit", MessageBoxButtons.YesNo) == DialogResult.Yes)
            {
                this.Close();
            }
        }
        #endregion
        #region Map Initialization
        /// <summary>
        /// Inizializza il controllo con la mappa
        /// </summary>
        private void InitializeMap() 
        {
           //TestMapForm.Default.AbuDabhi_administrative
            Map map = new Map(splitContainer1.Panel2.ClientSize);

            string mapFileName;
            mapFileName = TestMapForm.Default.Default_AbuDabhi_ShpFile;
            if (!File.Exists(mapFileName)) 
            {
                MessageBox.Show("Cannot find map file: " + mapFileName, "Application error", MessageBoxButtons.OK, MessageBoxIcon.Error);
                return;
            }
            if (mapFileName != null)
            {

                map.Layers.Clear();


                //SHP VECTOR LAYER
                VectorLayer vectorLayerBase = new VectorLayer(mapFileName);
                vectorLayerBase.Style.EnableOutline = true;
                vectorLayerBase.Style.Outline = Pens.Black;
                vectorLayerBase.DataSource = new SharpMap.Data.Providers.ShapeFile(mapFileName, true);
             

                #region Open Street Map
                //TILE LAYER - OPENSTREETMAP
                //MODO FACILE: SENZA CACHE 
                //TileLayer osmTileLayer = new TileLayer(new BruTile.Web.OsmTileSource(), "OSM");
                #region cached on file layer                
                //MODO DIFFICILE: CON CACHE
                FileCache fc = new FileCache(TestMapForm.Default.tileCacheDir, "png");

                //CredentialCache.DefaultCredentials.GetCredential("172.26.254.254:3128", 3128, "basic", new NetworkCredential("fcolledani", "kudu"));
                //TODO: ADD AUTH HERE
                TmsRequest req = new TmsRequest(new Uri("http://b.tile.openstreetmap.org"), "png");
                
                
                WebTileProvider provider = new WebTileProvider(req,fc);
                ITileSource tileSource = new TileSource(provider, new SphericalMercatorInvertedWorldSchema());

                TileLayer osmTileLayer = new TileLayer(tileSource, "OSM");
                #endregion

                map.Layers.Add(osmTileLayer);

                //COORDINATE TRANSFORM

                //The SRS for this datasource is EPSG:4326, therefore we need to transfrom it to OSM projection
                /*    
                 var shpProvider = new SharpMap.Data.Providers.ShapeFile(mapFileName, true);
                 var ctf = new ProjNet.CoordinateSystems.Transformations.CoordinateTransformationFactory();
                 var cf = new ProjNet.CoordinateSystems.CoordinateSystemFactory();
                 var epsgSHP = shpProvider.CoordinateSystem;
                 var epsg3785 = cf.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\",\"3785\"]]");
                 var ct = ctf.CreateFromCoordinateSystems(epsgSHP, epsg3785);
                 vectorLayerBase.CoordinateTransformation = ct;
                 */
                vectorLayerBase.CoordinateTransformation = GetCoordinateTransformation();

                //ALT WAY:

                #endregion
                

                /*
                #region BingMaps
                TileLayer bingTileLayer = new TileLayer(new BingTileSource(BingRequest.UrlBingStaging, "", BingMapType.Roads), "Bing Tiles");
                map.Layers.Add(bingTileLayer);
                #endregion
                */

                /*
                #region Google Maps
                TileLayer
                #endregion*/
                
             
                



                map.Layers.Add(vectorLayerBase);


                //TODO: HERE WE SHALL CREATE THE LAYERS FOR THE PORTALS
                //VectorLayer[] layersArray = new VectorLayer(PortaliDB.count)


                mapImage1.Map = map;
                map.MaximumZoom = TestMapForm.Default.maxZoomLevel;
                map.MinimumZoom = TestMapForm.Default.minZoomLevel;
                SharpMap.Geometries.Point center = new SharpMap.Geometries.Point(TestMapForm.Default.mapCenterLatitude, TestMapForm.Default.mapCenterLongitude);
                map.Center = center;
                mapImage1.AllowDrop = true;
                mapImage1.Map.Zoom = TestMapForm.Default.maxZoomLevel;
                mapImage1.Refresh();



                
            }
        
        }

        /// <summary>
        /// Trasformazione da WGS84 in Mercator
        /// </summary>
        /// <returns></returns>
        private static ICoordinateTransformation GetCoordinateTransformation()
        {

            //The SRS for this datasource is EPSG:4326, therefore we need to transfrom it to OSM projection
            var ctf = new ProjNet.CoordinateSystems.Transformations.CoordinateTransformationFactory();
            var cf = new ProjNet.CoordinateSystems.CoordinateSystemFactory();
            var epsg4326 = cf.CreateFromWkt("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\"]]");
            var epsg3785 = cf.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\",\"3785\"]]");
            return ctf.CreateFromCoordinateSystems(epsg4326, epsg3785);
        }
        #endregion

        private void OnResizeEnd(object sender, EventArgs e)
        {
            mapImage1.Refresh();
        }

        private void selectMapTool(object sender, EventArgs e)
        {
            if (sender == rbtnPan)
            {
                mapImage1.ActiveTool = SharpMap.Forms.MapImage.Tools.Pan;
            }
            else if (sender == rbtnZoomIn)
            {
                mapImage1.ActiveTool = SharpMap.Forms.MapImage.Tools.ZoomIn;
            }
            else if (sender == rbtnZoomOut) 
            {

                mapImage1.ActiveTool = SharpMap.Forms.MapImage.Tools.ZoomOut;
            }
        }

        private void onPickPosition(object sender, EventArgs e)
        {

        }

        private void onMapCenterChanged(SharpMap.Geometries.Point center)
        {
            toolStripStatusLabel1.Text = "MAP CENTER: " + center.AsText();
        }

        private void onMapZoomChanged(double zoom)
        {
            toolStripStatusLabel1.Text = "MAP Zoom: " + zoom.ToString();
        }
    }
}

Coordinator
Feb 23, 2011 at 12:03 PM

Hello FColle,

you have three options:

  • Transform your shapefile with ogr2ogr to another spatial reference system.
  • Customize your epsg3785 wkt string with False_Easting and False_Northing values (try and error)
  • Use DotSpatial.Projections for coordinate transformation (use ReleaseDS configuration). This implies your project must build on .Net40 Framework.

If the shapefile in question is fixed, I'd go for the first.

Hth FObermaier

Feb 23, 2011 at 12:15 PM

Thank you I will try as soon!

Apr 19, 2012 at 6:08 AM

I am very sorry I have to open an old thread...

Is there a way to create a Layer.CoordinateTransformation (ProjNet.CoordinateSystems.Transformations.CoordinateTransformation) that correctly converts EPSG:4326 (WGS-84) coordinates to EPSG:3857 (WGS 84 / Pseudo-Mercator)?

In SharpMap's examples - WinFormSamples example - it shows a way to convert from EPSG:4326 coordinates to EPSG:3857... This transformation produces, results that are not correct; I have seen results that y axis values are off up to 25km!

 

Thanks,
George J.

Developer
Apr 19, 2012 at 10:10 AM

If you take a look at the trunk, inside SharpMap.Demo.Wms project there are some data that is projected correctly.

You can take the projection code from here and use it wherever you want.

Jan 30, 2015 at 11:02 AM
Hi to all,
ran into same Problem as fcolle

Got some shapefiles from customers with projektion = WGS_84_Pseudo_Mercator and a sharpmap project with its base on OSM tiles in Background ...
Doin all this official (EPSG 3587) WKT stings, but always got this 20-50 km out-of-axis-thing.

hmmm run up+down that hill to get the right transformation to work and found this great article :

the-google-maps-bing-maps-spherical-mercator-projection

(damned good job from Alastair ; thumbs up )

so even I now understand the chaos about the different EPSG`s and this tricky projections :-)

So the following transformation will do the job for me :

(sorry vb.net -->C# replace "" with \")
                        Dim ctFact As New ProjNet.CoordinateSystems.Transformations.CoordinateTransformationFactory()

                        Dim wktstring As String = "PROJCS[""Popular Visualisation CRS / Mercator"", GEOGCS[""Popular Visualisation CRS"", DATUM[""WGS 84"", SPHEROID[""WGS 84"", 6378137, 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""]]"

                        Dim PseudoWebMercator As GeoAPI.CoordinateSystems.ICoordinateSystem = cfs.CreateFromWkt(wktstring)

                        cTrans = ctFact.CreateFromCoordinateSystems(PseudoWebMercator, ProjNet.CoordinateSystems.ProjectedCoordinateSystem.WebMercator)
                        CTransReverse = ctFact.CreateFromCoordinateSystems(ProjNet.CoordinateSystems.ProjectedCoordinateSystem.WebMercator, PseudoWebMercator)

                        layer.CoordinateTransformation = cTrans
                        layer.ReverseCoordinateTransformation = CTransReverse
Jan 30, 2015 at 1:49 PM
Hi,

This thread may help you: https://projnet.codeplex.com/discussions/352813


Regards,
George J.