Tell me more ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.
up vote 1 down vote favorite

I am retrieving thousands of geometry lines in WKT format from a Microsoft SQL Server database. I would like to plot all of those lines in one map, but I am having some difficulty with the code. How should I modify the following code to reach my goal?

library(RODBC)

png(file="examplex.png", width=600, height=480)
setwd("C:/ArcR")

con <- odbcDriverConnect('driver={SQL Server};server=SERVERNAME;database=DBNAME;trusted_connection=true')

objects_1 <- sqlQuery(con, 'SELECT TOP (1) Shape.STAsText() as ShapeWKT FROM TABLENAME ;')

things <- vector("list", 1)

z = 0
for(line in objects_1$ShapeWKT)
{
  {
    things[[z+1]]<-readWKT(line)
  }
  z = z + 1
}

plot(things[[1]]) 

dev.off()

Yes, the code does refer to only one returned record, but I wanted to start the question with code in working order, for those who may not know the answer, but are still capable of testing possibilities.

share|improve this question
1  
1. What difficulties are you having with this code as you have written it here? Can't test without the equivalent database. 2. This will be very slow for large numbers of objects if you add (z+1) to the vector each time - set the size initially with things <- vector("list", nrow(objects_1)) – Simbamangu Apr 19 at 17:20
The code as presented works well for only one record --plot(things[[x]])-- I would really like to use --Select Top(10)--, and --plot(things)--. Thank you for suggesting the list size be identified at creation. – ike Apr 19 at 17:37

2 Answers

up vote 1 down vote accepted

Since readWKT won't read more than one object at a time, you'll still have to iterate over each row of the dataframe returned from sqlQuery, but you can do it all from within R instead of system calls to ogr2ogr. Here's an example using SpatialPoints:

library(RODBC)
library(rgeos)

con <- odbcDriverConnect('driver={SQL Server}; server=SERVERNAME; database=DBNAME; trusted_connection=true')
objects_1 <- sqlQuery(con, 'SELECT Shape.STAsText() AS ShapeWKT, OBJECTID FROM TABLE;')

Set up a SpatialPointsDataFrame:

point.sp <- SpatialPointsDataFrame(readWKT(objects_1$ShapeWKT[1]), data=data.frame(OBJECTID=objects_1$OBJECTID[1]))

Iterate through the remaining rows.

for (n in 2:length(objects_1$OBJECTID)) {
  point.sp <- rbind(point.sp, 
                    SpatialPointsDataFrame(readWKT(objects_1$ShapeWKT[n]), 
                    data.frame(OBJECTID=objects_1$OBJECTID[n])))
}

Lacking your dataset, I ran it with some sample code:

obj_wkt <- c("POINT (30 10)", "POINT (30 11)", "POINT (31 11)", "POINT (31 10)")
obj_id <- 1:4
objects_1 <- data.frame(ShapeWKT=obj_wkt, OBJECTID=obj_id)

plot(point.sp, axes = 1)

enter image description here

You should set up the initial SpatialPointsDataFrame to be the length of objects_1 (with zeroes) or this will slow down the iteration even more.

share|improve this answer
This tosses the error : (Error in readWKT(objects_1$WKT) : WKT must have length 1). So this works for single objects, and if I string single reads and plots together it works for several objects. I need this to work for thousands/millions/ (possibly)billions of objects. – ike Apr 19 at 14:50
And yes, it is unlikely that a temporally/characteristically invariable map will even be capable of showing billions of objects without being ridiculously large. I am more interested in exposing the difference between simulations of 50 million objects across a minimum of 30 separate models through a time-lapse interface. – ike Apr 19 at 15:25
That is quite a lot of memory though, so I'm ok with a subset of several thousand for now. – ike Apr 19 at 15:35
1  
OK - I've never had a chance to try that with a dB connection, didn't realise it worked one at a time (and that's why you were iterating over nrows in your question). – Simbamangu Apr 19 at 17:16
This is pretty good. Your code is doing my points perfectly, but my lines aren't doing so well. For now though, this is probably the best answer. – ike Apr 22 at 16:01
show 1 more comment
up vote 0 down vote

This answer requires shelling out ogr2ogr. I'm sure there has to be a better way, but for now, this one works.

library(RODBC)

xmin <- xxxxxxx  #insert your number here
xmax <- xxxxxxx  #insert your number here
ymin <- xxxxxx   #insert your number here
ymax <- xxxxxx   #insert your number here

png(file="examplex.png", width=600, height=480)
setwd("C:/ArcR")

con <- odbcDriverConnect('driver={SQL     Server};server=SERVERNAME;database=DBNAME;trusted_connection=true')

objects_1 <- sqlQuery(con, 'SELECT Shape.STAsText() AS ShapeWKT, OBJECTID FROM TABLE;')

cat("WKT;OBJECTID",file="outfile.csv",sep="\n")

for(i in 1:nrow(objects_1)) {
  row <- objects_1[i,]

  cat('"',file="outfile.csv",append=TRUE)
  cat(toString(row$ShapeWKT),file="outfile.csv",append=TRUE)
  cat('";',file="outfile.csv",append=TRUE)
  cat(row$OBJECTID,file="outfile.csv",sep="\n",append=TRUE)
}

shell('ogr2ogr -f "ESRI Shapefile" mydata.shp "C:\\ArcR\\Outfile.csv"')

Inspections.mp <- readShapeLines("mydata")

#I put a shapefile called 'boundary' in the working directory.  Its not important, but its how I made my results.
border.mp <- readShapePoly("Boundary")

plot.new()
par(mar=c(0,0,0,0))
plot(border.mp, col="#000000FF", xlim = c(xmin, xmax), ylim = c(ymin, ymax), border="#000000FF")
lines(Inspections.mp, col="blue",lwd=0.25)

dev.off()

I'm pretty sure I'm misusing 'cat' here, but um... eh.

share|improve this answer
What about readWKT from the rgeos package? Won't that read your ShapeWKT data? – Simbamangu Apr 19 at 4:09
ReadWKT will only read ShapeWKT data one object at a time, and does not seem to remain linked with the associated table data that comes with it. – ike Apr 19 at 14:57

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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