Python arcpy.SpatialReference() Examples

def create_feature_class(out_path,
                         out_name,
                         geom_type,
                         wkid,
                         fields,
                         objectIdField):
    """ creates a feature class in a given gdb or folder """
    if arcpyFound == False:
        raise Exception("ArcPy is required to use this function")
    arcpy.env.overwriteOutput = True
    field_names = []
    fc =arcpy.CreateFeatureclass_management(out_path=out_path,
                                            out_name=out_name,
                                            geometry_type=lookUpGeometry(geom_type),
                                            spatial_reference=arcpy.SpatialReference(wkid))[0]
    for field in fields:
        if field['name'] != objectIdField:
            field_names.append(field['name'])
            arcpy.AddField_management(out_path + os.sep + out_name,
                                      field['name'],
                                      lookUpFieldType(field['type']))
    return fc, field_names
#---------------------------------------------------------------------- 

def to_points(tbl, out_fc, xcol, ycol, sr, zcol='#', w=''):
    """Convert table to point feature class, return path to the feature class.

    Required:
    tbl -- input table or table view
    out_fc -- path to output feature class
    xcol -- name of a column in tbl that stores x coordinates
    ycol -- name of a column in tbl that stores y coordinates
    sr -- spatial reference for out_fc
        sr can be either arcpy.SpatialReference object or a well known id as int

    Optional:
    zcol -- name of a column in tbl that stores y coordinates, default is '#'
    w -- where clause to limit the rows of tbl considered, default is ''

    Example:
    >>> t = 'c:\\foo\\bar.shp'
    >>> o = 'c:\\foo\\bar_pts.shp'
    >>> table_to_points(t, o, "XC", "YC", 4326, zcol='#', w='"FID" < 10')
    >>> table_to_points(t, o, "XC", "YC", arcpy.SpatialReference(27700))
    >>> table_to_points(t, o, "XC", "YC", arcpy.describe(tbl).spatialReference)
    """
    lrnm = tstamp('lr', '%m%d%H%M%S', '')
    if type(sr) != arcpy.SpatialReference:
        sr = arcpy.SpatialReference(sr)
    lr = arcpy.MakeXYEventLayer_management(tbl, xcol, ycol, lrnm, sr, zcol).getOutput(0)
    if str(w) not in ('', '*'):
        arcpy.SelectLayerByAttribute_management(lr, "NEW_SELECTION", w)
    out_fc = arcpy.CopyFeatures_management(lr, out_fc).getOutput(0)
    dlt(lr)
    return (arcpy.Describe(out_fc).catalogPath) 

def createPolygon(self, lat, lon, extent, out_polygons, scratchWorkspace):
        """Create a Thiessen polygon feature class from numpy.ndarray lat and lon
           Each polygon represents the area described by the center point
        """
        buffer = 2 * max(abs(lat[0]-lat[1]),abs(lon[0] - lon[1]))
        # Extract the lat and lon within buffered extent (buffer with 2* interval degree)
        lat0 = lat[(lat >= (extent.YMin - buffer)) & (lat <= (extent.YMax + buffer))]
        lon0 = lon[(lon >= (extent.XMin - buffer)) & (lon <= (extent.XMax + buffer))]
        # Spatial reference: GCS_WGS_1984
        sr = arcpy.SpatialReference(4326)

        # Create a list of geographic coordinate pairs
        pointGeometryList = []
        for i in range(len(lon0)):
            for j in range(len(lat0)):
                point = arcpy.Point()
                point.X = float(lon0[i])
                point.Y = float(lat0[j])
                pointGeometry = arcpy.PointGeometry(point, sr)
                pointGeometryList.append(pointGeometry)

        # Create a point feature class with longitude in Point_X, latitude in Point_Y
        out_points = os.path.join(scratchWorkspace, 'points_subset')
        result2 = arcpy.CopyFeatures_management(pointGeometryList, out_points)
        out_points = result2.getOutput(0)
        arcpy.AddGeometryAttributes_management(out_points, 'POINT_X_Y_Z_M')

        # Create Thiessen polygon based on the point feature
        result3 = arcpy.CreateThiessenPolygons_analysis(out_points, out_polygons, 'ALL')
        out_polygons = result3.getOutput(0)

        return out_points, out_polygons 

def execute(in_datasets, out_fc):
    # use gcs as output sr since all extents will fit in it
    out_sr = arcpy.SpatialReference("WGS 1984")

    in_datasets = in_datasets.split(";")

    arcpy.CreateFeatureclass_management(os.path.dirname(out_fc),
                                        os.path.basename(out_fc),
                                        "POLYGON",
                                        spatial_reference=out_sr)

    arcpy.AddField_management(out_fc, "dataset", "TEXT", 400)

    # add each dataset's extent & the dataset's name to the output
    with arcpy.da.InsertCursor(out_fc, ("SHAPE@", "dataset")) as cur:

        for i in in_datasets:
            d = arcpy.Describe(i)
            ex = d.Extent
            pts = arcpy.Array([arcpy.Point(ex.XMin, ex.YMin),
                              arcpy.Point(ex.XMin, ex.YMax),
                              arcpy.Point(ex.XMax, ex.YMax),
                              arcpy.Point(ex.XMax, ex.YMin),
                              arcpy.Point(ex.XMin, ex.YMin),])

            geom = arcpy.Polygon(pts,  d.SpatialReference)

            if d.SpatialReference != out_sr:
                geom = geom.projectAs(out_sr)
            cur.insertRow([geom, d.CatalogPath]) 

def maskCoastlineConflicts(prod_db, desktop_fldr):
    arcpy.AddMessage("\tMasking coastline and bridges")
    # Subtype field used in where clause to access bridges in CulturalFeaturesA
    subtype_fld = arcpy.AddFieldDelimiters(prod_db, "FCSubtype")
    # Get subtype of Bridge
    bridge = "5"
    # Define spatial reference
    sr = arcpy.SpatialReference(4326)

    # Get CoastlineL and CulturalFeaturesA layers
    coastlinel_fc = getFC(prod_db, "CoastlineL", NAUT_FDS)
    culturalfeaturesa_fc = getFC(prod_db, "CulturalFeaturesA", NAUT_FDS)

    # Only continue if CoastlineL and CulturalFeaturesA layers are in the TOC
    if coastlinel_fc != "" and culturalfeaturesa_fc != "":
        # Make feature layer form CoastlineL
        arcpy.MakeFeatureLayer_management(coastlinel_fc, "coastlinel_lyr")
        # Make feature layer of bridge features
        where = subtype_fld + " = " + bridge
        arcpy.MakeFeatureLayer_management(culturalfeaturesa_fc, "bridges", where)
        # Check if there are any bridge features in the layer
        if int(arcpy.GetCount_management("bridges").getOutput(0)) > 0:
            # Run Intersecting Layers Mask GP tool to create mask poly where coastline intersect bridges
            mask_fc = os.path.join(prod_db, CARTO_FDS, "MASK_CoastlineL")
            arcpy.IntersectingLayersMasks_cartography("bridges", "coastlinel_lyr", mask_fc, REF_SCALE, sr, "0.01 POINTS")

    return 

def project_coordinates(xys, in_sr, out_sr, datum_transformation=None):
    """Project list of coordinate pairs (or triplets).
        xys -- list of coordinate pairs or triplets to project one by one
        in_sr -- input spatial reference, wkid, prj file, etc.
        out_sr -- output spatial reference, wkid, prj file, etc.
        datum_transformation=None -- datum transformation to use
            if in_sr and out_sr are defined on different datums,
            defining appropriate datum_transformation is necessary
            in order to obtain correct results!
            (hint: use arcpy.ListTransformations to list valid transformations)

    Example:
    >>> dtt = 'TM65_To_WGS_1984_2 + OSGB_1936_To_WGS_1984_NGA_7PAR'
    >>> coordinates = [(240600.0, 375800.0), (245900.0, 372200.0)]
    >>> project_coordinates(coordinates, 29902, 27700, dtt)
    """

    if not type(in_sr) is arcpy.SpatialReference:
        in_sr = arcpy.SpatialReference(in_sr)
    if not type(out_sr) is arcpy.SpatialReference:
        out_sr = arcpy.SpatialReference(out_sr)

    xyspr = []
    for xy in xys:
        pt = arcpy.Point(*xy)
        hasz = True if pt.Z is not None else False
        ptgeo = arcpy.PointGeometry(pt, in_sr)
        ptgeopr = ptgeo.projectAs(out_sr, datum_transformation)
        ptpr = ptgeopr.firstPoint
        if hasz:
            xypr = (ptpr.X, ptpr.Y, ptpr.Z)
        else:
            xypr = (ptpr.X, ptpr.Y)
        xyspr.append(xypr)

    return xyspr 

def asShape(self):
        """Returns JSON as arcpy.Geometry() object."""
        if self.geometryType != ESRI_ENVELOPE:
            return arcpy.AsShape(self.json, True)
        else:
            ar = arcpy.Array([
                arcpy.Point(self.json[XMIN], self.json[YMAX]),
                arcpy.Point(self.json[XMAX], self.json[YMAX]),
                arcpy.Point(self.json[XMAX], self.json[YMIN]),
                arcpy.Point(self.json[XMIN], self.json[YMIN])
            ])
            return arcpy.Polygon(ar, arcpy.SpatialReference(self.spatialReference)) 

def defineFeatureProjection(fileInfo):
    rootPath=list(fileInfo.keys())  #待读取数据文件的根目录
#    print(rootPath)
    dataName=flatten_lst(list(fileInfo.values()))  #待读取数据文件的文件名列表
    outCS=arcpy.SpatialReference('WGS 1984 UTM Zone 49N') #定义投影
    for fName in dataName: 
        in_features=os.path.join(rootPath[0],fName)
        out_raster=fName+"_prj.shp"
        arcpy.Project_management(in_features, out_raster, outCS) 

def defineFeatureProjection(fileInfo):
    rootPath=list(fileInfo.keys())  #待读取数据文件的根目录
#    print(rootPath)
    dataName=flatten_lst(list(fileInfo.values()))  #待读取数据文件的文件名列表
    # outCS=arcpy.SpatialReference('WGS 1984 UTM Zone 49N') #定义投影 
    outCS=arcpy.SpatialReference('WGS 1984 UTM zone 16N') #定义投影 
    for fName in dataName: 
        in_features=os.path.join(rootPath[0],fName)
        out_raster=fName+"_prj.shp"
        arcpy.Project_management(in_features, out_raster, outCS) 

def compareDatum(fc):
    # Return True if fc datum is either WGS84 or NAD83

    try:
        # Create Spatial Reference of the input fc. It must first be converted in to string in ArcGIS10
        # otherwise .find will not work.
        fcSpatialRef = str(arcpy.CreateSpatialReference_management("#",fc,"#","#","#","#"))
        FCdatum_start = fcSpatialRef.find("DATUM") + 7
        FCdatum_stop = fcSpatialRef.find(",", FCdatum_start) - 1
        fc_datum = fcSpatialRef[FCdatum_start:FCdatum_stop]

        # Create the GCS WGS84 spatial reference and datum name using the factory code
        WGS84_sr = arcpy.SpatialReference(4326)
        WGS84_datum = WGS84_sr.datumName

        NAD83_datum = "D_North_American_1983"

        # Input datum is either WGS84 or NAD83; return true
        if fc_datum == WGS84_datum or fc_datum == NAD83_datum:
            del fcSpatialRef, FCdatum_start, FCdatum_stop,fc_datum,WGS84_sr,WGS84_datum,NAD83_datum
            return True

        # Input Datum is some other Datum; return false
        else:
            del fcSpatialRef, FCdatum_start, FCdatum_stop,fc_datum,WGS84_sr,WGS84_datum,NAD83_datum
            return False

    except arcpy.ExecuteError:
        AddMsgAndPrint(arcpy.GetMessages(2),2)
        return False

    except:
        print_exception()
        return False

## =============================================================================================================== 

def compareDatum(fc):
    # Return True if fc datum is either WGS84 or NAD83

    try:
        # Create Spatial Reference of the input fc. It must first be converted in to string in ArcGIS10
        # otherwise .find will not work.
        fcSpatialRef = str(arcpy.CreateSpatialReference_management("#",fc,"#","#","#","#"))
        FCdatum_start = fcSpatialRef.find("DATUM") + 7
        FCdatum_stop = fcSpatialRef.find(",", FCdatum_start) - 1
        fc_datum = fcSpatialRef[FCdatum_start:FCdatum_stop]

        # Create the GCS WGS84 spatial reference and datum name using the factory code
        WGS84_sr = arcpy.SpatialReference(4326)
        WGS84_datum = WGS84_sr.datumName

        NAD83_datum = "D_North_American_1983"

        # Input datum is either WGS84 or NAD83; return true
        if fc_datum == WGS84_datum or fc_datum == NAD83_datum:
            del fcSpatialRef, FCdatum_start, FCdatum_stop,fc_datum,WGS84_sr,WGS84_datum,NAD83_datum
            return True

        # Input Datum is some other Datum; return false
        else:
            del fcSpatialRef, FCdatum_start, FCdatum_stop,fc_datum,WGS84_sr,WGS84_datum,NAD83_datum
            return False

    except arcpy.ExecuteError:
        AddMsgAndPrint(arcpy.GetMessages(2),2)
        return False

    except:
        print_exception()
        return False

## =============================================================================================================== 

def create_empty_schema(feature_set, out_fc):
    """Creates empty schema in a feature class.

    Args:
        feature_set: Input feature set.
        out_fc: Output feature class path.

    Returns:
        A feature class.
    """

    # make copy of feature set
    fs = feature_set.getEmptyCopy()
    try:
        try:
            # this tool has been very buggy in the past :(
            tmp = fs.dump(tmp_json_file(), indent=None)
            arcpy.conversion.JSONToFeatures(tmp, out_fc)
        except:
            # this isn't much better..
            gp = arcpy.geoprocessing._base.Geoprocessor()

            # create arcpy.FeatureSet from raw JSON string
            arcpy_fs = gp.fromEsriJson(fs.dumps(indent=None))
            arcpy_fs.save(out_fc)

    except:
        # manually add records with insert cursor, this is SLOW!
        print('arcpy conversion failed, manually writing features...')
        outSR = arcpy.SpatialReference(fs.getSR())
        path, fc_name = os.path.split(out_fc)
        g_type = G_DICT.get(fs.geometryType, '').upper()
        arcpy.management.CreateFeatureclass(path, fc_name, g_type,
                                        spatial_reference=outSR)

        # add all fields
        for field in fs.fields:
            if field.type not in [OID, SHAPE] + list(SKIP_FIELDS.keys()):
                if '.' in field.name:
                    if 'shape.' not in field.name.lower():
                        field_name = field.name.split('.')[-1] #for weird SDE fields with periods
                    else:
                        field_name = '_'.join([f.title() for f in field.name.split('.')]) #keep geometry calcs if shapefile
                else:
                    field_name = field.name


                # need to filter even more as SDE sometimes yields weird field names...sigh
                restricted = ('fid', 'shape', 'objectid')
                if (not any(['shape_' in field.name.lower(),
                            'shape.' in field.name.lower(),
                            '(shape)' in field.name.lower()]) \
                            or isShp) and field.name.lower() not in restricted:
                    field_length = field.length if hasattr(field, 'length') else None
                    arcpy.management.AddField(out_fc, field_name, FTYPES[field.type],
                                                field_length=field_length,
                                                field_alias=field.alias)
        return out_fc 

def compareDatum(fc):
    # Return True if fc datum is either WGS84 or NAD83

    try:
        # Create Spatial Reference of the input fc. It must first be converted in to string in ArcGIS10
        # otherwise .find will not work.
        fcSpatialRef = str(arcpy.CreateSpatialReference_management("#",fc,"#","#","#","#"))
        FCdatum_start = fcSpatialRef.find("DATUM") + 7
        FCdatum_stop = fcSpatialRef.find(",", FCdatum_start) - 1
        fc_datum = fcSpatialRef[FCdatum_start:FCdatum_stop]

        # Create the GCS WGS84 spatial reference and datum name using the factory code
        WGS84_sr = arcpy.SpatialReference(4326)
        WGS84_datum = WGS84_sr.datumName

        NAD83_datum = "D_North_American_1983"

        # Input datum is either WGS84 or NAD83; return true
        if fc_datum == WGS84_datum or fc_datum == NAD83_datum:
            del fcSpatialRef
            del FCdatum_start
            del FCdatum_stop
            del fc_datum
            del WGS84_sr
            del WGS84_datum
            del NAD83_datum

            return True

        # Input Datum is some other Datum; return false
        else:
            del fcSpatialRef
            del FCdatum_start
            del FCdatum_stop
            del fc_datum
            del WGS84_sr
            del WGS84_datum
            del NAD83_datum

            return False

    except arcpy.ExecuteError:
        AddMsgAndPrint(arcpy.GetMessages(2),2)
        return False

    except:
        print_exception()
        return False

## =============================================================================================================== 

def compareDatum(fc):
    # Return True if fc datum is either WGS84 or NAD83

    try:
        # Create Spatial Reference of the input fc. It must first be converted in to string in ArcGIS10
        # otherwise .find will not work.
        fcSpatialRef = str(arcpy.CreateSpatialReference_management("#",fc,"#","#","#","#"))
        FCdatum_start = fcSpatialRef.find("DATUM") + 7
        FCdatum_stop = fcSpatialRef.find(",", FCdatum_start) - 1
        fc_datum = fcSpatialRef[FCdatum_start:FCdatum_stop]

        # Create the GCS WGS84 spatial reference and datum name using the factory code
        WGS84_sr = arcpy.SpatialReference(4326)
        WGS84_datum = WGS84_sr.datumName

        NAD83_datum = "D_North_American_1983"

        # Input datum is either WGS84 or NAD83; return true
        if fc_datum == WGS84_datum or fc_datum == NAD83_datum:
            del fcSpatialRef
            del FCdatum_start
            del FCdatum_stop
            del fc_datum
            del WGS84_sr
            del WGS84_datum
            del NAD83_datum

            return True

        # Input Datum is some other Datum; return false
        else:
            del fcSpatialRef
            del FCdatum_start
            del FCdatum_stop
            del fc_datum
            del WGS84_sr
            del WGS84_datum
            del NAD83_datum

            return False

    except arcpy.ExecuteError:
        AddMsgAndPrint(arcpy.GetMessages(2),2)
        return False

    except:
        print_exception()
        return False

## =============================================================================================================== 

def createOutputFC(dictOfFields,preFix="Temp",shape="POLYGON"):
    """ This function will creat an empty polygon feature class within the scratch FGDB.  The feature class will be in WGS84
        and will have have 2 fields created: mukey and mupolygonkey.  The feature class name will have a prefix of
        'nasisProject.' This feature class is used by the 'requestGeometryByMUKEY' function to write the polygons associated with a NASIS
        project.

        fieldDict ={field:(fieldType,fieldLength,alias)
        fieldDict = {"mukey":("TEXT",30,"Mapunit Key"),"mupolygonkey":("TEXT","",30)}

        Return the new feature class  including the path.  Return False if error ocurred."""

    try:
        epsgWGS84 = 4326 # EPSG code for: GCS-WGS-1984
        outputCS = arcpy.SpatialReference(epsgWGS84)

        newFC = arcpy.CreateScratchName(preFix, workspace=arcpy.env.scratchGDB,data_type="FeatureClass")

        # Create empty polygon featureclass with coordinate system that matches AOI.
        arcpy.CreateFeatureclass_management(env.scratchGDB, os.path.basename(newFC), shape, "", "DISABLED", "DISABLED", outputCS)

        for field,params in dictOfFields.iteritems():
            try:
                fldLength = params[1]
                fldAlias = params[2]
            except:
                fldLength = 0
                pass

            arcpy.AddField_management(newFC,field,params[0],"#","#",fldLength,fldAlias)

##            if len(params[1]) > 0:
##                expression = "\"" + params[1] + "\""
##                arcpy.CalculateField_management(helYesNo,field,expression,"VB")

##        arcpy.AddField_management(nasisProjectFC,"mukey", "TEXT", "", "", "30")
##        arcpy.AddField_management(nasisProjectFC,"mupolygonkey", "TEXT", "", "", "30")   # for outputShp

        if not arcpy.Exists(newFC):
            AddMsgAndPrint("\tFailed to create scratch " + newFC + " Feature Class",2)
            return False

        return newFC

    except:
        errorMsg()
        AddMsgAndPrint("\tFailed to create scratch " + newFC + " Feature Class",2)
        return False

## ===================================================================================