Is That Georeferencing Really Wrong?

TL;DR: What looks like a georeferencing error in a file may really be a software bug. Check independently.

I’ve gotten a few questions lately claiming that some of our geospatial data in the wrong place by several miles. I’ve also seen similar questions on a forum related to another agency’s data. While it is certainly possible that we’ve made an error somewhere, it’s always a good idea to check more than one software before deciding the georeferencing is wrong. In those recent cases, it’s a software bug and the data is fine. In this post, I’m going to run through a couple examples and illustrate how you might check with free tools. I do want to emphasis that while my examples are from one bug in one software, it’s not an attack on that software. All software has bugs.

Example 1: Raster

Our first case is about the digital elevation models in Tiff format for Berkeley County, SC. Someone wrote the following kind note to alert us of a problem.

“I wanted to make you aware the TIFs currently available for download are projected incorrectly. The data comes in several miles south of it’s actual location. I will be happy to send a screen capture if that will be beneficial. Thank you!”

Naturally, putting out bad data is not what we want to do, so we quickly took a look at our files. They are all in NAD83(2011) South Carolina State Plane. Running the gdalsrsinfo command from the open source GDAL package provides the following information:

PROJCS["NAD83(2011) / South Carolina (ft)",
    GEOGCS["NAD83(2011)",
        DATUM["NAD83_National_Spatial_Reference_System_2011",
            SPHEROID["GRS 1980",6378137,298.257222101,
                AUTHORITY["EPSG","7019"]],
            AUTHORITY["EPSG","1116"]],
        PRIMEM["Greenwich",0,
            AUTHORITY["EPSG","8901"]],
        UNIT["degree",0.0174532925199433,
            AUTHORITY["EPSG","9122"]],
        AUTHORITY["EPSG","6318"]],
    PROJECTION["Lambert_Conformal_Conic_2SP"],
    PARAMETER["standard_parallel_1",34.83333333333334],
    PARAMETER["standard_parallel_2",32.5],
    PARAMETER["latitude_of_origin",31.83333333333333],
    PARAMETER["central_meridian",-81],
    PARAMETER["false_easting",2000000],
    PARAMETER["false_northing",0],
    UNIT["foot",0.3048,
        AUTHORITY["EPSG","9002"]],
    AXIS["X",EAST],
    AXIS["Y",NORTH],
    AUTHORITY["EPSG","6570"]]

I’ve put a couple items I’ll refer back to later in bold. The main things we can see here is that GDAL thinks this is NAD83(2011) South Carolina data in feet, and uses the EPSG authority code of 6570. While you can’t tell from the gdalsrsinfo output, all that info about the parameters isn’t in the file at all. The only thing there about the projection is EPSG code 6570. You can see that with the listgeo program that comes with the GeoTiff libraries. So far, everything looks good. Those parameters match what they should be.

From prior experience I suspected that ArcGIS of some flavor might have been in use, which turned out to be true. I took a look at what ArcGIS Pro (version 2.4.0) says about the data’s projection.

Spatial Reference information as seen in ArcGIS Pro

That looks pretty good, except for the latitude of origin. Everything else agrees. In the report from gdalsrsinfo and in the information from EPSG, the latitude of origin should be 31.8333333 (or 31 degrees 50 minutes) for South Carolina State Plane. ArcGIS Pro is reporting 31.5 degrees. We would expect that difference in projection definition to throw everything off by around 20 miles in the north south direction.

What’s going on here? Why are we seeing two different latitude of origin when that’s not a parameter that is embedded in the file at all? I can only guess at the answer, but I think it’s something like a table for the state plane NAD83(2011) parameters for EPSG codes being entered in degrees minutes seconds and interpreted as decimal degrees. Thus, 31 degrees 50 minutes is being interpreted as 31.5 degrees.

Example 2: Lidar

For a second example, let’s take that case on the LAStools forum that I referred to at the top. There we have an LAZ file in NAD83(2011) CONUS Albers (EPSG 6350). This time we’ll use lasinfo from LAStools to look at the file, but we could also use the info command from PDAL. Here’s the full set of keys in the georeferencing header:

variable length header record 1 of 2:
reserved 0
user ID 'LASF_Projection'
record ID 34735
length after header 96
description 'GeoTiff GeoKeyDirectoryTag'
GeoKeyDirectoryTag version 1.1.0 number of keys 11
key 1024 tiff_tag_location 0 count 1 value_offset 1 - GTModelTypeGeoKey: ModelTypeProjected
key 1025 tiff_tag_location 0 count 1 value_offset 1 - GTRasterTypeGeoKey: RasterPixelIsArea
key 1026 tiff_tag_location 34737 count 43 value_offset 0 - GTCitationGeoKey: NAD83(2011) / Conus Albers + NAVD88 height
key 2049 tiff_tag_location 34737 count 12 value_offset 43 - GeogCitationGeoKey: NAD83(2011)
key 2054 tiff_tag_location 0 count 1 value_offset 9102 - GeogAngularUnitsGeoKey: Angular_Degree
key 3072 tiff_tag_location 0 count 1 value_offset 6350 - ProjectedCSTypeGeoKey: NAD83(2011) / Conus Albers
key 3076 tiff_tag_location 0 count 1 value_offset 9001 - ProjLinearUnitsGeoKey: Linear_Meter
key 4096 tiff_tag_location 0 count 1 value_offset 5703 - VerticalCSTypeGeoKey: NAVD88 height (Reserved EPSG)
key 4097 tiff_tag_location 34737 count 14 value_offset 55 - VerticalCitationGeoKey: NAVD88 height
key 4098 tiff_tag_location 0 count 1 value_offset 5103 - VerticalDatumGeoKey: Vertical Datum Codes 5103
key 4099 tiff_tag_location 0 count 1 value_offset 9001 - VerticalUnitsGeoKey: Linear_Meter

Yes, that’s kind of long. Things to note are 1) key 3072 (in bold) is the projection with EPSG code 6350 and 2) the standard parallels not there. If we look up 6350 on the EPSG registry, we’ll see that the standard parallels should be 29.5 and 45.5. Let’s look at what ArcGIS shows when we bring it in.

Projection information for LAS file in ArcGIS Pro

Here we see a few of interesting things. First, standard parallel 1 is wrong at 29.3 degrees. This matches the idea of 29.5 being entered as degrees minutes seconds (which is 29 degrees 30 minutes) and then interpreted as decimal degrees (29.3). Second, the second standard parallel is just fine. Third, even though an EPSG authority was in the file, it isn’t showing here. But, it’s picking up so much other information from that 6350 code that it can’t be that it didn’t understand it. As you’d guess, having a projection parameter wrong will put all the data in the wrong place.

Conclusion

My point here is not really about ArcGIS, which is a fine software package. I’m only using it as the example because so many people use it. I’ve seen other packages that have errors on projections too. I remember one had two of the Florida State Plane zones swapped. Everybody gets bugs in the mix and they can be subtle. For this particular bug, it appears that only the NAD83(2011) entries have a problem. The NAD83(NSRS2007) entries all seem to fine. My point is that you might want to double check what your software is telling you with something else.

Finally, how do you deal with this particular ArcGIS bug? If you run the project tool to define the projection as what you know it should be, it will set it. It’s not reprojecting, it’s just storing the information in a way that it understands correctly. In case you’re wondering, this bug was reported to ESRI and I found it was a bug they already had listed.