Rendering .shp from wgs84 to mercator spheric

Topics: Algorithms, SharpMap v0.9 / v1.x
Jan 14, 2008 at 7:04 PM
Hello,

I have street .shp data for a city, whose coordinates are in (the normal, ellipsoidal) WGS84. I've been trying to render the data in a standalone application using the same projection that Google Maps uses. I am utilizing Sharpmap v0.9 (I am using the Trunk/ code from changelist 26076).

I've read multiple threads and read through SharpGIS' July 2007 entry regarding the "spheric mercator" projection, but am confused about how to map the lat/lng coordinates from the data to the rendering widget. I have two questions:

1) Is it enough to just use the Mercator (mercator_2sp) projection, or do I need to use the long WKT SharpGIS cites in order to properly map the (ellipsoidal) WGS84 data to the spheric mercator projection?

2) Along those lines, I tried to create the transformation using the "spheric mercator" WKT string, but the transformation was not rendering anything. I explain what I did below. Could someone provide some assistance here, as to why the code is not working?

I did some testing with the DemoWinForm. I modified the loadLayer() code to apply a coordinate transform to the new layer, ran the program, loaded the .shp data, and then zoomed to extents. The changed code looked like this:

private void loadLayer()
    {
      DialogResult result = AddLayerDialog.ShowDialog( this );
 
      if ( result == DialogResult.OK )
      {
        foreach ( string fileName in AddLayerDialog.FileNames )
        {
          string extension = Path.GetExtension( fileName );
          ILayerFactory layerFactory = null;
 
          if ( !_layerFactoryCatalog.TryGetValue( extension, out layerFactory ) )
            continue;
 
          /** changed it from ILayer to VectorLayer to access transform */
          VectorLayer layer = layerFactory.Create( Path.GetFileNameWithoutExtension( fileName ), fileName ) as VectorLayer;
 
          addLayer( layer );
 
          /** additions to the original code begin here **/
          CoordinateSystemFactory cFac = new SharpMap.CoordinateSystems.CoordinateSystemFactory();
          ICoordinateSystem cs = ( layer.DataSource as ShapeFile ).CoordinateSystem;
 
          /** SPECIFY TRANSFORM **/
          List<ProjectionParameter> parameters = new List<ProjectionParameter>();
          parameters.Add( new ProjectionParameter( "latitude_of_origin", 0 ) );
          parameters.Add( new ProjectionParameter( "central_meridian", 0 ) );
          parameters.Add( new ProjectionParameter( "false_easting", 0 ) );
          parameters.Add( new ProjectionParameter( "false_northing", 0 ) );
          IProjection projection = cFac.CreateProjection( "Mercator", "Mercator_2SP", parameters );
          
          IProjectedCoordinateSystem coordsys = cFac.CreateProjectedCoordinateSystem( "Mercator",
            cs as IGeographicCoordinateSystem,
            projection, LinearUnit.Metre, new AxisInfo( "East", AxisOrientationEnum.East ), new AxisInfo( "North", AxisOrientationEnum.North ) );
          /** END of TRANSFORM **/
 
          layer.CoordinateTransformation = new CoordinateTransformationFactory().CreateFromCoordinateSystems( cs, coordsys );
        }
 
        changeUIOnLayerSelectionChange();
 
        MainMapImage.Refresh();
      }
    }

When I zoomed to extents, the map's envelope was set to:

{-8259763.31220257,4937989.70537966 -8209016.92735487,4966387.17730453}

and when the transformation was applied to the envelope, the resulting bounding box turned to:

{-74.1987162656125,40.6840500000287 -73.7428537343875,40.8779600000284}

which makes sense. (Although, again, I do not know if there should be another transform applied to convert the data's WGS84 coords to the spheric coords.)

I then tried to replace the area marked /** SPECIFY TRANSFORM **/ ... /** END of TRANSFORM **/ with this:

          ICoordinateSystem coordsys = cFac.CreateFromWkt(
            @"PROJCS[""Mercator Spheric"", GEOGCS[""WGS84basedSpheric_GCS""," +
            @"DATUM[""WGS84basedSpheric_Datum"", SPHEROID[""WGS84based_Sphere"", 6378137, 0], TOWGS84[0, 0, 0, 0, 0, 0, 0]]," +
            @"PRIMEM[""Greenwich"", 0, AUTHORITY[""EPSG"", ""8901""]]," +
            @"UNIT[""degree"", 0.0174532925199433, AUTHORITY[""EPSG"", ""9102""]]," +
            @"AXIS[""E"", EAST], AXIS[""N"", NORTH]]," +
            @"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]]" );

When I did this, I got an exception when I tried to load the .shp data as a layer. The error was in Point's constructor; it stated that "Only 2 dimensions are supported for points". So is the WKT string wrong, since it's creating a 3D point?

I decided to comment the sanity check out, since it seemed like the third coordinate would always be 0. When I did "zoom to extents", the map's envelope got set to:

{-0.000666057768367583,0.000367934130223892 -0.000662923966508798,0.000369687793082242}

and the transformed bounding box:

{-5.9832987343114E-09,3.32748241226463E-09 -5.95514731323991E-09,3.34333866684081E-09}

Any assistance would be appreciated.

Thanks,
Jon
Nov 14, 2008 at 3:36 PM
Jon,

Have you had any success converting shapes from WGS84 to Spherical Mercator?

Thanks,
Corey
Nov 15, 2008 at 8:16 AM
Edited Nov 15, 2008 at 8:27 AM
Corey,

I am successfully using the following code:

 

public static ICoordinateTransformation wgs84toGoogle()
{           
  CoordinateSystemFactory csFac = new SharpMap.CoordinateSystems.CoordinateSystemFactory();
  CoordinateTransformationFactory ctFac = new CoordinateTransformationFactory();

  IGeographicCoordinateSystem wgs84 = csFac.CreateGeographicCoordinateSystem(
    "WGS 84", AngularUnit.Degrees, HorizontalDatum.WGS84, PrimeMeridian.Greenwich,
    new AxisInfo("north", AxisOrientationEnum.North), new AxisInfo("east", AxisOrientationEnum.East)
  );

  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("Google Mercator", "mercator_1sp", parameters);

  IProjectedCoordinateSystem epsg900913 = csFac.CreateProjectedCoordinateSystem(

    "Google Mercator", wgs84, projection, LinearUnit.Metre,
new AxisInfo("East", AxisOrientationEnum.East), new AxisInfo("North", AxisOrientationEnum.North));

  return ctFac.CreateFromCoordinateSystems(wgs84, epsg900913);
}

It gives me a perfect match with Google Maps.

Hope this helps,
Marcin Bulandra
Nov 18, 2008 at 2:14 AM
Thanks Marcin,

My shapes now display better on my map.

Thanks,
Corey


Nov 18, 2008 at 9:24 AM
Hi. 

I have similar woes.  I am very new to GIS, Sharpmap et al, and have to develop a standalone app to render SHP files with coordinate system
shp.CoordinateSystem == GEOGCS["GCS_WGS_1984", DATUM["D_WGS_1984", SPHEROID["WGS_1984", 6378137, 298.257223563]], PRIMEM["Greenwich", 0], UNIT["Degree", 0.0174532925199433]]
using the UTM, how would i go about this.

however, when i try and display the map, all i get is a great blank space.

I have tried to tell the Map object to ZoomToExtents(), but still, nothing.

Here is the bits of the code that is not doing what i thought it would be doing:

<code>
        private IProjectedCoordinateSystem CreateUtmProjection(int utmZone)
        {
            CoordinateSystemFactory cFac = new SharpMap.CoordinateSystems.CoordinateSystemFactory();

            //Create geographic coordinate system based on the WGS84 datum
            IEllipsoid ellipsoid = cFac.CreateFlattenedSphere("WGS 84", 6378137, 298.257223563, LinearUnit.Metre);
            IHorizontalDatum datum = cFac.CreateHorizontalDatum("WGS_1984", DatumType.HD_Geocentric, ellipsoid, null);
            IGeographicCoordinateSystem gcs = cFac.CreateGeographicCoordinateSystem("WGS 84", AngularUnit.Degrees, datum,
              PrimeMeridian.Greenwich, new AxisInfo("Lon", AxisOrientationEnum.East),
              new AxisInfo("Lat", AxisOrientationEnum.North));

            //Create UTM projection
            List<ProjectionParameter> parameters = new List<ProjectionParameter>();
            parameters.Add(new ProjectionParameter("latitude_of_origin", 0));
            parameters.Add(new ProjectionParameter("central_meridian", -183 + 6 * utmZone));
            parameters.Add(new ProjectionParameter("scale_factor", 0.9996));
            parameters.Add(new ProjectionParameter("false_easting", 500000));
            parameters.Add(new ProjectionParameter("false_northing", 0.0));
            IProjection projection = cFac.CreateProjection("transverse_mercator", "transverse_mercator", parameters);

            return cFac.CreateProjectedCoordinateSystem("WGS 84 / UTM zone " + utmZone.ToString() + "N", gcs,
               projection, LinearUnit.Metre, new AxisInfo("East", AxisOrientationEnum.East),
               new AxisInfo("North", AxisOrientationEnum.North));
        }

public Add(List<ILayer> layers,ShapeFile shp)
{
    newlayer.LayerName =  Path.GetFileNameWithoutExtension(shp.Filename);
            newlayer.Style.Line = new System.Drawing.Pen(Color.White);
            newlayer.Style.Outline = new System.Drawing.Pen(Color.White);
            newlayer.Style.Fill = System.Drawing.Brushes.LightYellow;

            //Create zone UTM 32N projection
            IProjectedCoordinateSystem utmProj = CreateUtmProjection(32);

           //Create transformation
           CoordinateTransformationFactory ctFac = new CoordinateTransformationFactory();
            ICoordinateTransformation transform = ctFac.CreateFromCoordinateSystems(shp.CoordinateSystem, utmProj);

            //Apply transformation to a vectorlayer
             newlayer.CoordinateTransformation = transform;
             layers.Add(newlayer);
}
<code>
...
And then
...
<code>
Map _Map = new Map(new Size(width, height));

List<ILayer> layers = new List<ILayer>();
Add(layers,SHP1);
Add(layers, SHP2);
         ...
_Map.Layers.AddRange(layers);
_Map.ZoomToExtends();
Img = _Map.GetMap();
</code>
If i remove the coordinate transformation, my map displays perfectly, although the projection is wrong.

If anyone can help me, or even point me in the right direction i would appreciate it very much!

Regards

Erick

Nov 20, 2008 at 5:14 PM
Hi Erick,

While writing the above code I had some similiar problems. First the coordinate transform was wrong and the effect of this was map.GetExetens() returning four NaN values. The other time I was setting layer.MaxVisible in degrees and not in meters. Example: If you set MaxVisible to 360, then in wgs84 coordinate system the layer is visible up to whole-world view, but in projected coordinate system like ours it is visible only if the width of the rendered map is less than 360 meters. I had to multiply MaxVisible in all my layers by 6378137 * 2 * pi() or roughly 40 000 000 meters.

Regards,

Marcin Bulandra
Jul 22, 2009 at 8:58 AM
Edited Jul 22, 2009 at 8:59 AM

Hi mbu, I met accross the same kind of questions like yours. I found that because of the zoom and transformation , part of the map is distorted very much .Besides, some part of it is invisible.

I wanders what you did about this. I am currently want to change the zoom and MaxVisible to try again, wish me good luck.

Dec 15, 2009 at 3:55 AM

 

Hi mbu,
I am using Sharpmap v0.9.
I am rendering an image with the following data PROJCS["British_National_Grid",GEOGCS["GCS_OSGB_1936",DATUM["D_OSGB_1936",SPHEROID["Airy_1830",6377563.396,299.3249646]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Transverse_Mercator"],PARAMETER["False_Easting",400000.0],PARAMETER["False_Northing",-100000.0],PARAMETER["Central_Meridian",-2.0],PARAMETER["Scale_Factor",0.999601272],PARAMETER["Latitude_Of_Origin",49.0],UNIT["Meter",1.0]]
I have tried to use the wgs84toGoogle() method to get a perfect overlap of image on Google Maps. I got NaN for my center, envelope.
Any suggestions on how to recity this.
Thank You
Roche

Hi mbu,

I am using Sharpmap v0.9.

I am rendering an image with the following data PROJCS["British_National_Grid",GEOGCS["GCS_OSGB_1936",DATUM["D_OSGB_1936",SPHEROID["Airy_1830",6377563.396,299.3249646]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Transverse_Mercator"],PARAMETER["False_Easting",400000.0],PARAMETER["False_Northing",-100000.0],PARAMETER["Central_Meridian",-2.0],PARAMETER["Scale_Factor",0.999601272],PARAMETER["Latitude_Of_Origin",49.0],UNIT["Meter",1.0]]

I have tried to use the wgs84toGoogle() method to get a perfect overlap of image on Google Maps. I got NaN for my center, envelope.

Any suggestions on how to recify this.

Thank You

Roche

 

Dec 30, 2009 at 9:39 AM

Hi all

I am experiencing exactly the same issue that Roche describes above, does anyone have any suggestions on this point?

Regards

Mark 

Mar 27, 2010 at 9:55 AM
Edited Mar 27, 2010 at 9:59 AM
coreycoto wrote:
Thanks Marcin,

My shapes now display better on my map.

Thanks,
Corey


Sorry but i've tryed the same code...   it would be the projection returned by this code that is wrong!?!? :

 

public static ICoordinateTransformation wgs84toGoogle()
{           
  CoordinateSystemFactory csFac = new SharpMap.CoordinateSystems.CoordinateSystemFactory();
  CoordinateTransformationFactory ctFac = new CoordinateTransformationFactory();

  IGeographicCoordinateSystem wgs84 = csFac.CreateGeographicCoordinateSystem(
    "WGS 84", AngularUnit.Degrees, HorizontalDatum.WGS84, PrimeMeridian.Greenwich,
    new AxisInfo("north", AxisOrientationEnum.North), new AxisInfo("east", AxisOrientationEnum.East)
  );

  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("Google Mercator", "mercator_1sp", parameters);

  IProjectedCoordinateSystem epsg900913 = csFac.CreateProjectedCoordinateSystem(

    "Google Mercator", wgs84, projection, LinearUnit.Metre,
new AxisInfo("East", AxisOrientationEnum.East), new AxisInfo("North", AxisOrientationEnum.North));

  return ctFac.CreateFromCoordinateSystems(wgs84, epsg900913);
}