Projected coordinates do not match the shapefile coordinates

[Nitzan Peri]

Hi,

I am using a shapefile for local authority districts in the UK and trying to add data points from another source. I have the geographic longitude and latitude for each data point and I used Geo2xy to obtain the projected coordinates. However, the projected coordinates do not match the shapefile and the points appear far away from the map itself.

These are the details of the shapefile coordinates:

PROJCS["OSGB_1936_British_National_Grid",GEOGCS["GCS_OSGB 1936",DATUM["D_OSGB_1936",SPHEROID["Airy_1830",6377563.396,299.3249646]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",49],PARAMETER["central_meridian",-2],PARAMETER["scale_factor",0.9996012717],PARAMETER["false_easting",400000],PARAMETER["false_northing",-100000],UNIT["Meter",1]]


This is the code I used:

use "location_coords.dta", clear

geo2xy Latitude Longitude, generate(y_lat x_lon) projection(mercator, 6377563.396 299.3249646 0)

save "New_coords.dta", replace

use "LAD.dta", clear

spmap using "LAD_shp.dta", ///
id(_ID) ocolor(gs10 ..) osize(0.04 ..) ///
point(data("New_coords.dta") x(x_lon) y(y_lat) fcolor(red))


Thanks!
Nitzan

[Jared Greathouse]

Where is the shapefile located? What's the url?

[Nitzan Peri]

This is where the shapefile is from: https://geoportal.statistics.gov.uk/....960432%2C8.00

[Daniel Schaefer]

Hi Nitzan,

My guess is that you have two slightly different projections algorithms: one for the point data and one for the shape file. I don't know exactly what's going wrong here, but you might find Robert Picard's excellent discussion of geo2xy projections found in this thread useful.

[Nitzan Peri]

Thanks, I saw it, but the issue is that the shape file I'm using already has the projected coordinates, so I cannot recreate the projection. I'm uploading a screen shot of what I got.

[Daniel Schaefer]

Well, no. You can, in principle, change the projection of a shape file. In practice I feel this is relatively easy in QGIS and ArcGIS, but perhaps not so easy in Stata. I don't love that I'm becoming "that guy", but as convenient and well designed as these functions are, I feel this would be much more straightforward in ArcGIS, and easy enough in the open source alternative QGIS.

I suppose if you can't change the projection of the shapefile, you want to know how to get your point projection that you generate with geo2xy to match the shapefile projection? I'm not sure (maybe Jared knows), but I would start by pointing out the most obvious difference to me: mercator and transverse mercator are not the same thing.

[Scott Merryman]

The transverse Mercator is not the same as the Mercator projection. I am not sure what the LAT and LONG are referring to in the dbf file, but using those as an example you can use python within Stata to transform them to the same projection as the shapefile:

clear*
spshape2dta Local_Authority_Districts__April_2019__UK_BFE, saving(map) replace
use map

python:
from sfi import Data
x = Data.get('LONG')
y = Data.get('LAT')

import pyproj
#OSGB 1936 / British National Grid / Transverse_Mercator
transformer = pyproj.Proj("+proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy +datum=OSGB36 +units=m +no_defs ")
x1, y1 = transformer(x, y)
Data.addVarDouble("x1")
Data.addVarDouble("y1")
Data.store("x1", None, x1)
Data.store("y1", None, y1)
end
list in 1

. list in 1

     +-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+
     | _ID         _CX         _CY   OBJECTID     LAD19CD      LAD19NM   LAD19NMW    BNG_E    BNG_N                 LONG                  LAT                 Shape__Are              Shape__Len          x1          y1 |
     |-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------|
  1. |   1   448107.72   530696.61          1   E06000001   Hartlepool              447157   531476   -1.270230000000000   54.676200000000001   98441691.620788604021072   65270.260128191199328   447052.53   531494.13 |
     +-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+

[Jared Greathouse]

I'm also "that guy" now, but I would do (and have done) this in Python.

python: 
from osgeo import ogr from osgeo import osr  
   
   
  pointX = -11705274.6374 
  pointY = 4826473.6922
   
  # Spatial Reference System
  inputEPSG = 3857
  outputEPSG = 4326
   
  # create a geometry from coordinates
  point = ogr.Geometry(ogr.wkbPoint)
  point.AddPoint(pointX, pointY)
   
  # create coordinate transformation
  inSpatialRef = osr.SpatialReference()
  inSpatialRef.ImportFromEPSG(inputEPSG)
   
  outSpatialRef = osr.SpatialReference()
  outSpatialRef.ImportFromEPSG(outputEPSG)
   
  coordTransform = osr.CoordinateTransformation(inSpatialRef, outSpatialRef)
   
  # transform point
  point.Transform(coordTransform)
end

I think we'll do what you want, but note that you WILL have to change some of this code because the projection here isn't what yours is (I don't think).

[Nitzan Peri]

Thank you everyone!
I asked my colleague to run the projection code in Python and send me the coordinates. Now it works!
So glad I joined this forum