Announcement

Collapse
No announcement yet.

lat long wgs84 to british national grid

Collapse
X
 
  • Filter
  • Time
  • Show
Clear All
new posts

  • #21
    Thanks Stuart, with that adjustment things run but the results are 0,0,0 on the export coordinate test. I have SourceSRID=4936 which should be UTM cartesian, TargetSRID=27700 which should be BNG, SourceRevision=2015 and TargetRevision=2015. Last detail?


    #COMPILE EXE
    #DIM ALL
    #INCLUDE "win32api.inc"

    TYPE Coordinates
    X AS DOUBLE
    Y AS DOUBLE
    Z AS DOUBLE
    END TYPE

    DECLARE FUNCTION ConvertCoordinates _
    LIB "GIQ.dll" ALIAS "ConvertCoordinates" _
    (BYVAL SourceSRID AS INTEGER,_
    BYVAL TargetSRID AS INTEGER,_
    BYVAL SourceRevision AS INTEGER,_
    BYVAL TargetRevision AS INTEGER,_
    BYREF InputCoordinates AS Coordinates,_
    BYREF OutputCoordinates AS Coordinates,_
    BYREF Datum AS INTEGER) AS LONG


    FUNCTION PBMAIN () AS LONG

    LOCAL SourceSRID AS INTEGER
    LOCAL TargetSRID AS INTEGER
    LOCAL SourceRevision AS INTEGER
    LOCAL TargetRevision AS INTEGER
    LOCAL Datum AS INTEGER

    DIM InputCoordinates AS coordinates
    DIM OutputCoordinates AS coordinates

    SourceSRID=4936
    TargetSRID=27700
    SourceRevision=2015
    TargetRevision=2015
    datum=0

    'in utm
    inputcoordinates.x=529405.163
    inputcoordinates.y=5700001.532
    inputcoordinates.z=0

    CALL convertcoordinates (SourceSRID,TargetSRID,SourceRevision,TargetRevision,InputCoordinates,OutputCoordinates,Datum)
    MSGBOX STR$(outputcoordinates.x)+STR$(outputcoordinates.y)+STR$(outputcoordinates.z)

    END FUNCTION

    Comment


    • #22
      Two things:
      1. Your Revisions were reversed according to the Access module code
      2.X and Y coordinates should be radians, not degrees:

      Output:
      ---------------------------
      PowerBASIC
      ---------------------------
      438481 567433 -49
      ---------------------------
      OK
      ---------------------------




      '
      Code:
      #COMPILE EXE
      #DIM ALL
      #INCLUDE "win32api.inc"
      
      TYPE Coordinates
          X AS DOUBLE
          Y AS DOUBLE
          Z AS DOUBLE
      END TYPE
      
      DECLARE FUNCTION ConvertCoordinates  _
          LIB "GIQ.dll" ALIAS "ConvertCoordinates" (BYVAL SourceSRID AS INTEGER,_
          BYVAL TargetSRID AS INTEGER,_
          BYVAL SourceRevision AS INTEGER,_
          BYVAL TargetRevision AS INTEGER,_
          BYREF InputCoordinates AS Coordinates,_
          BYREF OutputCoordinates AS Coordinates,_
          BYREF Datum AS INTEGER) AS LONG
      
      
      FUNCTION PBMAIN () AS LONG
      
          LOCAL sourcesrid AS INTEGER
          LOCAL targetsrid AS INTEGER
          LOCAL sourcerevision AS INTEGER
          LOCAL targetrevision AS INTEGER
          LOCAL datum AS INTEGER
          LOCAL lngResult AS LONG
          DIM inputcoordinates AS coordinates
          DIM outputcoordinates AS coordinates
      
          sourcesrid=4937
          targetsrid=27700
          sourcerevision=0
          targetrevision=2015
          datum=1
      
          inputcoordinates.x= DegToRad(-1.4)
          inputcoordinates.y= DegToRad(55)
          inputcoordinates.z= 0
      
          lngResult = convertcoordinates (sourcesrid,targetsrid,sourcerevision,targetrevision,inputcoordinates,outputcoordinates,datum)
          ? USING$ ("#  #  #",outputcoordinates.x,outputcoordinates.y,outputcoordinates.z)
      
      END FUNCTION
      
      FUNCTION DegToRad(d AS DOUBLE) AS DOUBLE
          FUNCTION = 3.1415926 * d / 180
      END FUNCTION
      
      FUNCTION RadToDeg(r AS DOUBLE) AS DOUBLE
          FUNCTION = r * 180 / 3.1415926
      END FUNCTION
      '
      Last edited by Stuart McLachlan; 12 Apr 2020, 09:56 PM. Reason: Edited: Output not Rad so applying RadToDeg gave incorrect values

      Comment


      • #23
        Thanks Stuart...Still examining as those numbers on the test are out of bounds...

        Comment


        • #24
          Originally posted by dean goodman View Post
          Thanks Stuart...Still examining as those numbers on the test are out of bounds...
          Numbers were out because I was not familiar with the coordinate systems and I incorrectly applied RadToDeg to the output coordinates

          Code above updated.

          Comment


          • #25
            Stuart, thats it!!! You are very talented and see all the scenarios.... Thank you so much for spending this extra time on this. The above code you modified others can now use (along with the zip library and other necessary files) to get things into British National Grid. This was very unusual to have a dynamic georeference system that requires periodic updates and tries to account for seafloor spreading and other phenomena. But the British do things their own ways! Thanks again...

            Comment


            • #26
              Just for something to do, I converted the various functions in the .mdb example to a PB include file.

              They compile, but I haven't actually tested them

              '
              Code:
              'GIQFunctions.inc
              $OSGridLetters  = "ABCDEFGHJKLMNOPQRSTUVWXYZ"
              %OSGridLeft = 0
              %OSGridBottom = 0
              %OSGridRight = 1000000
              %OSGridTop = 1500000
              
              FUNCTION CoordinatesWithinRect(X AS DOUBLE, Y AS DOUBLE, L AS DOUBLE, B AS DOUBLE, T AS DOUBLE, R AS DOUBLE) AS LONG
                IF Y < B THEN
                  CoordinatesWithinRect = %False
                ELSEIF Y > T THEN
                  CoordinatesWithinRect = %False
                ELSEIF X < L THEN
                  CoordinatesWithinRect = %False
                ELSEIF X > R THEN
                  CoordinatesWithinRect = %False
                ELSE
                  CoordinatesWithinRect = %True
                END IF
              END FUNCTION
              
              FUNCTION CoordinatesToOSGridLetters(X AS DOUBLE, Y AS DOUBLE) AS STRING
                LOCAL EastingIndex AS INTEGER
                LOCAL NorthingIndex AS INTEGER
                LOCAL GridLetters AS STRING
                IF CoordinatesWithinRect(X, Y, %OSGridLeft, %OSGridBottom, %OSGridRight, %OSGridTop) THEN
                  ON ERROR RESUME NEXT
                  ' Determine the 500Km letter.
                  EastingIndex = 2 + INT(X * 0.000002)
                  NorthingIndex = 3 - INT(Y * 0.000002)
                  GridLetters = MID$($OSGridLetters, 1 + EastingIndex + (NorthingIndex * 5), 1)
                  ' Add the 100Km letter.
                  EastingIndex = INT(X * 0.0001) MOD 5
                  NorthingIndex = 4 - (INT(Y * 0.0001) MOD 5)
                  GridLetters = GridLetters + MID$($OSGridLetters, 1 + EastingIndex + (NorthingIndex * 5), 1)
                  ON ERROR GOTO 0
                  IF LEN(GridLetters) = 2 THEN
                    FUNCTION = GridLetters
                  ELSE
                    FUNCTION = ""
                  END IF
                ELSE
                  FUNCTION = ""
                END IF
              END FUNCTION
              
              FUNCTION OSGridLettersToCoordinates(GridLetters AS STRING,Result AS Coordinates) AS LONG
                LOCAL TestOSGridLetters AS STRING
                LOCAL Index AS INTEGER
                ' Ensure upper case grid letters.
                TestOSGridLetters = UCASE$(TRIM$(GridLetters))
                ' Convert 500Km letter to coordinates.
                Index = ASC(LEFT$(TestOSGridLetters, 1)) - ASC("A")
                IF Index > 8 THEN
                  Index = Index - 1 ' Skip I in alphabetical sequence.
                END IF
                Result.X = ((Index MOD 5) - 2) * 500000
                Result.Y = (3 - (Index \ 5)) * 500000
                ' Add the 100Km letter to them.
                Index = ASC(MID$(TestOSGridLetters, 2, 1)) - ASC("A")
                IF Index > 8 THEN
                  Index = Index - 1 ' Skip I in alphabetical sequence.
                END IF
                Result.X = Result.X + ((Index MOD 5) * 100000)
                Result.Y = Result.Y + ((4 - (Index \ 5)) * 100000)
              END FUNCTION
              
              FUNCTION OSGridReferenceToCoordinates(GridReference AS STRING, Result AS Coordinates) AS LONG
                LOCAL TestOSGridReference AS STRING
                LOCAL CoordinateLength AS INTEGER
                LOCAL EastingPart AS STRING
                LOCAL NorthingPart AS STRING
                LOCAL lngResult AS LONG
                ' Ensure the removal of any white space.
                TestOSGridReference = REMOVE$(GridReference, " ")
                ' Convert the grid letters into the 100Km square coordinate.
                lngResult = OSGridLettersToCoordinates(LEFT$(TestOSGridReference, 2),Result)
                ' Determine the length of the given coordinates.
                CoordinateLength = (LEN(TestOSGridReference) - 2) \ 2
                ' Split the reference into the coordinate parts.
                EastingPart = MID$(TestOSGridReference, 3, CoordinateLength)
                NorthingPart = MID$(TestOSGridReference, CoordinateLength + 3, CoordinateLength)
                ' If the length is less than 5, pad the parts to 5 digits.
                IF CoordinateLength < 5 THEN
                  EastingPart = EastingPart & STRING$(5 - CoordinateLength, "0")
                  NorthingPart = NorthingPart & STRING$(5 - CoordinateLength, "0")
                END IF
                ' If the length is greater than 5, insert a decimal point after the first five.
                IF CoordinateLength > 5 THEN
                  EastingPart = MID$(EastingPart, 1, 5) & "." & MID$(EastingPart, 6, CoordinateLength - 5)
                  NorthingPart = MID$(NorthingPart, 1, 5) & "." & MID$(NorthingPart, 6, CoordinateLength - 5)
                END IF
                Result.X = Result.X + VAL(EastingPart)
                Result.Y = Result.Y + VAL(NorthingPart)
              END FUNCTION
              
              FUNCTION CoordinatesToOSGridReference(X AS DOUBLE, Y AS DOUBLE, OPTIONAL Figures AS INTEGER) AS STRING
                IF ISMISSING(Figures) THEN Figures = 6
                LOCAL CoordinateLength AS INTEGER
                LOCAL Multiplier AS DOUBLE
                LOCAL Result AS STRING
                Result = CoordinatesToOSGridLetters(X, Y)
                IF LEN(Result) = 2 THEN
                  CoordinateLength = Figures \ 2
                  Multiplier = (10 ^ (5 - CoordinateLength))
                  Result = Result & MID$(FORMAT$(Multiplier * ROUND(X / Multiplier,0), "00000000000"), 7, CoordinateLength) & _
                                    MID$(FORMAT$(Multiplier * ROUND(Y / Multiplier,0), "00000000000"), 7, CoordinateLength)
                ELSE
                  Result = ""
                END IF
                FUNCTION = Result
              END FUNCTION
              
              FUNCTION IsValidOSGridReference(GridReference AS STRING) AS LONG
                LOCAL TestOSGridReference AS STRING
                LOCAL TextLength AS INTEGER
                LOCAL Result AS LONG
                LOCAL TestValue AS DOUBLE
                TestOSGridReference = REMOVE$(GridReference, " ")
                TestOSGridReference = UCASE$(REMOVE$(TestOSGridReference, ","))
                TextLength = LEN(TestOSGridReference)
                Result = ((TextLength MOD 2) = 0) AND (TextLength > 0) AND (TextLength < 14)
                IF Result THEN
                  IF (INSTR("HJNOST", LEFT$(TestOSGridReference, 1)) = 0) OR (INSTR("ABCDEFGHJKLMNOPQRSTUVWXYZ", MID$(TestOSGridReference, 2, 1)) = 0) THEN
                    Result = %False
                  END IF
                  IF Result THEN
                    ON ERROR RESUME NEXT
                    TestValue = VAL(MID$(TestOSGridReference, 3, TextLength))
                    Result = (ERR = 0)
                    ON ERROR GOTO 0
                  END IF
                END IF
                FUNCTION = result
              END FUNCTION
              '

              Comment


              • #27
                And....should someone need the utmtolatlong operation here is a routine i use.....where longitude1,latitude1and zonenumber1 are global variables..

                FUNCTION utmtolatlong AS LONG

                LOCAL northing AS DOUBLE
                LOCAL easting AS DOUBLE
                LOCAL longitude AS DOUBLE
                LOCAL latitude AS DOUBLE
                LOCAL cm AS DOUBLE
                LOCAL ol AS DOUBLE
                LOCAL pi AS DOUBLE
                LOCAL a AS DOUBLE
                LOCAL b AS DOUBLE
                LOCAL e12 AS DOUBLE
                LOCAL e AS DOUBLE
                LOCAL f AS DOUBLE
                LOCAL e1 AS DOUBLE
                LOCAL ep2 AS DOUBLE
                LOCAL m AS DOUBLE
                LOCAL n AS DOUBLE
                LOCAL mu AS DOUBLE
                LOCAL mo AS DOUBLE
                LOCAL j1 AS DOUBLE
                LOCAL j2 AS DOUBLE
                LOCAL j3 AS DOUBLE
                LOCAL j4 AS DOUBLE
                LOCAL fp AS DOUBLE
                LOCAL c1 AS DOUBLE
                LOCAL t1 AS DOUBLE
                LOCAL r1 AS DOUBLE
                LOCAL n1 AS DOUBLE
                LOCAL d AS DOUBLE
                LOCAL q1 AS DOUBLE
                LOCAL q2 AS DOUBLE
                LOCAL q3 AS DOUBLE
                LOCAL q4 AS DOUBLE
                LOCAL q5 AS DOUBLE
                LOCAL q6 AS DOUBLE
                LOCAL q7 AS DOUBLE
                LOCAL k AS DOUBLE
                LOCAL zonenumber AS LONG

                zonenumber=ABS(zonenumber1)

                pi = 4 * ATN(1)
                cm = (6*zonenumber-183) / 180 * pi
                ol = 0 / 180 * pi

                'ellipsoid parameters k, a, f (currently wgs-84)
                k = 0.9996
                a = 6378137
                f = 1/298.257223563

                b = (1-f)*a
                e = SQR((2-f)* f)
                ep2 = e^2/(1-e^2)
                e1 = (1 - SQR(1 - e^2))/(1 + SQR(1 - e^2))

                mo = a * ( (1 - e^2/4 - 3*e^4/64 - 5*e^6/256) * ol _
                - (3*e^2/8 + 3*e^4/32 + 45*e^6/1024) * SIN(2 * ol) _
                + (15*e^4/256 + 45*e^6/1024) * SIN(4 * ol) _
                - (35*e^6/3072) * SIN(6 * ol))

                northing = utmnorthing
                easting = 500000 - utmeasting

                IF zonenumber1<0 THEN
                Northing =utmnorthing-10000000' //10000000 meter offset FOR southern hemisphere
                END IF


                m = mo + northing / k
                mu = m / (a * (1 - e^2/4 - 3*e^4/64 - 5*e^6/256))


                j1 = 3*e1/2 - 27*e1^3/32
                j2 = 21*e1^2/16 - 55*e1^4/32 '
                j3 = 151*e1^3/96 '
                j4 = 1097*e1^4/512 '

                fp = mu + j1 * SIN(2*mu) + j2 * SIN(4*mu) + j3 * SIN(6*mu) + j4 * SIN(8*mu) '


                c1 = ep2 * COS(fp) ^ 2 '
                t1 = TAN(fp) ^ 2 '
                n1 = a / (1 - e^2 * SIN(fp)^2) ^ 0.5 '
                r1 = a * (1 - e^2) / (1 - e^2 * SIN(fp)^2) ^ 1.5 '
                d = easting / (n1 * k) '

                q1 = (n1 * TAN(fp)) / r1 '
                q2 = (d^2 / 2) '
                q3 = (5 + 3 * t1 + 10 * c1 - 4 * c1^2 - 9 * ep2) * d ^ 4 / 24
                q4 = (61 + 90 * t1 + 298 * c1 + 45 * t1^2 - 3 * c1^2 - 252 * ep2) * d ^ 6 / 720

                q5 = d
                q6 = (1 + 2*t1 + c1) * d ^ 3 / 6
                q7 = (5 - 2*c1 + 28*t1 - 3*c1^2 + 8*ep2 + 24*t1^2) * d^5/120

                latitude = fp - q1 * (q2 - q3 + q4) '
                longitude = cm - (q5 - q6 + q7) / COS(fp) '

                longitude1 = longitude * 180 / pi
                latitude1 = latitude * 180 / pi

                FUNCTION = 1

                END FUNCTION

                Comment

                Working...
                X