1

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
4
  • How many million point features in this table? Maybe you should look at the efficiency of the Python code.
    – Vince
    Commented Mar 6, 2015 at 18:24
  • You are right. There are about 40,000 points for now, but I expect it to get up to a few hundred thousand. I'm using arcpy.PointGeometry.projectAs for each point as I'm going through the table. It was a simple implementation until we figured out a better way to structure the database in the backend. I am assuming that having a database table ready to go in the right projection will still be faster than finding another way to implement the transformation on the user's side. Is this true?
    – TSJ
    Commented Mar 6, 2015 at 19:05
  • 1
    Are you using a dictionary to store the source coordinate references? What version of ArcGIS are you using? Are you using DA cursors? I doubt you'd see a performance difference between the most efficient implementations with FGDB and EGDB targets.
    – Vince
    Commented Mar 6, 2015 at 20:52
  • I'll update the answer with the code I'm using right now.
    – TSJ
    Commented Mar 6, 2015 at 20:55

1 Answer 1

3

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.

1
  • This makes sense. I am in another job now, sorry I won't be able to take this further. Thanks for the answer
    – TSJ
    Commented Apr 27, 2016 at 0:37

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.