[GIS] How to Transform Many Projections Into One With SQL Server

arcpycoordinate systemsql server

I know very little SQL Server, but am quite familiar with ArcGIS and Python. The data I want to transform is point data which is in a single table [ID], [X] [Y], [COORDSYS]. There are 5 possible coordinate systems. NAD83 UTM Zone 13 and 14, NAD27 UTM Zone 13 and 14, WGS84.

Right now I have a Python tool that accesses the table, transforms the points based on the coordinate system and outputs a SHP file that is in WGS84. It takes about 5 minutes when the user starts it which is less than ideal.

It seems it is possible to do this with SQL Server spatial functions so that the same thing can be done in the background and be ready for the user to just display in ArcGIS. But I'm not sure how to best change which transformation to use depending on a value in the table using SQL. I need to know what it's called that I'm trying to do so I can look it up.

Following is the code I'm using to implement this right now. I suppose that the part which might be taking the longest is creating Polygons out of each distinct set of points, rather than the transformation of each point into WGS84

def addCorner(self, x, y, z):
        pointgeo = arcpy.PointGeometry(arcpy.Point(x,y,z), self.origSR, True)

        try:
            # transforms the points to WGS84 as they are added to the Loop object
            projectedPointGeo = pointgeo.projectAs(self.wgs84SR,self.transformation)
        except Exception as e:
            print("Please ensure that you have all the necessary transformation tables."
                           "You are likely missing the 'sk83-98.gsb' in your ArcGIS installation."
                           "Anyone in the GIS group can help.")

        # Adds each point as an arcpy.Point with Z values
        self.__corners.add(projectedPointGeo.firstPoint)

def CreateLoopsFromAcquire(loopsources, curdir, statusbar):
    if not curdir:
        if not os.path.exists(r"C:\temp"):
            os.mkdir(r"C:\temp")
        curdir = r"C:\temp"

    shapefiles = {}
    for loopsource in loopsources.values():
        print("-- From {0}".format(loopsource.tablename))
        cursor = arcpy.da.SearchCursor(loopsource.tablename,
                                        [loopsource.surveyfield, #0
                                        loopsource.corner_idfield, #1
                                        loopsource.coordtypefield, #2
                                        loopsource.xfield, #3
                                        loopsource.yfield, #4
                                        loopsource.zfield]) #5
        loops = {}

        # Gather Loops (group of corners) into Dict
        i = 0
        for row in cursor: # each row is a loop corner
            loopid = "{0}_{1}".format(row[0],row[2])
            if not loops.has_key(loopid):
                loops[loopid] = Sources.Loop(row[0], row[2])
                i+=1
            loops[loopid].addCorner(row[3],row[4],row[5])
        del cursor

        print("Found {0} loops".format(len(loops.items())))

        newfilepath = os.path.join(curdir, loopsource.outputname)
        if os.path.exists(loopsource.outputname):
            arcpy.Delete_management(newfilepath)

        newfc = arcpy.CreateFeatureclass_management(curdir, loopsource.outputname, "POLYGON", loopsource.templatepath, "SAME_AS_TEMPLATE", "SAME_AS_TEMPLATE", loopsource.templatepath)
        cursor = arcpy.da.InsertCursor(newfc,["SHAPE@","Lp_Survey"])

        for loop in loops.values():
            cursor.insertRow([loop.getPolygonOfLoop(),loop.uniquename])
            # statusbar.StepProgressBar()
        del cursor

        # return a shapefile object so that symbology can be made
        shpobjectname = os.path.splitext(loopsource.outputname)[0]
        shapefileobject = Sources.ShapeFile(
                                    newfilepath,
                                    shpobjectname,
                                    display_expression="Lp_Survey",
                                    symbology=loopsource.symbology)

        shapefiles[shpobjectname] = shapefileobject

    return shapefiles

Best Answer

SQL Server spatial functions implementation is in its current state not able to reproject or transform between different coordinate reference systems. The spatial data types are nothing more than a CLR (Coomon Language Runtime) data type implementation of the OGC Simple Feature Access specification plus some extended functions plus the spatial index. You can find its documentation for the geometry data type here for example.

The only way to bring reprojection capability to a SQL Server database is to write your own SQL Server user defined functions (the CLR way, which is another word for .NET). So basically every .NET library which is able to reproject can be integrated into SQL server by writing some glue code between the library and the SQL Server function call.

You can do this with ArcObjects .NET Edition or with dotSpatial for example. We have done this with dotSpatial in a Project, unfortunately our implementation is only able to reproject between two specific austrian CRS as this was the only Thing we needed to reach our Project Goal.

Related Question