GEOREF Support

Question

GEOREF Support

matt-stomper opened this issue 2 years ago · 11 comments

Is your feature request related to a problem? Please describe.
This request is not part of a problem with CoordinateSharp, it is simply to add support for GEOREF https://en.wikipedia.org/wiki/World_Geographic_Reference_System

Describe the solution you'd like

Adaption of https://github.com/NationalSecurityAgency/qgis-latlontools-plugin/blob/master/georef.py
GeoRef.ToGeoRef(Coordinate Coordinate)
GeoRef.ToString(int precision) - outputs 2-12 characters
Coordinate.Parse(string georef) - returns the center of the GeoRef pending precision of the GeoRef string
GeoRef.ToGeoFence(string georef)
GeoFence.Get_List_GeoRef(int precision) - returns a List within a GeoFence area

The ToGeoFence() function aims to GeoFence a GEOREF based on its precision. A 4 character GEOREF will equal a 111km x 111km GeoFence area. This allows for the use of IsPointInPolygon().

Describe alternatives you've considered
I have not considered any alternative. CoordinateSharp is a fantastic library. I would love to see this library continue to grow and at best, replace the Microsoft Location Class.

Additional context
The below code is a quick 5 minute adaption of qgis-latlontools-plugin. I have not researched the licensing if adapting this code to C# is allowed. I will continue researching this.

More research needs to go into how to store the additional precision provided by a decimal place 1 minute value.

namespace CoordinateSharp.GeoRef
{
    public partial class GeoRef 
    {
        public string DegreeQuad15Latitude { get; set; }

        public string DegreeQuad15Longitude { get; set; }

        public string DegreeQuad1Latitude { get; set; }

        public string DegreeQuad1Longitude { get; set; }

        public string Minute1Latitude { get; set; }

        public string Minute1longitude { get; set; }
    }
}

namespace CoordinateSharp.GeoRef
{
    public partial class GeoRef
    {
        private const string _digits = "0123456789";
        private const string _lontile = "ABCDEFGHJKLMNPQRSTUVWXYZ";
        private const string _lattile = "ABCDEFGHJKLMM"; // Repeat the last M for 90 degrees which rounds up - Prevents extra checks in the code
        private const string _degrees = "ABCDEFGHJKLMNPQ";
        private const int _tile = 15;
        private const int _lonorig = -180;
        private const int _latorig = -90;
        private const int _base = 10;
        private const int _baselen = 4;
        private const int _maxprec = 11;
        private const int _maxlen = _baselen + 2 * _maxprec;

        /// <summary>
        /// Converts a GeoRef to a Geofence
        /// </summary>
        public static void ToGeoFence()
        {

        }

        public static bool ToGeoRef(Coordinate coordinate, int precision, out GeoRef georef)
        {
            georef = new GeoRef();

            var latitude = coordinate.Latitude.DecimalDegree;
            var longitude = coordinate.Longitude.DecimalDegree;

            if (latitude == 90d)
            {
                latitude -= float.Epsilon;
            }

            precision = Math.Abs(precision);
            if (precision == 1)
            {
                //Disallow prec = 1
                precision += 1;
            }

            var m = 60000000000;
            var x = (int)Math.Floor(longitude * m) - _lonorig * m;
            var y = (int)Math.Floor(latitude * m) - _latorig * m;
            var ilon = (int)x / m;
            var ilat = (int)y / m;
            georef.DegreeQuad15Longitude = _lontile[(int)(ilon / _tile)].ToString();
            georef.DegreeQuad15Latitude = _lattile[(int)(ilat / _tile)].ToString();
            if (precision >= 0)
            {
                georef.DegreeQuad1Longitude = _degrees[(int)(ilon % _tile)].ToString();
                georef.DegreeQuad1Latitude = _degrees[(int)(ilat % _tile)].ToString();
                var c = 0;
                if (precision > 0)
                {
                    x = (int)(x - m * ilon);
                    y = (int)(y - m * ilat);
                    var d = Math.Pow(_base, _maxprec - precision);
                    x = (int)(x / d);
                    y = (int)(y / d);
                    c = precision;
                }

                //var georef1 = new char[c];
                //while (c > 0)
                //{
                //    georef.Minute1longitude +=  1[c] = _digits[(int)(x % _base)];
                //    x = (int)(x / _base);
                //    georef1[(int)(c + precision)] = _digits[(int)(y % _base)];
                //    y = (int)(y / _base);
                //    c = c - 1;
                //}
            }

            //return (''.join(georef1))
            return true;
        }

    }

}

I'm happy to collaborate further.

Answer 1 · 2022-09-13T20:19:01.000Z

I think this is a great idea, and I appreciate the provided code (and your kind words). The QGIS adaption carries a GPL license which would be compatible with our AGPL, but it is not compatible with our Commercial License waiver. Their code appears to be adapted from Charles Karney's work however, which carries a compatible MIT license. We would most likely adapt from his research https://geographiclib.sourceforge.io/doc/research.html#geodesics.

If it's good enough for the NSA it's good enough for us.

We will look into feasibility and performance impacts of this implementation.

Answer 2 · 2022-09-30T18:03:58.000Z

Sample algorithms created. Performance impact of conversions appear minimal. Will get to working adding to CoordinateSharp soon.

//Example class
class GEOREF
{
        public string Quad_15 { get; set; }
        public string Quad_1 { get; set; }
        public string Easting { get; set; } //Storing as string due to zero padding in system
        public string Northing { get; set; } //Storing as string due to zero padding in system
        public int Precision { get; set; } = 8;
        
        public override string ToString()
        {
            return $"{Quad_15}{Quad_1}{Easting}{Northing}";
        }
}

//Example Program
 static string digits = "0123456789";
 static string lngTile = "ABCDEFGHJKLMNPQRSTUVWXYZ";
 static string latTile = "ABCDEFGHJKLMM"; // Repeat the last M for 90 degree check efficiency
 static string degrees = "ABCDEFGHJKLMNPQ";
 static int tile = 15;
 static double lngorig = -180;
 static double latorig = -90;
 static double based = 10;
 static double maxprec = 11;

 static void Main(string[] args) {
   var georef = ToGEOREF(45.075869, 65.691258);
   ToLatLong(georef);
 }

 static GEOREF ToGEOREF(double lat, double lng) {
   GEOREF georef = new GEOREF();
   string easting = "";;
   string northing = "";

   if (lat == 90) {
     lat = lat - double.Epsilon;
   }
   if (lat == -90) {
     lat = lat + double.Epsilon;
   }

   long m = 60000000000;

   long x = (long)(Math.Floor(lng * m) - lngorig * m);
   long y = (long)(Math.Floor(lat * m) - latorig * m);
   long ilon = (long)(x / m);
   long ilat = (long)(y / m);

   georef.Quad_15 = lngTile[(int)(ilon / tile)].ToString() + latTile[(int)(ilat / tile)].ToString();

   georef.Quad_1 = degrees[(int)(ilon % tile)].ToString() + degrees[(int)(ilat % tile)].ToString();

   x = (long)(x - m * ilon);
   y = (long)(y - m * ilat);
   double d = Math.Pow(based, maxprec - georef.Precision);
   x = (long)(x / d);
   y = (long)(y / d);
   int c = georef.Precision;

   while (c > 0) {
     easting = easting.Insert(0, digits[Math.Abs((int)(x % based))].ToString());
     x = (long)(x / based);

     northing = northing.Insert(0, digits[Math.Abs((int)(y % based))].ToString());
     y = (long)(y / based);

     c -= 1;
   }

   georef.Easting = easting;
   georef.Northing = northing;

   Console.WriteLine(georef);

   return georef;
 }

 public static void ToLatLong(GEOREF georef) {
   double lon1 = lngTile.IndexOf(georef.Quad_15[0]) + lngorig / tile;
   double lat1 = latTile.IndexOf(georef.Quad_15[1]) + latorig / tile;

   double unit = 1;

   if (georef.Precision > 0) {
     unit = unit * tile;

     lon1 = lon1 * tile + degrees.IndexOf(georef.Quad_1[0]);
     lat1 = lat1 * tile + degrees.IndexOf(georef.Quad_1[1]);

     int i = 0;

     while (i < georef.Precision) {
       double m = based;
       if (i == 0) {
         m = 6;
       }

       unit = unit * m;
       int x = digits.IndexOf(georef.Easting[i].ToString());
       int y = digits.IndexOf(georef.Northing[i].ToString());

       lon1 = m * lon1 + x;
       lat1 = m * lat1 + y;

       i++;
     }

     double lat = (tile * lat1) / unit;
     double lon = (tile * lon1) / unit;

     Console.WriteLine(lat + " " + lon);
   }

 }

Answer 3 · 2022-11-09T21:22:22.000Z

Quick status update. This is still on our radar, but we've had to reprioritize our focus to handle some time sensitive projects. Might be a couple months before we can finish this one up. The conversions are built, but the georef to geofence still needs some work. We may do a split release and get the conversions out first before the geofence features.

Will keep you posted.

Answer 4 · 2022-11-14T00:29:33.000Z

Hi,

Thanks for the update. If the release was split, we could work on the GeoFence functionality for you. I am working on an add function GeoFence.Add(Coordinate coordinate) that takes a point and sorts the points from bottom left -> top left -> top right -> bottom right. The aim is to check if a point is in GeoFence, if not, add it and sort the new list.

Thanks,
Matt.

Answer 5 · 2022-11-19T00:29:00.000Z

@matt-stomper this sounds like a great plan.

I will work on documenting the GEOREF and try to get the initial part out ASAP (should have a beta out in the next week).

Answer 6 · 2022-11-21T01:45:27.000Z

@matt-stomper

BETA pushed to nuget if you want to start testing it out: https://www.nuget.org/packages/CoordinateSharp/2.15.1.1-beta

Includes GEOREF conversions, but not GeoFence features as previously discussed.

Examples

Coordinate c = new Coordinate(45, 64);
//Default output set at precision 6
Console.WriteLine(c.GEOREF.ToString()); //SKEA000000000000

//Override output precision via override
Console.WriteLine(c.GEOREF.ToString(2)); //SKEA0000

You may parse from a string or build a GEOREF using the four parts of the coordinate.

Coordinate c = Coordinate.Parse("SKEA425698152658"); //Recommend TryParse, but using parse for brevity

//Or

GEOREF geo = new GEOREF("SK","EA","425698","152658");
Coordinate c = GEOREF.ConvertGEOREFtoLatLong(geo); //N 45º 15' 15.948" E 64º 42' 34.188"

Will release a stable version after an initial testing period.

Answer 7 · 2022-12-13T01:56:09.000Z

GEOREF conversion have been included in v2.15.2.1 released on nuget.org.

This release completes the conversion part of this feature request. This issue will remain open until the geofence feature has been implemented or a determination is made to not implement.

As you stated, it sounds like you are working to the geofence logic potentially. Such a contribution (should you share the code) would warrant free commercial licensing per our licensing policy.

We may decide to work on it as well in the near future time permitting, so we will keep you posted on our efforts via this issue so as to not duplicate them.

Thank you for the feature suggestion!

Answer 8 · 2023-05-06T16:59:28.000Z

Status update:

This should be complete with the next release (2.18). Just need to implement testing and correct any issues that appear with drawn boxes crossing in to the previous precision level (maxing out seconds, minutes or quad 1 parameters).

Examples:

int precision = 6;
//Get GeoFence
georef.ToGeoFence(precision);

//Get Corners
georef.Get_BottomLeftCorner(precision);
georef.Get_TopRightCorner(precision);

Pole tips may not be reliable (which should be fine as long as it's documented).

We also discovered GEOREFs could overproject beyond 60 minutes during creation in the easting and northing levels. Some systems automatically adjust when this occurs, but we are going to restrict to 59.999999 minutes / seconds to coincide with the standards of the library.

Will close this issue once the update is published

Answer 9 · 2023-05-07T22:11:22.000Z

All features released with v2.18.1.1.

Let us know if anything behaves unexpectedly.

Answer 10 · 2023-05-07T23:58:32.000Z

Hi Team,

Sorry for being on radio silence and not responding with your updates. The work you have achieved is amazing. We have done some work on the GeoFence component but needs further testing. We will haven't a PR for that feature by the end of the month.

Thanks again team. Your efforts are very appreciated.

Answer 11 · 2023-05-08T01:20:04.000Z

@matt-stomper no worries. Just an FYI, we built out the GeoFence component if you want to test it out with 2.18.1.1. If you plan to build your own however, consider a PR still as you may have considered things we didn't.