Shapefile polygons not positioned correctly in VirtualEarth

Topics: SharpMap Project, SharpMap v2.0
May 8, 2009 at 9:21 PM
Edited May 8, 2009 at 10:51 PM

Hi, I have some Shapefiles that I've converted from NAD83 to WGS84 and loaded into SQL Server 08. I'm able to view them in Virtual Earth via my web app. However, they are positioned about 10 feet too far to the north. Is there any way I can nudge them all 10 feet to the south when I convert my coordinates from one spatial reference system to another?

Thanks,

Brad

Coordinator
May 9, 2009 at 9:11 PM

Hi Brad I think VE uses EPSG:900913 hth jd

May 11, 2009 at 6:15 PM

No, I tried that but it didn't work. Actually they're closer to 40 feet off, too far north.

This is my code. I boldfaced where I changed the SRID values.  I've got a coordinate transformation assigned to the ShapeFile that parses out the WKT for WGS84, but I didn't replace it with a 900913 WKT because I didn't think it was being used anyway. Do you think I missed anything?

    private void LoadShapeFileToSQL(string filePath)
    {
        InitSRIDMap();
        SqlConnection conn = new SqlConnection();
        conn.ConnectionString = @"Data Source=....";
        SqlCommand cmd = new SqlCommand("Insert into [QA_GIS]..Parcel (MBA,ParcelNo,SeqNo,Geom,Area,Perimeter,TaxYear,dtCreate) Values (@MBA,@ParcelNo,@SeqNo,@Geometry,@Area,@Perimeter,@TaxYear,@dtCreate)", conn);
        cmd.Parameters.Add("@MBA", SqlDbType.Decimal);
        cmd.Parameters.Add("@Area", SqlDbType.Real);
        cmd.Parameters.Add("@Perimeter", SqlDbType.Real);
        cmd.Parameters.Add("@ParcelNo", SqlDbType.VarChar);
        cmd.Parameters.Add("@SeqNo", SqlDbType.TinyInt);
        cmd.Parameters.Add("@TaxYear", SqlDbType.Int);
        cmd.Parameters.Add("@dtCreate", SqlDbType.DateTime);
        SqlParameter sp = new SqlParameter();
        sp.ParameterName = "@Geometry";
        sp.UdtTypeName = "Geometry";
        cmd.Parameters.Add(sp);
        Random random = new Random((int)DateTime.Now.Ticks);
        // SharpMap's Shapefile reader / writer
        ShapeFileProvider sf = null;
        GeometryServices gstest = new GeometryServices();
        GeometryServices gs = new GeometryServices();
        sf = new ShapeFileProvider(filePath + ".shp", gs.DefaultGeometryFactory, gs.CoordinateSystemFactory);
        //this keeps it from loading the entire spatial index
        sf.IsSpatiallyIndexed = false;
        sf.Open(true);

        //ICoordinateSystem srs = gs.CoordinateSystemFactory.CreateFromWkt("PROJCS[\"NAD_1983_StatePlane_Texas_North_Central_FIPS_4202_Feet\",GEOGCS[\"GCS_North_American_1983\",DATUM[\"D_North_American_1983\",SPHEROID[\"GRS_1980\",6378137.0,298.257222101]],PRIMEM[\"Greenwich\",0.0],UNIT[\"Degree\",0.0174532925199433]],PROJECTION[\"Lambert_Conformal_Conic\"],PARAMETER[\"False_Easting\",1968500.0],PARAMETER[\"False_Northing\",6561666.666666666],PARAMETER[\"Central_Meridian\",-98.5],PARAMETER[\"Standard_Parallel_1\",32.13333333333333],PARAMETER[\"Standard_Parallel_2\",33.96666666666667],PARAMETER[\"Latitude_Of_Origin\",31.66666666666667],UNIT[\"Foot_US\",0.3048006096012192]]");
        StreamReader sr = new StreamReader(filePath + ".prj");
        String wkt = sr.ReadLine();
        ICoordinateSystem srs = gs.CoordinateSystemFactory.CreateFromWkt(wkt);
        ICoordinateSystem wgs84 = gs.CoordinateSystemFactory.CreateFromWkt("GEOGCS[\"GCS_WGS_1984\",DATUM[\"D_WGS_1984\",SPHEROID[\"WGS_1984\",6378137.0,298.257223563]],PRIMEM[\"Greenwich\",0.0],UNIT[\"Degree\",0.0174532925199433]]");
        ICoordinateTransformation ics =
            gs.CoordinateTransformationFactory.CreateFromCoordinateSystems(srs, wgs84);
        sf.CoordinateTransformation = ics;
        //IGeometryFactory gfVE = gs["EPSG:4326"];
        IGeometryFactory gfVE = gs["EPSG:900913"];

        if (sf.ShapeType != ShapeType.Polygon)
            throw new Exception();
        //IFeatureDataReader ifdr = sf.GetReader();

        IGeometry geom;
        ICoordinate2D i2d;

        conn.Open();
        String prcl;
        String lastprcl = "";
        int seq = 1;
        int errorcnt = 0;
        int linenumber = 1;
        bool morerecs = true;
        IFeatureDataRecord ifdr;
        uint a;
        int shapesloaded = 0;
        int rem;
        int quot;
        for (a = 1; a <= sf.FeatureCount;a++ )
        {
            try
            {
                ifdr = sf.GetFeatureByOid(a);
               
                if (!String.IsNullOrEmpty(ifdr.GetString(ifdr.GetOrdinal("Geo")+1)))
                {
                    try
                    {
                        SqlGeometryBuilder gb = new SqlGeometryBuilder();
                        gb.SetSrid(900913);
                        //gb.SetSrid(4326);
                        geom = ifdr.Geometry;
                        SqlBytes b = new SqlBytes(geom.AsBinary());
                        //SqlGeometry g = SqlGeometry.STGeomFromWKB(b, 4326);
                        SqlGeometry g = SqlGeometry.STGeomFromWKB(b, 900913);
                        cmd.Parameters["@MBA"].Value = "4825100000";
                        cmd.Parameters["@SeqNo"].Value = seq;
                        cmd.Parameters["@Area"].Value = ifdr.GetFloat(ifdr.GetOrdinal("SHAPE_Area") + 1);
                        cmd.Parameters["@Perimeter"].Value = ifdr.GetFloat(ifdr.GetOrdinal("SHAPE_Leng") + 1);
                        cmd.Parameters["@ParcelNo"].Value = ifdr.GetString(ifdr.GetOrdinal("Geo") + 1);
                        cmd.Parameters["@TaxYear"].Value = 2008;
                        cmd.Parameters["@dtCreate"].Value = System.DateTime.Now;
                        sp.Value = g;
                        cmd.ExecuteNonQuery();
                        seq = 1;
                        shapesloaded++;
                        quot = (int)(Math.DivRem((int)a, 100, out rem));
                        if (rem == 0) { Console.WriteLine(a.ToString()); }
                    }
                    catch (SqlException e) { Console.WriteLine("Error " + e.Message + " at ID:" + a.ToString()); errorcnt++; seq++; }
                    catch (Exception e) { Console.WriteLine("Error " + e.Message + " at ID:" + a.ToString()); errorcnt++; }
                }
            }
            catch (Exception e) { a++; Console.WriteLine("Error " + e.Message + " at ID:" + a.ToString()); errorcnt++; }
            linenumber++;
        }

        conn.Close();

        Console.WriteLine("Total Features in ShapeFile:" + sf.FeatureCount.ToString());
        Console.WriteLine(shapesloaded.ToString() + " Shapes loaded");
        sf.Close();
        Console.WriteLine("Total Errors:" + errorcnt.ToString());
    }

 

May 11, 2009 at 6:23 PM

BTW, I found this WKT for the 900913 projection system. Do you think I should use it?

900913=PROJCS["WGS84 / Simple Mercator",GEOGCS["WGS84", DATUM["WGS_1984",SPHEROID["WGS_1984", 6378137.0, 298.257223563]],PRIMEM["Greenwich", 0.0],UNIT["degree", 0.017453292519943295],AXIS["Longitude", EAST],AXIS["Latitude", NORTH]],PROJECTION["Mercator_1SP_Google"],PARAMETER["latitude_of_origin", 0.0],PARAMETER["central_meridian", 0.0],PARAMETER["scale_factor", 1.0],PARAMETER["false_easting", 0.0],PARAMETER["false_northing", 0.0], UNIT["m", 1.0],AXIS["x", EAST], AXIS["y",NORTH], AUTHORITY["EPSG","900913"]]

Coordinator
May 11, 2009 at 6:41 PM

I was going to suggest almost the same : http://spatialreference.org/ref/sr-org/6627/ogcwkt/ It is possible that the one included in proj4utilities is incorrect

then

create a coordinate system from the shapefiles prj file

create another one from the wkt above

create a transformation object between the two coordinate systems

try transforming the geometries by calling transform.Tranform(geometry, gs.DefaultGeometryFactory)

 

hth jd

May 11, 2009 at 11:53 PM

It wasn't able to parse that WKT string. It bombed in GeoAPI...WktTokenizer because of "expected ',' found AXIS at..."

Coordinator
May 12, 2009 at 11:17 AM

the parser doesn't like spaces after the comma (,). remove all the spaces outside of quotes and you should be all set. Unfortunately the WKT property of the crs suggests using a space after the comma.

Hth

FObermaier

May 12, 2009 at 11:25 PM

I removed all the embedded spaces but I'm still getting this error. Does the WKT look ok?

"Expecting (',') but got a 'AXIS' at line 1 column 200."

Column 200 is after the word "AXIS" in the bold area.

PROJCS["WGS84 / Simple Mercator",GEOGCS["WGS84",DATUM["WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["degree",0.017453292519943295],AXIS["Longitude",EAST],AXIS["Latitude",NORTH]],PROJECTION["Mercator_1SP_Google"],PARAMETER["latitude_of_origin",0.0],PARAMETER["central_meridian",0.0],PARAMETER["scale_factor",1.0],PARAMETER["false_easting",0.0],PARAMETER["false_northing",0.0],UNIT["m",1.0],AXIS["x",EAST],AXIS["y",NORTH],AUTHORITY["EPSG","900913"]]

   at GeoAPI.IO.WellKnownText.WktTokenizer.ReadToken(String expectedToken) in C:\Code\VisualStudio\projects\sharpmapv2\GeoAPI\src\GeoAPI\IO\WellKnownText\WktTokenizer.cs:line 81
   at GeoAPI.IO.WellKnownText.CoordinateSystemWktReader.readProjectedCoordinateSystem(WktTokenizer tokenizer, ICoordinateSystemFactory factory) in C:\Code\VisualStudio\projects\sharpmapv2\GeoAPI\src\GeoAPI\IO\WellKnownText\CoordinateSystemWktReader.cs:line 398
   at GeoAPI.IO.WellKnownText.CoordinateSystemWktReader.readCoordinateSystem(WktTokenizer tokenizer, ICoordinateSystemFactory factory) in C:\Code\VisualStudio\projects\sharpmapv2\GeoAPI\src\GeoAPI\IO\WellKnownText\CoordinateSystemWktReader.cs:line 232
   at GeoAPI.IO.WellKnownText.CoordinateSystemWktReader.Parse(TextReader reader, ICoordinateSystemFactory factory) in C:\Code\VisualStudio\projects\sharpmapv2\GeoAPI\src\GeoAPI\IO\WellKnownText\CoordinateSystemWktReader.cs:line 97
   at GeoAPI.IO.WellKnownText.CoordinateSystemWktReader.Parse(String wkt, ICoordinateSystemFactory factory) in C:\Code\VisualStudio\projects\sharpmapv2\GeoAPI\src\GeoAPI\IO\WellKnownText\CoordinateSystemWktReader.cs:line 61
   at GeoAPI.IO.WellKnownText.WktReader`1.ToCoordinateSystemInfo(String wktData, ICoordinateSystemFactory`1 factory) in C:\Code\VisualStudio\projects\sharpmapv2\GeoAPI\src\GeoAPI\IO\WellKnownText\WktReader.cs:line 55
   at ProjNet.CoordinateSystems.CoordinateSystemFactory`1.CreateFromWkt(String Wkt) in C:\Code\VisualStudio\projects\sharpmapv2\Proj.Net\ProjNet\CoordinateSystems\CoordinateSystemFactory`1.cs:line 95
   at ProjNet.CoordinateSystems.CoordinateSystemFactory`1.GeoAPI.CoordinateSystems.ICoordinateSystemFactory.CreateFromWkt(String wkt) in C:\Code\VisualStudio\projects\sharpmapv2\Proj.Net\ProjNet\CoordinateSystems\CoordinateSystemFactory`1.cs:line 868
   at SQLGIS.LoadShapeFileToSQL(String filePath) in C:\Code\VisualStudio\projects\LoadShapeFilesToSQL08\SQLGIS.cs:line 88

Coordinator
May 13, 2009 at 11:29 AM

I submitted a patchfile solving your issue with AXIS definition in WKT. See:

http://sharpmap.codeplex.com/SourceControl/PatchList.aspx ID #2879

Hth

FObermaier

May 13, 2009 at 3:02 PM

What is the best way to merge this patch into the source? Should I take a Subversion update of the CoordinateSystemWktReader.cs class, or copy and paste the changes? Also, there are three CoordinateSystemWktReader.cs files, one in ProjNet, one in GeoAPI, and one in NTS. Do they all need to be updated?

May 13, 2009 at 3:27 PM

OK, I figured that last one out. I downloaded the patch file into the GEOAPI folder, right-clicked on it and selected Tortoise-->Apply Patch. Then I rebuilt the whole v2.0 solution, updated the references to said solution in my project, and re-ran.

Now I'm getting this error:

{"Projection Mercator_1SP_Google is not supported."}   Something tells me this isn't good.

"   at ProjNet.CoordinateSystems.Transformations.CoordinateTransformationFactory`1.createCoordinateOperation(IProjection projection, IEllipsoid ellipsoid, ILinearUnit unit) in C:\\Code\\VisualStudio\\projects\\sharpmapv2\\Proj.Net\\ProjNet\\CoordinateSystems\\Transformations\\CoordinateTransformationFactory`1.cs:line 597\r\n   at ProjNet.CoordinateSystems.Transformations.CoordinateTransformationFactory`1.geographicToProjected(IGeographicCoordinateSystem`1 source, IProjectedCoordinateSystem`1 target) in C:\\Code\\VisualStudio\\projects\\sharpmapv2\\Proj.Net\\ProjNet\\CoordinateSystems\\Transformations\\CoordinateTransformationFactory`1.cs:line 286\r\n   at ProjNet.CoordinateSystems.Transformations.CoordinateTransformationFactory`1.CreateFromCoordinateSystems(ICoordinateSystem`1 source, ICoordinateSystem`1 target) in C:\\Code\\VisualStudio\\projects\\sharpmapv2\\Proj.Net\\ProjNet\\CoordinateSystems\\Transformations\\CoordinateTransformationFactory`1.cs:line 108\r\n   at ProjNet.CoordinateSystems.Transformations.CoordinateTransformationFactory`1.projectedToProjected(IProjectedCoordinateSystem`1 source, IProjectedCoordinateSystem`1 target) in C:\\Code\\VisualStudio\\projects\\sharpmapv2\\Proj.Net\\ProjNet\\CoordinateSystems\\Transformations\\CoordinateTransformationFactory`1.cs:line 250\r\n   at ProjNet.CoordinateSystems.Transformations.CoordinateTransformationFactory`1.CreateFromCoordinateSystems(ICoordinateSystem`1 source, ICoordinateSystem`1 target) in C:\\Code\\VisualStudio\\projects\\sharpmapv2\\Proj.Net\\ProjNet\\CoordinateSystems\\Transformations\\CoordinateTransformationFactory`1.cs:line 179\r\n   at ProjNet.CoordinateSystems.Transformations.CoordinateTransformationFactory`1.GeoAPI.CoordinateSystems.Transformations.ICoordinateTransformationFactory.CreateFromCoordinateSystems(ICoordinateSystem source, ICoordinateSystem target) in C:\\Code\\VisualStudio\\projects\\sharpmapv2\\Proj.Net\\ProjNet\\CoordinateSystems\\Transformations\\CoordinateTransformationFactory`1.cs:line 80\r\n   at SQLGIS.LoadShapeFileToSQL(String filePath) in C:\\Code\\VisualStudio\\projects\\LoadShapeFilesToSQL08\\SQLGIS.cs:line 89"

May 13, 2009 at 5:02 PM

I'm pretty new to GIS but I'm wondering what could cause my shapes to be positioned 40 feet too far north. Could the false northing value be incorrect?

 

May 13, 2009 at 8:43 PM

I loaded another shapefile from a different county to see if the shapes were positioned incorrectly on that one, too. I converted it from NAD83 to WGS84 using the same parameters as the other file. The shapes on that one are also shifted about 40 feet to the north. Both counties use the same projection system, but they are completely different entities.

May 15, 2009 at 4:18 PM

I'm still having this problem. All my shapes are still a good 40 or 50 feet off. Any other suggestions, please?

Thank you,

Brad

Coordinator
May 18, 2009 at 7:55 AM

If you don't get any further with converting the SRS, you could try OGR2OGR of the FWTools collection to convert your shapefiles to the desired SRS and then importing the modified shapefile with SharpMap omitting the coordinate transformation.

In your case it could be:

ogr2ogr -f "ESRI Shapefile" original.shp 900913.shp -s_srs EPSG:2276 -t_srs EPSG:900913

Instead of the EPSG-codes you can specify filenames containing the wkt representation of the desired spatial reference systems.

Hth

FObermaier

May 19, 2009 at 4:00 PM

Thx but I can't get that thing to work. It says it can't create the output file.

I have a question. Are the equations used to convert between SRS's standard, or did the SharpMap project write its own? Also, is a 40 to 50 foot margin of error normal when converting from a regional SRS to WGS84?

I noticed something odd when I did the CoordinateTransformation last time. I'll post something here shortly to see what you think.

Thank you

Coordinator
May 20, 2009 at 11:22 AM

Hi Brad, the equations are standard. 40-50 foot sounds like quite a big error to me. http://projnet.codeplex.com/Thread/View.aspx?ThreadId=45600 may help. hth jd

May 20, 2009 at 4:30 PM

Thx John, I posted my question to that thread.

http://projnet.codeplex.com/Thread/View.aspx?ThreadId=45600

Brad