Android Question Spatialite

[Terradrones]

Hi All

I have read through Warwound's article, but I have no flippen idea how to use Spatialite to convert Lats\Longs to different Coordinate systems. Must be oldage or stupid...maybe both!

Any help or push in the right direction please?

 

[drgottjr]

i don't thing spatialite is used for the purpose you describe,
although i am happy to sit with you in the car reserved for the
old and stupid on the train and look out the window.

you may have an easier time searching for a specific conversion formula,
assuming you're not looking to convert the many thousands of
coordinate systems available. (apparently there are at least 6000!)

if you want to be able to convert from system to system at will,
you're going to have a hard time. especially for android.
there are some libraries out there, but links to them and/or to their
dependencies were disappointing. plus, you'd have to find someone to
wrap them, unless you're up to the task. out of curiosity i followed
the trail of a few of these libraries. i'm as interested in lat/lon as the
next guy, but only to a point. you might need to expand on exactly what
you're looking to do.

 

[Brian Dean]

If your objective is to convert Lat/Lon to a local Cartesian coordinate system then you are looking for the Transverse Mercator Projection. This is pretty heavy trigonometry, but I have a routine here from an old project of mine. You need to know the parameters of the Cartesian system you are targeting. There are hundreds of different systems because the world is not a true sphere and is anyway covered in lumps. Each Cartesian system is optimised fot the particular lump in its locality. Alternatively you can use UTM (Univeral Transverse Mercator) which does treat the Earth as smooth, but is not a widely used system.

Anyway, here is the code ...

Private Sub TRFtoGrid(Degrees As LatLon, a As Double, b As Double, N0 As Int, E0 As Int, F0 As Double, phi0 As Double, lambda0 As Double) As GridRef
'    TRF Lat/Lon To Grid East/North by Transverse Mercator projection.

'    <Phi> is latitude - negative if south of equator (all angles in radians)
'    <Lambda> is longitude - negative if west of prime meridian
'    <a> And <b> are major And minor axes of ellipsoid
'    <N0> And <E0> are easting And northing of grid origin
'    <F0> is scale factor on central meridian of grid
'    <phi0> And <lambda0> are latitude And longitude of grid origin
'    <East> And <North> are easting And northing of converted grid reference

'    Tested And confirmed 11th Oct 2018

    Dim Phi, Lambda As Double
    Dim n, n2, n3, e2, nu, rho, eta2 As Double                    ' Geometric terms
    Dim M, I, II, III, IIIa, IV, V, VI As Double          ' Intermediate terms
    Dim sin1, sin2, cos1, cos2, cos3, cos5 As Double
    Dim tan1, tan2, tan4 As Double                                            ' Trig's of <phi> and their powers
    Dim phidiff, phi2diff, phi3diff As Double                        ' (Phi - phi0), (Phi + phi0) and
    Dim phisum, phi2sum, phi3sum As Double                            '    their multiples
    Dim L, L2, L3, L4, L5, L6 As Double                                    ' (Lambda - lambda0) And its powers
    Dim x As Double                                                                            ' Temporary intermediate value term
    
    Dim result As GridRef

    Phi = Degrees.Latitude * RADIAN                                            ' Convert to radian measure
    Lambda = Degrees.Longitude * RADIAN

  n = (a - b) / (a + b)
  n2 = n * n
  n3 = n2 * n
  e2 = (a * a - b * b) / (a * a)                                            ' Eccentricity squared value

  ' Calculate trigonometric factors
    sin1 = Sin(Phi)    : cos1 = Cos(Phi)
  sin2 = sin1 * sin1
  cos2 = cos1 * cos1    : cos3 = cos2 * cos1    : cos5 = cos3 * cos2
  tan1 = sin1 / cos1    : tan2 = tan1 * tan1    : tan4 = tan2 * tan2

  x = (1 - e2 * sin2)
  nu = a * F0 / Sqrt(x)
  rho = a * F0 * ( 1 - e2) / Sqrt(x * x * x)

  eta2 = nu / rho - 1.0

  phidiff = (Phi - phi0)    : phi2diff = phidiff + phidiff
  phi3diff = phi2diff + phidiff
  phisum = (Phi + phi0)    : phi2sum = phisum + phisum
  phi3sum = phi2sum + phisum

  M = b*F0 * ( (1 + n + 1.25 * (n2 + n3)) * phidiff _
                - (3 * (n + n2) + 2.625 * n3) * Sin(phidiff) * Cos (phisum) _
                + 1.875 * (n2 + n3) * Sin(phi2diff) * Cos(phi2sum) _
                - 35/24 * n3 * Sin(phi3diff) * Cos(phi3sum))
  I = M + N0
  x = nu * sin1
  II = x * cos1 * 0.5
  III = x * cos3 * (5 - tan2 + 9*eta2) / 24
  IIIa = x * cos5 * (61 - 58*tan2 + tan4) / 720
  IV = nu * cos1
  V = nu * cos3 * (nu / rho - tan2) / 6
  VI = nu * cos5 * (5 - 18*tan2 + tan4 + 14*eta2 - 58*tan2*eta2) / 120

  L = (Lambda - lambda0)
  L2 = L * L    : L3 = L2 * L    : L4 = L3 * L    : L5 = L4 * L    : L6 = L5 * L

    result.Northing = I + II * L2 + III * L4 + IIIa * L6
    result.Easting = E0 + IV * L + V * L3 + VI * L5
    Return result   
End Sub

If you think this is what you are looking for then let me know, and tell me what part of the world you are living in so that I might be able to give you more help.

 

[Terradrones]

Thank you so much Brain. This is what I am looking for.

I stay in Paarl, about 50km north of Cape Town in South Africa.

 

[Terradrones]

Do you perhaps have a routine from Grid to Lat\Long?

 

[Terradrones]

Brian, is there a library missing? Degrees as LatLon?

Sorry, I spelt your name wrong!

 

[Brian Dean]

Brian, is there a library missing? Degrees as LatLon? Sorry, I spelt your name wrong!

Don't worry about getting my name wrong; it happens quite a lot and nobody would object to being called "Brain".

My earlier post was just checking that I understood what you wanted before I spent too much time on it. The code was taken from a class. The class includes a routine for Reverse Transverse Mercator. It also defines two Types to hold Lat/Lon values and Grid Coordinate pairs.

I am attaching a modified class module. The original was designed for the National Grid used here in the UK, and I have left in some sample entry points that. These should should show you what is required. If you are using WGS84 then the TRF factors (a and b) will be the same, but you will have to supply your local grid parameters.

Let me know if you need anything more.

 

[Terradrones]

Hi Brian, thank you very much for that. I will need more help and am willing to pay for it via Paypal.

 

[Brian Dean]

Sorry @Terradrones - I missed your post. Just catching up now.

You don't say what sort of help you want, but I would guess that it is about grid systems rather than about B4X. The advent of GPS has disrupted traditional mapping systems, which are generally based on older geodetic datums, all over the world. I have looked briefly at the South Africa situation and that seems typical.

Give me some idea of what help you think you need. If it is about mapping then this forum is not the place to discuss it. Maybe we can work something out. Anyway, don't worry about Paypal; anything I can tell you you would be able to find out for yourself if you knew where to look.

 

[Terradrones]

Hi Brian, do you perhaps have a database (SQLite) that contains all the constants for the different coordinate systems?

 

[Brian Dean]

No, I do not know of any such database. There are many hundreds, maybe thousands, of Cartesian coordinate systems, and projections other than Transverse Mercator. I don't think it would be logical to try to create a single database - one database per country or region might be hard enough.

Assuming that you are mainly interested in South Africa the document here is worth a read. It explains the current reference system using Gauss Conform coordinates which is a modification of Transverse Mercator and produces southings and westings rather than eastings and northings. The document provides a full set of conversion formulae.

In the introduction to the document is this sentence - "Numerous map projections and coordinate systems are used in South Africa, especially for mapping purposes." I guess that there are many maps around that do not use the current standard.

Here is another reasonably readable document that covers much the same ground. Again, it provides full formulae to convert between Gauss Conform coordinates and Latitude/Longitude. I hope this is helpful.

 

[Terradrones]

Thanks Brian. I have a SQLite Database that covers most of the countries that I used in VB.net compact, but I cannot remember what I did there. Old age!!

I tried to attach the file here, but it would not allow me. Maybe I can email it to you?

 

[Brian Dean]

It seems to me that you probably know more about this topic than I do. The work I have done has been mainly in the British Isles with a bit in Australia where UTM got me through.

Thanks for the offer of the database, but I would not have a use for it. If you want to thank me just "Like" some of my posts - that is reward enough for me.

 

[Green1]

Spatialite can transform geometries, given input and output SRIDs. For a lat/lon point you can do something like

SELECT X(utmPt) AS X, Y(utmPt) AS Y
FROM (SELECT Transform(MakePoint([LON], [LAT], [SRID]), [outputSRID]) AS utmPt) AS Pt

Note order: LON, LAT : Easting, Northing (X,Y)

If you're only projecting to/from UTM and no other spatial functions, then Erel's UTM class might be easier to incorporate and considerably lighter and simpler

https://www.b4x.com/android/forum/t...onvert-lat-lon-utm-coordinates.30702/#content

 

[Terradrones]

Hi Green1, could you elaborate a bit more please. What is "SRID" and "outputSRID"?

 

[Green1]

You didn't provide any detail, so it was generic

If you are transforming between spatial reference (coordinate) systems, you need to know what systems they are. SRID is the Spatial Reference ID, sometimes known as EPSG and you can find these at epsg.io
Spatialite transforms geometries so it needs to create a geometry with a known spatial reference MakePoint(x, y, sr). The Transform function takes the resultant geometry as it's first argument and the output spatial reference as the second.
The X, Y query returns the coordinates from the resultant transformed geometry of the sub query

Only you know which systems you're using.
As mentioned, if you're only projecting WGS84 (GPS) to UTM, Erel's solution and example does it all for you

 

[gvoulg]

Continuing Green1's answer , this ia how I use Spatialite to make the transformations from lat,lon to x,y and back.

'''transform lat lon to local
Sub get_point_to_Local( lat As Double , lon As Double) As String
Dim sqStm As Spatialite_TableResult
Dim result As String
Dim sql1 As String
sql1 = "SELECT ST_AsText(ST_Transform(ST_GeomFromText('POINT("
sql1 = sql1 & lon & " " & lat & ")',4326),"
sql1 = sql1 & Starter.local_srid & ")) As trans_geom;"
sqStm = SpatialiteDatabase.GetTable(sql1)
result = sqStm.Rows(0,0)
Return result
End Sub

'''transform χ y local to wgs84
Sub get_point_from_local( x As Double , y As Double) As String
Dim sqStm As Spatialite_TableResult
Dim result As String
Dim sql1 As String
sql1 = "SELECT ST_AsText(ST_Transform(ST_GeomFromText('POINT("
sql1 = sql1 & x & " " & y & ")'," & Starter.local_srid & "),4326)) As trans_geom; "
Try
sqStm = SpatialiteDatabase.GetTable(sql1)
result = sqStm.Rows(0,0)
Catch
Log(LastException.Message)
MsgboxAsync(LastException.Message & CRLF & "wrong EPSG code" ,"ATTENTION")
Starter.local_srid="2100"
Return "0000000000000000000"
End Try
If result=Null Or result="" Then
Return "0000000000000000000"
End If
Return result
End Sub

you need a spatialite database in your project to make the calculations and dim it like

Private SpatialiteDatabase As Spatialite_Database

to use it you have to initialize it like this

Initialize(data_dir ,data_file,"READONLY")

and the init function is

Public Sub Initialize (DPATH As String ,DBASE As String , OPEN_STYLE As String)
If Starter.debugmode Then Log ("accessing " & DPATH &"/" & DBASE & " as " & OPEN_STYLE)
SpatialiteDatabase.Initialize
Select OPEN_STYLE
Case "READWRITE"
SpatialiteDatabase.Open(DPATH, DBASE, SpatialiteConstants.SQLITE_OPEN_READWRITE)
Case "READONLY"
SpatialiteDatabase.Open(DPATH, DBASE, SpatialiteConstants.SQLITE_OPEN_READONLY)
End Select
End Sub

local_srid is the epsg code of yor system.
2100 is the system in Greece.
Hope that helps
George

 

[Terradrones]

Continuing Green1's answer , this ia how I use Spatialite to make the transformations from lat,lon to x,y and back.



you need a spatialite database in your project to make the calculations and dim it like


to use it you have to initialize it like this

and the init function is

local_srid is the epsg code of yor system.
2100 is the system in Greece.
Hope that helps
George

Hi George

Thank you for the help.

Just one more question: where do I get the SpatialiteDatabase?

 

[gvoulg]

here is a lot of info

The Gaia-SINS federated project home-page

GAIA - Spatial Is Not Special - libraries and tool

www.gaia-gis.it

you need the spatialite-gui.exe to handle this type of sqlite databases (spatial) (too big to attach)
George

 

[Terradrones]

Hi All

I need a push again please.

Gvoulg has helped me a lot, but I am doing something wrong.

In the "Addcoords" Activity I have a routine to plot the coordinates on Google Map:

[/
Sub PlotPoints(A As Int)
    Dim AveY, AveX As Double
    Dim cp As CameraPosition
    Dim Co As CircleOptions
  
    File.MakeDir(File.DirRootExternal & "/CEASER/","Data/Projections")
  
    Co.Initialize
    'Dim BitMap1 As Bitmap  = LoadBitmapResize(File.DirRootExternal & "/CEASER/DATA/", "Point.png",8dip,8dip,True)
    Try
        gmap.Clear
  
        If CGlobals.CoordCode=1 Then
            'Site Coords
            Dim rs As ResultSet = CGlobals.sql1.ExecQuery("SELECT Name, East, North, Elevation, Description FROM SCoords")
        else if CGlobals.CoordCode=2 Then
            'Global Coords
            Dim rs As ResultSet = CGlobals.sql1.ExecQuery("SELECT Name, East, North, Elevation, Description FROM GCoords")
        Else
            'Job Coords
            Dim rs As ResultSet = CGlobals.sql1.ExecQuery("SELECT Name, East, North, Elevation, Description FROM Coords")
        End If
        i=0: AveY=0: AveX=0
      
        Do While rs.NextRow
            Dim row(5) As Object
            row(0) = rs.GetString("Name")
            row(1) = NumberFormat2(rs.GetDouble("East"),1,3,3,False)
            row(2) = NumberFormat2(rs.GetDouble("North"),1,3,3,False)
            'row(3) = NumberFormat2(rs.GetDouble("Elevation"),1,3,3,False)
            Geo.get_point_from_local(row(1),row(2))
            Co.StrokeWidth(2)
'            co.Center2(Engine.Lat,Engine.Lon).Radius(Radius.SelectedItem).FillColor(Colors.White).StrokeColor(Colors.DarkGray)
'            gmapextra.AddCircle(gmap,co)
            'gmap.AddMarker(Engine.Lat,Engine.Lon, "New Marker")
            'gmap.AddMarker3(Engine.lat,Engine.Lon,row(0),Bitmap1)
            i=i+1
            AveY=AveY+row(1)
            AveX=AveX+row(2)
        Loop
        rs.Close
        If (AveY<>0 Or AveX<>0) And A=0 Then
            'Zoom in to center of points
            AveY=AveY/i
            AveX=AveX/i
            Convert(row(1),row(2),34)
'            Engine.Hart94_YX(AveY,AveX,19)
'            cp.Initialize(Engine.Lat,Engine.Lon, 10)
'            gmap.MoveCamera(cp)
            'ConvertLatLong
        End If
    Catch
        Log(LastException)
    End Try
  
End Sub


]

In my Geodesy Class (which I initialized as "Geo"), I have the following code:

[/

Sub Class_Globals
    Private Spa As Spatialite_Database

    Dim PI As Double = 3.1415926535897932385
    Dim RADIAN As Double = 2.0 * PI / 360    
    
    Dim Local_srid As String
    'Dim Data_File As String
    Dim spconstants As Spatialite_Constants

'    Initialize(File.DirRootExternal & "/CEASER/Data/" ,Data_File,"READONLY")
'    Spa.Open(File.DirInternal, "mm.sqlite", spconstants.SQLITE_OPEN_READONLY)

End Sub

Public Sub Initialize '(DPATH As String ,DBASE As String , OPEN_STYLE As String)
    Dim DPath, DBase, Open_Style As String

    Spa.Initialize
    
    DPath=File.DirRootExternal & "/CEASER/Data"
    DBase="Projections.sl3"
    Open_Style="ReadOnly"
    
    Spa.Open(DPath, DBase, spconstants.SQLITE_OPEN_READONLY)

    Select Open_Style
        Case "READWRITE"
            Spa.Open(DPath, DBase, spconstants.SQLITE_OPEN_READWRITE)
        Case "READONLY"
            Spa.Open(DPath, DBase, spconstants.SQLITE_OPEN_READONLY)
    End Select
End Sub

Sub get_point_from_local( x As Double , y As Double) As String
    Dim sqStm As Spatialite_TableResult
    Dim result As String
    Dim sql1 As String
    
    sql1 = "SELECT ST_AsText(ST_Transform(ST_GeomFromText('POINT("
    sql1 = sql1 & x & " " & y & ")'," & CGlobals.EPSG & "),4326)) As trans_geom; "
    Try
        sqStm = Spa.GetTable(sql1)
        result = sqStm.Rows(0,0)
    Catch
        Log(LastException.Message)
        MsgboxAsync(LastException.Message & CRLF & "wrong EPSG code" ,"ATTENTION")
        Local_srid="2048"
        Return "0000000000000000000"
    End Try
    If result=Null Or result="" Then
        Return "0000000000000000000"
    End If
    Return result
End Sub

'''transform lat lon to local
Sub get_point_to_Local( lat As Double , lon As Double) As String
    Dim sqStm As Spatialite_TableResult
    Dim result As String
    Dim sql1 As String
    sql1 = "SELECT ST_AsText(ST_Transform(ST_GeomFromText('POINT("
    sql1 = sql1 & lon & " " & lat & ")',4326),"
    sql1 = sql1 & CGlobals.EPSG & ")) As trans_geom;"
    sqStm = Spa.GetTable(sql1)
    result = sqStm.Rows(0,0)
    Return result
End Sub
]

The EPSG code is correct, as I load it from my Database.

It keeps returning 0, when I want to convert from a Local System to Lats\Longs.