Read SQL Server Geometry datatype in R

Question

I have two datasets that come from ArcGis, stored in SQL Server. I have a column called - Shape which is a geometry data type and the value looks like this:

0x110F0000010407000000DCD781517FD76F4140FAED5B2EFF46C190C2F53C9FD76F41A0703D5242FF46C1C8BAB873A1D76F416888632504FF46C1B0BFECF4A3D76F41801D38F7C5FE46C19031770984D76F41A080264AB5FE46C114D0448881D76F4100B37B6AF2FE46C1DCD781517FD76F4140FAED5B2EFF46C101000000020000000001000000FFFFFFFF0000000003

After executing this query for my Polygon dataset: SELECT Shape.STGeometryType() AS GeometryType FROM dbo.POLY

I get Polygon.

I have another dataset that is a POINT dataset which also has the same column Shape.

I use this query to do an Intersect between the two datasets:

update ABC 
 set ABC.ID = mb.ID
 from STG ABC 
 join POLY mb with (index([SPIDX_POLY_Shape]))
 on ABC.Shape.STIntersects(mb.Shape) = 1

My question is : How can I read this data into R, specially the Shape column, when I connect my SQL Server database.
And my second question is which function will help me execute the above spatial intersect query? I am good with R, but bad with Spatial Analysis. I am not looking for you all to write me a code, but help me understand some concepts.

Answer 1

There is functions wkb::readWKB and rgeos::readWKT that return Spatial*DataFrames in R.

The general idea is to cast your geometry to one of these types and build the objects after reading the right type with SQL. You can connect to SQL Server with RODBC, or with this package: https://cran.rstudio.com/web/packages/RSQLServer/index.html

(I haven't used RSQLServer. )

Your database will have converter functions to cast data to a useable type. In Manifold GIS I use

CGeomWKB(Geom(ID)) 

inside the SQL query which returns a column of WKB type, then I read that with wkb::readWKB. SQL Server would have similar facilities.

You can see this in action in this function: https://github.com/mdsumner/manifoldr/blob/master/R/manifoldr.r

Note that R's Spatial Classes don't support mixing of types, like polygons, lines and points in the same data set so you'll need to filter out the ones you want to build individually if you see layers with this. Also note that WKB and WKT don't include the CRS, though Manifold's Geom does (in their format), and SQL Server may do similar. My Manifold code gets that CRS via a different mechanism.

Ultimately you could parse the native (or Shape) SQL Server Geom rather than casting to these other types, perhaps with readBin but that's another level of complexity.

• CuriousBeing: Thanks. I am using the rgeos and sp package to create my SpatialPointDataFrame. But, I have run into another problem and that is my loop keeps running forever. I've got 2 million pts nd its been sitting for 3.5 hours on it. My R statements: point.sp <- SpatialPointsDataFrame(readWKT(df$WKT[1]),data=data.frame(OBJECTID=df$OBJECTID[1], LOT_NO=df$LOT_NO[1])) for (n in 2:as.integer(cnt)) { point.sp <- rbind(point.sp, SpatialPointsDataFrame(readWKT(df$WKT[n]), data.frame(OBJECTID=df$OBJECTID[n], LOT_NO=df$LOT_NO[n]))) }
• mdsumner: readWKT is vectorized, there's no need to loop - and especially not by "growing" the object incrementally - that is hopelessly inefficient. If you have another question, don't post it in comments, but probably all you need is to remove the index: SpatialPointsDataFrame(readWKT(df$WKT),data=data.frame(OBJECTID=df$OBJECTID, LOT_NO=df$LOT_NO)) - might need to handle the rownames too.
• CuriousBeing: I've posted a new question. Thank you and sorry my bad. [gis.stackexchange.com/questions/182497/…
• mdsumner: Cripes, sorry - readWKT is not vectorized. I don't know why that is. I'll see if can be. There'll be a much faster way to do it, and via WKB would be much better.