Problem: Disparity in distance values (they vary with latitude)

MapWinGIS
1

Hi

Firstly, I am new to using Mapwingis.ocx, so I apologise if I have overlooked something simple to resolve the problem outlined below. I have created a GIS application using Visual Basic (in Visual Studio 2010) – so far so good; all is functioning well. However, there is one problem I have not been able to resolve, so I am looking for some help to find a solution.

The problem relates to the distance values that are generated for polygons when interacting with the map. I am getting a disparity in distance values with latitude. For example, if I programmatically create a 1 hectare polygon (100m x 100m) at the equator, it is created fine; but if I create it at 51.7 degrees latitude (in the UK), the polygon created is 62m x 62m. A similar disparity in the distance occurs if I create a buffer around a shape (setting 100m buffer results in a 62m buffer on the map). Hence, I would like to know if there is something that I have overlooked that is causing this problem. I have tried changing the projection, geoprojection, mapunits, changing the tile provider (or having no tiles), and/or using the Utils ConvertDistance function, but none have worked; and I am struggling to find any further potential solutions online.

The settings I have in the application are:

  • AxMap1.Projection = MapWinGIS.tkMapProjection.PROJECTION_GOOGLE_MERCATOR
  • AxMap1.GeoProjection.SetGoogleMercator()
  • AxMap1.TileProvider = MapWinGIS.tkTileProvider.BingMaps
  • AxMap1…MapUnits = tkUnitsOfMeasure.umMeters

Any help to point me in the right direction would be much appreciated. Happy to provide any further information to help resolve this.

Thanks

John

2

Hello @JohnTz , and welcome.

Could you post the code you are using to create the polygon? That might help to see what’s happening.

Regards,
Jerry.

3

Hi Jerry

Thanks for the quick reply.

The steps/code for generating the polygon are (in summary):

Step 1. Set the map:

AxMap1.Projection = MapWinGIS.tkMapProjection.PROJECTION_GOOGLE_MERCATOR
AxMap1.GeoProjection.SetGoogleMercator()
AxMap1.TileProvider = MapWinGIS.tkTileProvider.BingMaps
AxMap1.MapUnits = tkUnitsOfMeasure.umMeters

Step 2: Create and save the shapefile:

sfile = save_shp_dialog(1, "") 'gets the filename to save the shapefile
Dim sf As New MapWinGIS.Shapefile()
sf.Projection = "+proj=merc +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs"
sf.CreateNewWithShapeID(sfile, MapWinGIS.ShpfileType.SHP_POLYGON)
sf.Save()
ih = AxMap1.NumLayers - 1
hndl = AxMap1.get_LayerHandle(ih)
current_layer = hndl 'Note: current_layer is the public variable that stores the layer currently selected.

Step 3: Create 100m x 100m square polygon:

'Note x and y derived from AxMap mousedown event (AxMap1.PixelToProj(e.x, e.y, x, y) which are then sent to the following routine:
Dim sf As Shapefile = AxMap1.get_Shapefile(current_layer)
Dim shp = New Shape
sf.StartEditingShapes(True)
ext = AxMap1.MaxExtents
xStep = 100
yStep = 100
shp.Create(ShpfileType.SHP_POLYGON)
shp.AddPoint(x, y)
shp.AddPoint(x + xStep, y)
shp.AddPoint(x + xStep, y - yStep)
shp.AddPoint(x, y - yStep)
shp.AddPoint(x, y)
shapeIndex = sf.EditAddShape(shp)
AxMap1.Redraw()

As described previously, this results in a square polygon ok, but the size varies with latitude, e.g. 62m at 51.7 degrees latitude; 100m at the equator.

If you think there is an issue in my code, then do please flag it up. My suspicion is that this is a projection issue; I came across the following articles over the weekend:

I have not been able to find a solution either in terms of a different projection or a mathematical conversion; albeit limitations in my GIS skills and knowledge of projections may also be contributing.

Any help or advice would be much appreciated. Please let me know if you need any further details.

Kind regards

John

4

Hi John, Sorry for the late reply.

I also think the problem is in the projection.
Try to set the projection of the shapefile to match the projection of the map:

// Pseudo code:
sf.GeoProjection = Axmap1.GeoProjection.Clone()
// Or
sf.GeoProjection.SetGoogleMercator()

Projection is deprecated, use Geoprojection instead.

5

Hi Paul,

Thanks for this suggestion. I have tried it, but unfortunately it did not resolve the problem.

For the moment, I have been using a fudge factor based on the latitude, which is far from ideal, but it’s just about ok for the current purpose of the application we are building. That said, I would of course prefer to find the root cause of this problem, so that our solution going forward is more robust. My suspicion is that this is some fundamental incompatibility in the projection used by Google/Bing, but I can’t find any solid evidence of that, it is just a suspicion. Any further thoughts would be very welcome.

Thanks

John

6

Hi John,

when looking at your problem, my first question is: How do you measure the not fitting sizes/distances?
Are you just using the measuring tool or do you query your shape-Object?

BTW:
Google Mercator Projection is not at all made, for what you’re doing here. EPSG 3857 tries to achieve the impossible: A carthesic coordinate system for the whole globe, - made ONLY to make it possible to DISPLAY square tiles all over the world. It was never made for accuracy: If you measure things on such a map in small scale, results may be not that bad, - but don’t dig holes to lay pipes with these -, but the larger the scale growes, the disproportional larger the miss-measurments will get,
and i wouldn’t rule out, that this can be one of the causes of your problem.
So whereever the problem lies, - for working with shapes on large scales you should definetly choose a different CRS/Projection, from my point of view EPSG 4326 should be a good candidate.
Maybe that can be the next step: try to generate your shape in 4326 and measure it on a 4326 Map,
and see, if all is correct there.

cheers
Stefan

7

Hi Stefan

Thanks for your thoughts on this. I think they are starting to confirm my suspicions.

The issue manifests itself when, for example, I get our software application to programmatically create a 1 hectare square polygon (100m x 100m). At the equator a 100m x 100m square is created, but this decreases northwards, for example in the south of the UK a 62m x 62m polygon is created on the map. The measuring tool seems to work fine on whatever polygon is created (i.e. even though I tell the map to create 100 x 100m square and it creates a 62 x 62m square, the measuring tool says it is 62m and not 100m; likewise when I go to edit the polygon, 62m is shown).

The problem is I need to display satellite imagery (using Bing maps) in my application (so that users can see landscape features). I have tried changing the projection so that the polygons are created to the correct size, but then the Bing tiles fail to display correctly (they do not line up). So I think this aligns with the description of the problem you have outlined. As mentioned previously, currently my only solution so far has been to develop a ‘fudge factor’ based on latitude, which although not ideal, seems be ok over small distances, which is generally what the application is designed to do.

One thing that remains confusing is clearly the map is able to measure the distance ok using the measuring tool (i.e. using my example above), so I am bit baffled as to why when I programmatically tell it to create a 100m square it creates a 62m square. My only guess there are two different routines at work here, one of which is coping with the projection and which is not.

I welcome any further thoughts or suggestions on this.

Thanks

John

8

Hi John,

a from my point of view promising try would be, to create your shapefile in EPSG 4326 or a different local carthesian projection, and let your map stay in 3857. When you display your shapefile, you need to reproject it, - temporary to a different file, or in memory, and when you need coordinates from your map to change objects in your shapefile, you need to reproject them from EPSG 3857 to your shapefile crs before using them.
I work with this setup in my application, displaying osm tiles in the background, and it looks good so far.

cheers
Stefan

9

Hi Stefan

Thanks for this suggestion. I have tried something similar before, but without success. However, it’s possible that the order in which I did the process was perhaps wrong, so I may try to develop the routine from scratch again. The order I was doing things was:

  1. Open the application, load the background map (Bing Satellite) using the EPSG 3857 (Google Mercator) projection.
  2. Create a new shapefile with the same projection.
  3. Start the process of adding a shape (polygon)
  4. Save the new shapefile to a temporary file, and reproject it to EPSG 4326.
  5. Add the polygon
  6. Save the shape and shapefile.
  7. Save the temporary shapefile back to the file created in step 2 (overwrite it) and reproject it to EPSG 3857.

That resulted in the same output as before, i.e. 100m ending up as 62m.

If anything looks wrong with that process, please let me know.

Thanks

John

10

Hello Gentlemen.

I know I’m late to this party, but if it’s alright, I’ll give a little input.

As Stefan has suggested, it is not uncommon to set the map to 3857 to support the tiles, and maintain the layers in 4326 for data integrity. And set your GlobalSettings to AllowProjectionMismatch, and ReprojectLayersOnAdding.

So, regarding your steps, John, I would change step 2 to Create the new shapefile in 4326, replace steps 3, 4, and 5 with a single step, Add the polygon using DMS coordinates, and change step 7 to Add the new layer to map, letting the map reproject it to 3857, and see if that gives a better result.

Respectfully,
Jerry.

11

Hi Jerry and Stefan

Thanks for your suggestions on this issue last week, much appreciated. Since then I have been trying to implement a solution in my code based on your suggestions, and I think I have made some progress. I have discovered a few things in the process, so I thought I would share these with you.

  • Firstly, I have found that ESPG 4326 does not seem to cope that well with programmatically calculating distances, up to 10% out (over a 100m distance) in some instances. I have found that ESPG 3035 works better (at least for Europe, which is the context of my application).
  • Secondly, when reprojecting to and from 3857 and 3035 (and 4326) (using the utils.ReprojectShapefile function), I am finding the text in the PRJ file is not quite as it should be, thus was causing the projection to fail. I have easily fixed this by rewriting the PRJ file.
  • Thirdly, I can project a 3035 shapefile ok while having Bing maps displayed (using ESPG 3857), but if I try to programmatically calculate distances the measurements generated are derived from the ESPG 3857 projection. I have to clear the background map and projection settings, reload the 3035 shapefile, and then do any distance calculations. If there is workaround for this, do please let me know.
  • Finally, I did also find that trying to add polygons to a 3035 shapefile while projected on a 3857 map, resulted in the polygons not being saved – but this is not a major problem as I do not intend to use the 3035 shapefile as the interface for building the map.

I think that’s everything I’ve discovered so far. My next task will be to implement this in my main application (at the moment I’ve only played with the approach in a test app). If I encounter any further issues or have any additional questions, I’ll let you know.

Kind regards

John

12

Did you see this example already: MapWinGIS: CalculateArea.cs

When you reproject on the fly, the layer contains an in-memory version of the projected shapefile.
If you want to alter it you need to get it from the layer: sf = axMap1.get_Shapefile(layerHandle);
Now you can save it, if you like. But remember that this shapefile is projected to the map projection so don’t overwrite your original shapefile.