DotSpatial Projection with GridShift is Extremely Slow for NAD27

[cwgrc2]

Bug?

Not sure if this is a bug or expected behavior, but using a Projected.StatePlaneNad1927.NAD1927StatePlaneWyomingEastCentralFIPS4902
yields correct results with unacceptable times using a GridShift.

Steps to reproduce

1. For source , create: ProjectionInfo src = ProjectionInfo.FromEpsgCode(4326); // WGS84
2. For destination1, create: KnownCoordinateSystems.Projected.StatePlaneNad1927.NAD1927StatePlaneWyomingEastCentralFIPS4902
3. For destination2, create:Projected.StatePlaneNad1983.NAD1983StatePlaneWyomingEastCentralFIPS4902
   4: create 1500 lat long data points, in the range of Casper Wyoming (see code below)
   5: Reproject these points to dest1 and dest2. (see code below)
   6: The GridShift has to be turned on and the one I am using is: ntv1_can.dat (see code below)

DotSpatial version: 2.0.0-rc1

Expected behaviour

When I run my C# main with a NAD83 projection vs a NAD27(requiring a grid shift), NAD27 takes about 15 seconds to Reproject.ReprojectPoints for 1500 points. NAD83 is instantaneous. The NAD27 does give the right answer, but the speed is unacceptable.

Expected Behaviour: Instantaneous like the NAD83

Actual behaviour

Time to convert 1500 locations: 14.797 secs. Using NAD27 WITH Gridshifts
Time to convert 1500 locations: 0 secs. Using NAD83

Code in C# if Applicable:

using DotSpatial.Projections;
using System;
using System.Collections.Generic;
using System.Diagnostics;
using Utilities.Coordinates;

namespace CRSProjections
{
    class LatLongTest
    {

        static void Main(string[] args)
        {

            // Directly access the projected known systems for wyoming east cerntal, NAD27 & NAD83
            GridShift.InitializeExternalGrids("GridShifts", true);

            ProjectionInfo wyomingNad27 = KnownCoordinateSystems.Projected.StatePlaneNad1927.NAD1927StatePlaneWyomingEastCentralFIPS4902;
            ProjectionInfo wyomingNad83 = KnownCoordinateSystems.Projected.StatePlaneNad1983.NAD1983StatePlaneWyomingEastCentralFIPS4902;


            int srcEPSG = 4326; // WGS84 Geographic
            ProjectionInfo src = ProjectionInfo.FromEpsgCode(srcEPSG);

            int n = 1500;
            double latMin = 43.26808129519191;
            double latMax = 43.32542743693507;

            double longMin = -106.20071411132812;
            double longMax = -106.21667861938475;

            double latDiff = latMax - latMin;
            double longDiff = longMax - longMin;

            double latInc = latDiff / (double)n;
            double longInc = longDiff / (double)n;

            double[] xy = new double[n*2];
            double[] z = new double[n];  // Defaults to 'n' zeros
            int c = 0;
            for (int i = 0; i < n; i++)
            {
                double newLat = latMin + i * latInc;
                double newLong = longMin + i * longInc;
                xy[c++] = newLong;
                xy[c++] = newLat;
            }
            Stopwatch sw1 = Stopwatch.StartNew();
            Reproject.ReprojectPoints(xy, z, src, wyomingNad27, 0, z.Length);
            Console.WriteLine("Time to convert " + n + " locations: " + sw1.ElapsedMilliseconds / 1000.0 + " secs. Using NAD27 WITH Gridshifts");

            // reset points
            c = 0;
            for (int i = 0; i < n; i++)
            {
                double newLat = latMin + i * latInc;
                double newLong = longMin + i * longInc;
                xy[c++] = newLong;
                xy[c++] = newLat;
            }

            Stopwatch sw2 = Stopwatch.StartNew();
            Reproject.ReprojectPoints(xy, z, src, wyomingNad83, 0, z.Length);
            Console.WriteLine("Time to convert " + n + " locations: " + sw2.ElapsedMilliseconds / 1000.0 + " secs. Using NAD83 WITH Gridshifts");


            Console.ReadLine();

        }
    }
}

[VectorZita]

Not entirely sure about what is going on, but try to set:

GridShifts.ThrowGridShiftMissingExceptions = false

and run the code again. If you happen to see a faster conversion, then it might have something to do with an exception being thrown for each and every point, because the required grid is missing. As far as I have seen, the table you want is the conus grid-shift table.

Find the ...Nad\GridShifts folder that is being used in your application and unzip all files within the folder. Last time I checked (which was today), the GridShift.InitializeExternalGrids method does not search within zip files, only immediately available files (checked for presence by their file extensions). The ...Nad\GridShifts folder is located within the original DotSpatial.Projections folder of the main repository, but you may have relocated it to go with your application, so make sure to unzip the files in the folder you are actually using.

[cwgrc2]

Thank you for responding Vector:

I spent some time this morning testing out various scenarios around your recommendations.

Scenario 1: GridShifts Folder is Empty:

	GridShift.InitializeExternalGrids("GridShifts", true);
	GridShift.ThrowGridShiftMissingExceptions = true;
	
	Result: Threw exception: Grid shift is missing. For details go to https://github.com/DotSpatial/DotSpatial/wiki/Grid-Shifts

Scenario 2: GridShifts Folder is Still Empty:

	GridShift.InitializeExternalGrids("GridShifts", true);
	GridShift.ThrowGridShiftMissingExceptions = false;
	
	Result: Reports tons of spew to the debug window.  This is reported below as "Debug Report 1"

Scenario 3: GridShifts Folder has 4 Grid shift files in it: "alaska.lla, conus.lla, ntv1_can.dat, ntv2_0.GSB"

 I use these four because in the debug report from step 2, it warned that it tried to find/search for those 4

    GridShift.InitializeExternalGrids("GridShifts", true);
    GridShift.ThrowGridShiftMissingExceptions = true;
	
	Result: Good results but unacceptably slow.  

Scenario 4: Use same GridShift files as Step #3 above.

    GridShift.InitializeExternalGrids("GridShifts", true);
    GridShift.ThrowGridShiftMissingExceptions = false;
	
	Result: Exactly same as 3, Good results but unacceptably slow.  No difference by turning off the exceptions. 

Note that both steps 3 & 4 both reported the debug window as shown in 'Debug Report 2' shown below.

Final Note:

I tried removing the different files one-by-one in the GridShifts folder with the exceptions flag set to true. The ONLY file it did not throw an exception for was "ntv1_can.dat". So your recommendation for the "conus.lla" unfortunately did not work.

Any help you and your team can provide is greatly appreciated. But for now, I cannot use GridShifts for NAD27 in Wyoming because something is wrong with performance. Right answers, but too slow.

Thanks for your help!

Chris
scubagrant60@gmail.com
cell: 303-725-8323

Debug Report 1:

...
Unrecognized parameter skipped axis.
Unrecognized parameter skipped gamma.
ProjectionInfo.ParseProj4String: command 'no_uoff' not understood in '+proj=omerc +lat_0=45.30916666666666 +lonc=-86 +alpha=337.25556 +k=0.9996 +x_0=2546731.496 +y_0=-4354009.816 +no_uoff +gamma=337.25556 +datum=NAD83 +units=m +no_defs'
Unrecognized parameter skipped gamma.
ProjectionInfo.ParseProj4String: command 'no_uoff' not understood in '+proj=omerc +lat_0=45.30916666666666 +lonc=-86 +alpha=337.25556 +k=0.9996 +x_0=2546731.496 +y_0=-4354009.816 +no_uoff +gamma=337.25556 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs'
Unrecognized parameter skipped gamma.
ProjectionInfo.ParseProj4String: command 'no_uoff' not understood in '+proj=omerc +lat_0=4 +lonc=102.25 +alpha=323.0257905 +k=0.99984 +x_0=40000 +y_0=0 +no_uoff +gamma=323.1301023611111 +a=6377295.664 +b=6356094.667915204 +to_meter=20.116756 +no_defs'
Unrecognized parameter skipped gamma.
ProjectionInfo.ParseProj4String: command 'no_uoff' not understood in '+proj=omerc +lat_0=4 +lonc=102.25 +alpha=323.0257905 +k=0.99984 +x_0=804670.24 +y_0=0 +no_uoff +gamma=323.1301023611111 +a=6377295.664 +b=6356094.667915204 +units=m +no_defs'
Unrecognized parameter skipped gamma.
ProjectionInfo.ParseProj4String: command 'no_uoff' not understood in '+proj=omerc +lat_0=4 +lonc=102.25 +alpha=323.0257964666666 +k=0.99984 +x_0=804671 +y_0=0 +no_uoff +gamma=323.1301023611111 +ellps=GRS80 +units=m +no_defs'
Unrecognized parameter skipped gamma.
ProjectionInfo.ParseProj4String: command 'no_uoff' not understood in '+proj=omerc +lat_0=4 +lonc=115 +alpha=53.31580995 +k=0.99984 +x_0=0 +y_0=0 +no_uoff +gamma=53.13010236111111 +ellps=GRS80 +units=m +no_defs'
Unrecognized parameter skipped gamma.
ProjectionInfo.ParseProj4String: command 'no_uoff' not understood in '+proj=omerc +lat_0=57 +lonc=-133.6666666666667 +alpha=323.1301023611111 +k=0.9999 +x_0=5000000 +y_0=-5000000 +no_uoff +gamma=323.1301023611111 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs'
Unrecognized parameter skipped gamma.
ProjectionInfo.ParseProj4String: command 'no_uoff' not understood in '+proj=omerc +lat_0=45.30916666666666 +lonc=-86 +alpha=337.25556 +k=0.9996 +x_0=2546731.496 +y_0=-4354009.816 +no_uoff +gamma=337.25556 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs'
Unrecognized parameter skipped gamma.
ProjectionInfo.ParseProj4String: command 'wktext' not understood in '+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs'
ProjectionInfo.ParseProj4String: command 'wktext' not understood in '+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs'
...
...

Unrecognized parameter skipped vunits.
Unrecognized parameter skipped geoidgrids.
Unrecognized parameter skipped vunits.
Unrecognized parameter skipped geoidgrids.
Unrecognized parameter skipped vunits.
Unrecognized parameter skipped vunits.
Unrecognized parameter skipped vunits.
Unrecognized parameter skipped vunits.
...
...

Actually one for each points, so 1500 times...
...
pj_apply_gridshift(): failed to find a grid shift Table for location: (-106.21661476135253, 43.3251980523681). Names: @conus,@ntv1_can,@ntv2_0,@alaska
pj_apply_gridshift(): failed to find a grid shift Table for location: (-106.21662540435788, 43.32523628312926). Names: @conus,@ntv1_can,@ntv2_0,@alaska
pj_apply_gridshift(): failed to find a grid shift Table for location: (-106.21663604736325, 43.32527451389042). Names: @conus,@ntv1_can,@ntv2_0,@alaska
pj_apply_gridshift(): failed to find a grid shift Table for location: (-106.21664669036862, 43.32531274465158). Names: @conus,@ntv1_can,@ntv2_0,@alaska
pj_apply_gridshift(): failed to find a grid shift Table for location: (-106.21665733337402, 43.32535097541274). Names: @conus,@ntv1_can,@ntv2_0,@alaska
pj_apply_gridshift(): failed to find a grid shift Table for location: (-106.21666797637938, 43.32538920617391). Names: @conus,@ntv1_can,@ntv2_0,@alaska
...
...

Debug Report 2

nrecognized parameter skipped axis.
Unrecognized parameter skipped axis.
Unrecognized parameter skipped axis.
Unrecognized parameter skipped axis.
Unrecognized parameter skipped axis.
Unrecognized parameter skipped axis.
Unrecognized parameter skipped axis.
Unrecognized parameter skipped axis.
Unrecognized parameter skipped axis.
Unrecognized parameter skipped gamma.
ProjectionInfo.ParseProj4String: command 'no_uoff' not understood in '+proj=omerc +lat_0=45.30916666666666 +lonc=-86 +alpha=337.25556 +k=0.9996 +x_0=2546731.496 +y_0=-4354009.816 +no_uoff +gamma=337.25556 +datum=NAD83 +units=m +no_defs'
Unrecognized parameter skipped gamma.
ProjectionInfo.ParseProj4String: command 'no_uoff' not understood in '+proj=omerc +lat_0=45.30916666666666 +lonc=-86 +alpha=337.25556 +k=0.9996 +x_0=2546731.496 +y_0=-4354009.816 +no_uoff +gamma=337.25556 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs'
Unrecognized parameter skipped gamma.
ProjectionInfo.ParseProj4String: command 'no_uoff' not understood in '+proj=omerc +lat_0=4 +lonc=102.25 +alpha=323.0257905 +k=0.99984 +x_0=40000 +y_0=0 +no_uoff +gamma=323.1301023611111 +a=6377295.664 +b=6356094.667915204 +to_meter=20.116756 +no_defs'
Unrecognized parameter skipped gamma.
ProjectionInfo.ParseProj4String: command 'no_uoff' not understood in '+proj=omerc +lat_0=4 +lonc=102.25 +alpha=323.0257905 +k=0.99984 +x_0=804670.24 +y_0=0 +no_uoff +gamma=323.1301023611111 +a=6377295.664 +b=6356094.667915204 +units=m +no_defs'
Unrecognized parameter skipped gamma.
ProjectionInfo.ParseProj4String: command 'no_uoff' not understood in '+proj=omerc +lat_0=4 +lonc=102.25 +alpha=323.0257964666666 +k=0.99984 +x_0=804671 +y_0=0 +no_uoff +gamma=323.1301023611111 +ellps=GRS80 +units=m +no_defs'
Unrecognized parameter skipped gamma.
ProjectionInfo.ParseProj4String: command 'no_uoff' not understood in '+proj=omerc +lat_0=4 +lonc=115 +alpha=53.31580995 +k=0.99984 +x_0=0 +y_0=0 +no_uoff +gamma=53.13010236111111 +ellps=GRS80 +units=m +no_defs'
Unrecognized parameter skipped gamma.

Unrecognized parameter skipped axis.
Unrecognized parameter skipped axis.
Unrecognized parameter skipped axis.
ProjectionInfo.ParseProj4String: command 'no_uoff' not understood in '+proj=omerc +lat_0=4 +lonc=102.25 +alpha=323.0257905 +k=0.99984 +x_0=804671.2997750348 +y_0=0 +no_uoff +gamma=323.1301023611111 +a=6377304.063 +b=6356103.038993155 +towgs84=-11,851,5,0,0,0,0 +to_meter=20.11678249437587 +no_defs'
Unrecognized parameter skipped gamma.
ProjectionInfo.ParseProj4String: command 'no_uoff' not understood in '+proj=omerc +lat_0=57 +lonc=-133.6666666666667 +alpha=323.1301023611111 +k=0.9999 +x_0=5000000.001016002 +y_0=-5000000.001016002 +no_uoff +gamma=323.1301023611111 +datum=NAD27 +units=us-ft +no_defs'
Unrecognized parameter skipped gamma.
ProjectionInfo.ParseProj4String: command 'no_uoff' not understood in '+proj=omerc +lat_0=57 +lonc=-133.6666666666667 +alpha=323.1301023611111 +k=0.9999 +x_0=5000000 +y_0=-5000000 +no_uoff +gamma=323.1301023611111 +datum=NAD83 +units=m +no_defs'
Unrecognized parameter skipped gamma.
Unrecognized parameter skipped axis.
Unrecognized parameter skipped axis.

[VectorZita]

I will try to think a bit more to see if I can come up with something. The thing is, I have tried your code and, given that I have unzipped the files in the GridShifts folder and everything is properly set up, the conversion is almost instantaneous for me in both cases, with or without GridShifts.

I will keep looking into this.

[VectorZita]

It seems that if you only use conus.lla (remove all other files from the GridShifts folder), it will run almost immediately. However, I was able to reproduce the behavior you describe with ntv1_can.dat! It does take a lot longer, a bit over ~15 sec for a test attempt.

One thing that doesn't "exactly" make sense is that the code seems to prioritize grid shifts based on what is found first. I will also try to ascertain that this is the case and I will take a look to try to find the cause for the absurd amount of time.

[VectorZita]

It appears that the NadGrids for this Datum are searched in alphabetical order and, for the NAD23 they are conus, ntv1_can, ntv2_0 and alaska. Thus, if conus is found, it is used instead because it is "lexicographically" before all other.

The good thing is that I have found the problem. It turns out that this was what I was "afraid" it would be. The grid shifts file is being constantly re-read for data retrieval, while this should properly not needed after the first read. The bad thing is that, if you would like to fix the code, you will have to either edit and recompile the source code, or somehow "hack" your way around it (only if absolutely necessary).

I will post the exact source of the problem in a little while.

[jany-tenaj]

@VectorZita Can you create a pull request with a fix?

[VectorZita]

The DotSpatial.Projections.Nad.NadTable class contains a virtual void FillData() method, which is empty and needs to be overridden in derived classes. Most derived classes do this properly, but two don't (almost):

Class GbsNadTable overrides the method but forgets to set Filled = true; after assigning the Cvs variable. In my version, the assignment is in line 156. Immediately after that, Filled = true; is needed.

Class DatNadTable has the same problem. Immediately after line 105, Filled = true; is needed as well.

Because Filled = true; is missing, the method DotSpatial.Projections.GridShift.Convert(...) (in line 131 in the latest version) checks if (!table.Filled) table.FillData(); and, because table.Filled never becomes true, it is reloaded from the file continuously, when DatNadTable or GbsNadTable classes are used (i.e. when *.gbs or *.dat grid-shift files are employed).

@jany-tenaj Unfortunately, I have no experience with making pull requests as yet. Would this be possible for someone else to do? If not, I can give it a shot.

[jany-tenaj]

@cwgrc2 @VectorZita Ok check out the change.

[VectorZita]

@jany-tenaj it seems to be fine. I have tested the code outlined above by @cwgrc2 after fixing in a local repository on my machine and it no longer takes 15+ seconds to run, it is executed almost instantly (< ~0.05 secs) for the NAD27 conversion. Good job!

@cwgrc2 has to test and check if the results are OK, but they should be, since the only problem was the constant re-reading of the file, the calculations remain otherwise unaffected.

[cwgrc2]

Thanks all!!

Vector I could reproduce your step:

It seems that if you only use conus.lla (remove all other files from the GridShifts folder),
it will run almost immediately.

But... I had to set:

        GridShift.ThrowGridShiftMissingExceptions = false;

else it threw an exception. I tested those results and they look good and it is fast.

To get the fix, I have a noob question. I cloned the entire DotSDpatial-master (per @jany-tenaj - #1333) and I have loaded as a seperate Solution in Visual Studio 2019. How do I get that imported into my project as a Nuget package? I show at the bottom (see picture) how I did it originally, but I still see 2.0.0.-rc1 and I presume this is old and does not incorporate the new changes?

Thanks!!

Chris

[VectorZita]

As far as I understand, I think that a Nuget package will have to be re-created first, I don't know how this works for DotSpatial, or who can do it, specifically.

However, this shouldn't be necessary for you. You only need to build the project and use the compiled dll files instead, as references in your project. If you open the entire solution, try to build the DotSpatial.Projections project (you probably would want to set it to release mode, depending on your requirements). After successfully building it, you can get the dll file from the \bin\Release folder of the solution (unless you change the output path, of course) and add it as a reference to your own project (maybe a new Application?)

[VectorZita]

else it threw an exception. I tested those results and they look good and it is fast.

Did you test also using only the ntv1_can.dat? Because this is the specific case that was triggering the problem at all. The conus.lla file never posed a problem, only *.dat and *.gbs grid-shift files did.

[jany-tenaj]

@cwgrc2 Check out this wiki entry.

We build a new nuget package with every comit. But this doesn't get uploaded to nuget.

[cwgrc2]

Got the new release dll included. Works great! Thank you so much for all of your help. Good job.

Chris

[VectorZita]

No problem!

There is another bug I have detected in the DotSpatial.Projections library, which is rather subtle but has bitten me at some point in the past. If I remember correctly, it had to do with the Azimuthal Equidistant projection. Now that the Projections library has come up, I recalled it and I think it may be important to share it.

I will post a new issue relatively soon.

[jany-tenaj]

@VectorZita I'm guessing there are quite a lot of bugs in DotSpatial.Projections based on the fact that there were a lot of issues raised and it is an old port of proj.4.

[VectorZita]

@jany-tenaj Yes, I have traced the error I am referring to in the porting process, back to the original proj4 c++ source code (which is correct). I will try to document all this in a new issue as soon as possible.

[cwgrc2]

@jany-tenaj

I want to follow your steps in the wiki you pointed me to. I think I follow everything you do but am getting errors. I am trying to point to the continuous build nuget in your wiki. See picture below.

I get the error (see below) when I try and install this Dotspatial continuous feed nuget. Any ideas? Our build is failing on Azure DevOps pipelines.

Thanks

Chris

[dotspatial] The V2 feed at 'https://github.com/DotSpatial/DotSpatial/Search()?$filter=IsLatestVersion&searchTerm='Dotspatial'&targetFramework=''&includePrerelease=false&$skip=0&$top=26&semVerLevel=2.0.0' returned an unexpected status code '404 Not Found'.
[dotspatial] The V2 feed at 'https://github.com/DotSpatial/DotSpatial/Search()?$filter=IsAbsoluteLatestVersion&searchTerm='Dotspatial'&targetFramework=''&includePrerelease=true&$skip=0&$top=26&semVerLevel=2.0.0' returned an unexpected status code '404 Not Found'.
[dotspatial] The V2 feed at 'https://github.com/DotSpatial/DotSpatial/tree/master/Source/.nuget/Search()?$filter=IsLatestVersion&searchTerm='Dotspatial'&targetFramework=''&includePrerelease=false&$skip=0&$top=26&semVerLevel=2.0.0' returned an unexpected status code '404 Not Found'.
[dotspatial] The V2 feed at 'https://github.com/DotSpatial/DotSpatial/tree/master/Source/.nuget/Search()?$filter=IsAbsoluteLatestVersion&searchTerm='Dotspatial'&targetFramework=''&includePrerelease=true&$skip=0&$top=26&semVerLevel=2.0.0' returned an unexpected status code '404 Not Found'.
[Dotspatial.Projections] The V2 feed at 'https://github.com/DotSpatial/DotSpatial/Search()?$filter=IsAbsoluteLatestVersion&searchTerm='Dotspatial'&targetFramework=''&includePrerelease=true&$skip=0&$top=26&semVerLevel=2.0.0' returned an unexpected status code '404 Not Found'.