10

Using ArcPy, how do you check if two feature classes have the same spatial reference?

Just checking if the two are equal doesn't work:

>>> import arcpy
>>> fc1 = r"C:\Users\e1b8\Desktop\E1B8\GIS_Stackexchange\data.gdb\test"
>>> sr1 = arcpy.Describe (fc1).spatialReference 
>>> sr2 = arcpy.Describe (fc1).spatialReference
>>> sr1 == sr2
False

factoryCode doesn't work, because custom projections don't have them.

>>> fc2 = r"C:\Users\e1b8\Desktop\E1B8\GIS_Stackexchange\data.gdb\customproj"
>>> sr2 = arcpy.Describe (fc2).spatialReference
>>> sr2.factoryCode
0

There's name, but names can be the same, but have different units:

>>> sr1 = arcpy.Describe (fc1).spatialReference
>>> sr2 = arcpy.Describe (fc2).spatialReference
>>> sr1.name
u'NAD_1983_UTM_Zone_10N'
>>> sr2.name
u'NAD_1983_UTM_Zone_10N'
>>> sr1.linearUnitCode
9003
>>> sr2.linearUnitCode
9001

So it gets a bit complicated. The best I've come up with is:

>>> def CompareSRs (inFc1, inFc2):
    sr1 = arcpy.Describe (inFc1).spatialReference
    sr2 = arcpy.Describe (inFc2).spatialReference
    if not sr1.name != sr2.name:
        return False
    srType = sr1.type
    if srType != sr2.type:
        return False
    if srType == "Geographic":
        return sr1.angularUnitCode == sr2.angularUnitCode
    return sr1.linearUnitCode == sr2.linearUnitCode

And I'm still not sure the above code is air tight.

Is there a better way?

PolyGeo
  • 65,136
  • 29
  • 109
  • 338
Emil Brundage
  • 13,859
  • 3
  • 26
  • 62

2 Answers2

14

Judging from comments, you might have it already :)

You could compare the Well-Known Text (WKT) descriptions of the spatial references.

sr1 = arcpy.Describe(dataset1).spatialReference
sr2 = arcpy.Describe(dataset2).spatialReference
sr1String = sr1.exportToString()
sr2String = sr2.exportToString()

matching = False

if sr1String == sr2String:
    # Exact string match
    matching = True
else:
    # difference
    pass
Paulo Raposo
  • 1,930
  • 14
  • 21
3

In case 2019, I had similar issues using ArcMap 10.3 and wanted to be as sure as possible as to whether projections matched. As in the above question/answers, you can get the spatial reference using arcpy.Describe(dataset).spatialReference. In a function library of mine, I then integrate this into a workflow, set-up to handle the comparison of 2 datasets.

The individual attributes of a geoprocessing spatial reference object are available here.

The following functions should help - feel free to use/modify of course. Worth checking what is omitted - some attributes of the spatial reference systems will be harmless if they don't match but that's up to you :)

import arcpy

def check_crs(dataset): """Return a coordinate reference system string

Get coordinate reference system of dataset
"""
crs = arcpy.Describe(dataset).spatialReference
return(crs)

def assert_crs_attribs(dataset1, dataset2, strict=False): """Returns Nothing

Asserts equality of all attributes of the provided geoprocessing spatial reference objects.
These are generated using arcpy.Describe(your_dataset).spatialReference.
Attributes of spatial reference object: https://pro.arcgis.com/en/pro-app/arcpy/classes/spatialreference.htm

dataset1 - a spatial dataset with projection info e.g. shp
dataset2 - a spatial dataset with projection info e.g. shp
strict - boolean - if True will compare every element (default: False)
"""
crs1=check_crs(dataset1)
crs2=check_crs(dataset2)

try:
    # Consider these
    assert(crs1.name==crs2.name) # The name of the spatial reference.
    assert(crs1.PCSCode==crs2.PCSCode) # The projected coordinate system code.1 
    assert(crs1.PCSName==crs2.PCSName) # The projected coordinate system name.1 
    assert(crs1.azimuth==crs2.azimuth) # The azimuth of a projected coordinate system.1 
    assert(crs1.centralMeridian==crs2.centralMeridian) # The central meridian of a projected coordinate system.1    
    assert(crs1.centralMeridianInDegrees==crs2.centralMeridianInDegrees) # The central meridian (Lambda0) of a projected coordinate system in degrees.1 
    assert(crs1.centralParallel==crs2.centralParallel) # The central parallel of a projected coordinate system.1
    assert(crs1.falseEasting==crs2.falseEasting) # The false easting of a projected coordinate system.1 
    assert(crs1.falseNorthing==crs2.falseNorthing) # The false northing of a projected coordinate system.1  
    assert(crs1.MFalseOriginAndUnits==crs2.MFalseOriginAndUnits) # The measure false origin and units.
    assert(crs1.MResolution==crs2.MResolution) # The measure resolution.
    assert(crs1.MTolerance==crs2.MTolerance) # The measure tolerance.
    assert(crs1.XYTolerance==crs2.XYTolerance) # The xy tolerance.
    assert(crs1.ZDomain==crs2.ZDomain) # The extent of the z domain.
    assert(crs1.ZFalseOriginAndUnits==crs2.ZFalseOriginAndUnits) # The z false origin and units.
    assert(crs1.factoryCode==crs2.factoryCode) # The factory code or well-known ID (WKID) of the spatial reference.
    assert(crs1.isHighPrecision==crs2.isHighPrecision) # Indicates whether the spatial reference has high precision set.
    assert(crs1.latitudeOf1st==crs2.latitudeOf1st) # The latitude of the first point of a projected coordinate system.1
    assert(crs1.latitudeOf2nd==crs2.latitudeOf2nd) # The latitude of the second point of a projected coordinate system.1    
    assert(crs1.latitudeOfOrigin==crs2.latitudeOfOrigin) # The latitude of origin of a projected coordinate system.1    
    assert(crs1.linearUnitCode==crs2.linearUnitCode) # The linear unit code.    
    assert(crs1.linearUnitName==crs2.linearUnitName) # The linear unit name.1
    assert(crs1.longitude==crs2.longitude) # The longitude value of this prime meridian.1
    assert(crs1.longitudeOf1st==crs2.longitudeOf1st) #The longitude of the first point of a projected coordinate system.1
    assert(crs1.longitudeOf2nd==crs2.longitudeOf2nd) # The longitude of the second point of a projected coordinate system.1
    assert(crs1.longitudeOfOrigin==crs2.longitudeOfOrigin) # The longitude of origin of a projected coordinate system.1
    assert(crs1.metersPerUnit==crs2.metersPerUnit) # The meters per linear unit.1
    assert(crs1.projectionCode==crs2.projectionCode) # The projection code.1
    assert(crs1.projectionName==crs2.projectionName) # The projection name.1
    assert(crs1.scaleFactor==crs2.scaleFactor) # The scale factor of a projected coordinate system.1
    assert(crs1.standardParallel1==crs2.standardParallel1) # The first parallel of a projected coordinate system.1
    assert(crs1.standardParallel2==crs2.standardParallel2) # The second parallel of a projected coordinate system.1
    assert(crs1.angularUnitCode==crs2.angularUnitCode) # The angular unit code.2
    assert(crs1.angularUnitName==crs2.angularUnitName) # The angular unit name.2
    assert(crs1.datumCode==crs2.datumCode) # The datum code.2
    assert(crs1.datumName==crs2.datumName) # The datum name.2
    assert(crs1.flattening==crs2.flattening) # The flattening ratio of this spheroid.2
    assert(crs1.longitude==crs2.longitude) # The longitude value of this prime meridian.2
    assert(crs1.primeMeridianCode==crs2.primeMeridianCode) # The prime meridian code.2

    ## Prob can be ignored
    if strict:
        assert(crs1.ZResolution==crs2.ZResolution) # The z resolution property.
        assert(crs1.ZTolerance==crs2.ZTolerance) # The z-tolerance property.
        assert(crs1.hasMPrecision==crs2.hasMPrecision) # Indicates whether m-value precision information has been defined.
        assert(crs1.hasXYPrecision==crs2.hasXYPrecision) # Indicates whether xy precision information has been defined.
        assert(crs1.hasZPrecision==crs2.hasZPrecision) # Indicates whether z-value precision information has been defined.
        assert(crs1.XYResolution==crs2.XYResolution) # The xy resolution.
        assert(crs1.domain==crs2.domain) # The extent of the xy domain.
        assert(crs1.MDomain==crs2.MDomain) # The extent of the measure domain.
        assert(crs1.remarks==crs2.remarks) # The comment string of the spatial reference.
        assert(crs1.type==crs2.type) # The type of the spatial reference. Geographic: A geographic coordinate system. Projected: A projected coordinate system.
        assert(crs1.usage==crs2.usage) # The usage notes.   
        assert(crs1.classification==crs2.classification) # The classification of a map projection.1 
        assert(crs1.GCSCode==crs2.GCSCode) # The geographic coordinate system code.2
        assert(crs1.GCSName==crs2.GCSName) # The geographic coordinate system name.2
        assert(crs1.primeMeridianName==crs2.primeMeridianName) # The prime meridian name.2
        assert(crs1.radiansPerUnit==crs2.radiansPerUnit) # The radians per angular unit.2
        assert(crs1.semiMajorAxis==crs2.semiMajorAxis) # The semi-major axis length of this spheroid.2
        assert(crs1.semiMinorAxis==crs2.semiMinorAxis) # The semi-minor axis length of this spheroid.2
        assert(crs1.spheroidCode==crs2.spheroidCode) # The spheroid code.2
        assert(crs1.spheroidName==crs2.spheroidName) # The spheroid name.2
    return(True)
except:
    output_message="CRS differs between datasets."#\ncrs1: %s\ncrs2 : %s" %(crs1.exportToString(), crs2.exportToString())
    print(output_message)
    return(False)
## Differs to the falseEasting and falseNorthingUnits are odd on occasion but false eastings and northings make sense
# crs.falseOriginAndUnits # The false origin and units.

## Not required
#crs.GCS # A projected coordinate system returns a SpatialReference object for the geographic coordinate system it is based on. A geographic crs.coordinate system returns the same SpatialReference.
#crs.SpatialReference
#crs.VCS # If the coordinate system has a vertical coordinate system, it returns a VCS object for the vertical coordinate system it is based on.
#crs.abbreviation # The abbreviated name of the spatial reference.
#crs.alias # The alias of the spatial reference.

Given the above, you might use them like:

dataset1="your_vector_1.shp"
dataset2="your_vector_2.shp"

assert_crs_attribs(dataset1, dataset2)

Given your use case, hopefully the asserts won't fail.

I integrate these functions now into many processes such as where I have a range of spatial datasets that I'm joining and want to eliminate any doubt that things were misaligned.

PolyGeo
  • 65,136
  • 29
  • 109
  • 338
ChrisWills
  • 135
  • 1
  • 7