Search code examples
sql-serverleafletspatialmap-projectionssqlgeography

Change projection in MSSQL for web mapping (Leaflet,Openlayer, OpenStreetMaps, GoogleAPI, ...) to WSG48 or any other format


I have some WKT/WKB data in the MSSQL server like this and would like to show them on the map with the help of leaflet, Openlayer, OpenStreetMaps, or GoogleAPI. My data look likes this:

POLYGON ((1736946.0983 5923253.9175, 1736895.6852 5923333.9451, 1736936.0082 5923356.6991, ......))

in this format

 EPSG:2193

as shown below

enter image description here and would like to convert them to

  WGS48 or EPSG:4326

I am trying to add them to leaflet, and seems that I need to convert this data to an appropriate format like this:

[[42.353770, -71.103606], [42.355447, -71.104475], [42.362681, -71.089830], [42.361829, -71.079230]]

Unfortuanlty, I am not sure how to do this conversion.

I have tries some approaches like these

[1.]. Error:

 A .NET Framework error occurred during execution of user-defined routine or aggregate "geography": 
 System.FormatException: 24201: Latitude values must be between -90 and 90 degrees.

[2.] Query:

select  GEOGRAPHY::STPolyFromText ([stastxt],4326) from mytable
Error:
 A .NET Framework error occurred during execution of user-defined routine or aggregate "geography": 
 System.FormatException: 24201: Latitude values must be between -90 and 90 degrees.

[3.] No success

[4.] I am able to see this data as a layer in GeoServer when I set the Declared SRS as EPSG:2193

After doing more study on this issue, the question is:

Can I do this transformation merely in MSSQL server or Leaflet or I need to use some other tools like Proj4net and replace '(' with '['?


Solution

  • Finally, I developed an interesting solution, I will explain my solution step by step as I am sure other faces similar problems.

    The first point is to know which projection the data is (src) and which projection you would like to have (dst). Usually, dst is EPSG:4326, EPSG:3857, or WGS48. For this solution, I needed to find the correct mathematics behind that, so I used this website https://mygeodata.cloud/cs2cs/ to find the correct format for src and dst (If you are familiar with R it has a function for this also called spTransform). Another reason is that I gonna use Proj4 component for this conversion. Another critical step is conversion to GeoJson as these web maps can read GeoJson files. I did not want to write my data to physical GeoJson files, so I did the conversion on demand in MSSQL (No need to write it somewhere as a GeoJson file).

    1. Make a library in C#, and get a dll.
    2. Import this to MSSQL
    3. Enjoy it

    The code for the first part: Install DotSpatial via Nugget

     using Microsoft.SqlServer.Server;
     using System;
     using System.Collections;
     using System.Collections.Generic;
     using System.Linq;
     using System.Text;
     using System.Threading.Tasks;
    
    
    public class classval
    {
        public string line;
        public string src;
        public string dst;
    
        public classval(string line, string src, string dst)
        {
            this.line = line;
            this.src = src;
            this.dst = dst;
        }
    
    
    
    }
    
    public class CLRProjection
    {
        private static IEnumerable<classval> ConvertedEnumerable(string line, string src, string dst)
        {
            return new List<classval> { new classval(line, src, dst) };
        }
    
        [SqlFunction(FillRowMethodName = "FillRow")]
        public static IEnumerable ToLatLong(string Geometry, string src, string dst)
        {
            return ConvertedEnumerable(Geometry, src, dst);
        }
    
        private static void FillRow(Object classvalobj, out string Geometry, out string srcprj, out string dstprj)
        {
        
           classval geomobj = (classval)classvalobj;
            string _geometry = geomobj.line; //"POLYGON ((1736946.0983 5923253.9175,....))";
        string proj4_src = geomobj.src; //"+proj=tmerc +lat_0=0 +lon_0=173 +k=0.9996 +x_0=1600000 +y_0=10000000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs ";
        string proj4_dst = geomobj.dst;//"+proj=longlat +datum=WGS84 +no_defs";
        _geometry = _geometry.Replace(",", " , ");
        _geometry = _geometry.Remove(0, _geometry.IndexOf('('));
        _geometry = _geometry.Replace("(", "[ ");
        _geometry = _geometry.Replace(")", " ]");
        string[] splitbycomma = _geometry.Split(',');
        
        foreach (var itembycomma in splitbycomma)
            {
    
                string tmpitem = itembycomma;
                tmpitem = tmpitem.Replace('[', ' ');
                tmpitem = tmpitem.Replace(']', ' ');
                tmpitem = tmpitem.Trim();
                string[] splitbyspace = tmpitem.Split(' ');
                for (int ibs = 0; ibs < splitbyspace.Length - 1; ibs++)
                {
    
                    double[] x = { double.Parse(splitbyspace[ibs]) };
                    double[] y = { double.Parse(splitbyspace[ibs + 1]) };
                    double[] z = new double[x.Length];
                    //rewrite xy array for input into Proj4
                    double[] xy = new double[2 * x.Length];
                    int ixy = 0;
                    for (int i = 0; i <= x.Length - 1; i++)
                    {
                        xy[ixy] = x[i];
                        xy[ixy + 1] = y[i];
                        z[i] = 0;
                        ixy += 2;
                    }
                    double[] xy_geometry = new double[xy.Length];
                    Array.Copy(xy, xy_geometry, xy.Length);
    
    
                    
                    DotSpatial.Projections.ProjectionInfo src =
                        DotSpatial.Projections.ProjectionInfo.FromProj4String(proj4_src);
                    DotSpatial.Projections.ProjectionInfo trg =
                        DotSpatial.Projections.ProjectionInfo.FromProj4String(proj4_dst);
    
                    DotSpatial.Projections.Reproject.ReprojectPoints(xy, z, src, trg, 0, x.Length);
    
    
                    ixy = 0;
                    for (int i = 0; i <= x.Length - 1; i++)
                    {
                    _geometry = _geometry.Replace(xy_geometry[ixy].ToString() + " ", "[" + xy[ixy + 1].ToString() + " , ");
                    _geometry = _geometry.Replace(xy_geometry[ixy + 1].ToString() + " ", xy[ixy].ToString() + " ] ");
                    _geometry = _geometry.Replace("- ", "-");
                        string tt = (i + 1 + " " + xy[ixy] + " " + xy[ixy + 1]);
    
                        ixy += 2;
                    }
    
                }
            }
        _geometry = _geometry.Replace("  ", " ");
        _geometry = _geometry.Replace(" [ ", "[");
        _geometry = _geometry.Replace(" ] ", "]");
        _geometry = _geometry.Replace(" , ", ",");
        srcprj = proj4_src;
        dstprj = proj4_dst;
        Geometry = _geometry;
        }
    
    }
    

    The code for the second part (Inside MSSQL)

      ALTER DATABASE test SET trustworthy ON
     CREATE ASSEMBLY CLRFunctionAssem
     FROM N'C:\Users\...\bin\Debug\Convertor_Projection.dll'
     WITH PERMISSION_SET = UNSAFE
     GO
    
    
     CREATE FUNCTION dbo.ToLatLong(@Geometry nvarchar(max), @src nvarchar(max),@dst nvarchar(max))
     RETURNS TABLE
     ( _geom  nvarchar(max) ,srcprj  nvarchar(max) ,dstprj  nvarchar(max) 
    ) with execute as caller
    AS
     EXTERNAL NAME CLRFunctionAssem.[CLRProjection].[ToLatLong]
    

    MSSQL code

     SELECT       
      [parcelid]
      ,[Geom1]
      ,[stastxt]
       ,conv._geom
       
     FROM [test].[dbo].[TEST_JSON] as a  CROSS APPLY dbo.ToLatLong (a.
    [stastxt],'+proj=tmerc +lat_0=0 +lon_0=173 +k=0.9996 +x_0=1600000 +y_0=10000000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs','+proj=longlat +datum=WGS84 +no_defs') as conv
     where [stastxt] is not null
    

    MSSQL2016 has JSON functiality, while older versions do not have this ability.

    Output enter image description here

    Resources which helped me are:

    1,2,3