coordinate projection in shapefile

Topics: Algorithms, Data Access, SharpMap Project, SharpMap v0.9 / v1.x
Feb 3, 2011 at 5:18 PM

Hi there!

I have a shapefile that includes the administrative divisions of Portugal. I want to know if a given WGS84 coordinate is within a certain division...

The shapefile comes with a prj file like this:

PROJCS["ETRS_1989_TM06-Portugal",GEOGCS["GCS_ETRS_1989",DATUM["D_ETRS_1989",SPHEROID["GRS_1980",6378137.0,298.257222101]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Transverse_Mercator"],PARAMETER["False_Easting",0.0],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",-8.133108333333333],PARAMETER["Scale_Factor",1.0],PARAMETER["Latitude_Of_Origin",39.66825833333333],UNIT["Meter",1.0]]

 

I am using the following code to get the info from the division where the point is:

 

 

    private string GetStateName(double x, double y)
    {
        double ptX = 0;
        double ptY = 0;

        if (x != 0 || y != 0)
        {
            ptX = x;
            ptY = y;
        }
        string sHASC = "";
        SharpMap.Data.Providers.ShapeFile oShape = new SharpMap.Data.Providers.ShapeFile(Server.MapPath(@"~\gis_data\Cont_AAD_CAOP2010.shp"));
        oShape.Open();
        uint i = 0;
        uint iFeature = (uint)oShape.GetFeatureCount();

        SharpMap.Geometries.Point oPoint = new SharpMap.Geometries.Point(ptX, ptY);
        GisSharpBlog.NetTopologySuite.Geometries.Geometry oPt = SharpMap.Converters.NTS.GeometryConverter.ToNTSGeometry(oPoint, new GisSharpBlog.NetTopologySuite.Geometries.GeometryFactory());

        // find a record
        for (i = 0; i < iFeature; i++)
        {
            GisSharpBlog.NetTopologySuite.Geometries.Geometry oPoly = SharpMap.Converters.NTS.GeometryConverter.ToNTSGeometry(oShape.GetFeature(i).Geometry, new GisSharpBlog.NetTopologySuite.Geometries.GeometryFactory());
            
            if (oPt.Within(oPoly))
            {
                sHASC += oShape.GetFeature(i).ItemArray.GetValue(1).ToString() + " - " + oShape.GetFeature(i).ItemArray.GetValue(2).ToString() + " - " + oShape.GetFeature(i).ItemArray.GetValue(3).ToString() + " |||||||| ";
                //break;
            }
        }
        return sHASC;
    }

I call the function like this: 

 

GetStateName(-9.159333, 38.716833)

But I always get the same division no matter the coordinates I insert... I think it has to do with the projection of the shapefile that is used... but I don't know how to change that! I am trying with the program MapWindows GIS and the reproject utility, but no luck so far! Does anyone have an idea on how to change this? I've tested with another shapefile that has other divisions, and it works fine so it really has something to do with the projection in use...
If someone can help me I would appreciate it!

Thanks!

 

Coordinator
Feb 3, 2011 at 8:04 PM

Hello Hulkman

you are heading the right direction, you either need to transform your geometries to epsg:4326 or your input coordinates to that of the shapefile.

I'd go for the second.

What you need to do is the following:

//Setting up coordinate transformation factory
var ctFac = new CoordinateTransformationFactory();

//Setting up coordinate system fractory
var csFac = new CoordinateSystemFactory();

//Setting up EPSG:4326
var epsg4326 = csFac.CreateFromFromWkt(...);

//Setting up  transformation
var ct = ctFac.CreateFromCoordinateSystems(epsg4326, oShape.CoordinateSystem);

//Perform coordinate transformation
var c = new double[] { x, y };
var cTrans = ct.MathTransform.Transform(c);

var pt = new Point( cTrans );

//Seek for possible features
var fds = new FeatureDataSet();
oShape.ExecuteIntersectionQuery( pt.GetBoundingBox(), fds);

foreach ( var fdr in fds.Tables[0].Rows)
{
//Perform true intersection test with NTS
}

Note: This code has not seen a compiler, there might be typos, but I hope you get the idea

FObermaier

Feb 10, 2011 at 3:01 PM

Hi again!

Ive finally managed to test your sample code and got the same results I had after sucessfully reprojecting the shapefiles... Thanks for the help!