Search code examples
rpostgresqlpostgisspatialogr

How to import OSM polygons data from PostGIS to R?


I am new to the world of spatial analysis using R. Using this link I have downloaded OSM data in .osm.pbf format. Then I used osm2pgsql tool to get data in PostgreSQL (PostGIS extension). Now I have several tables in my database and I want to access the polygons table in R and then perform spatial analysis on the polygon vector data. I have been searching around allot but am not able to import required data in R. I found this tutorial quite similar to what I am looking for but its in Python. I want to access polygon data from PostGIS using R.

Therefore, essentially I would like to know interaction of R with PostGIS. Can anybody recommend me any book on this topic? Since I couldn't find a blog or tutorial so far that works for me on my Windows 10 64-bit machine.

Thanks for your time and looking forward for the suggestions.


Solution

  • I have still not found a way to get required data form PostGIS using rgdal package available in R. Probably it is because of my OS issues. (I am not exactly sure as I am not an expert). But I have found an alternative to rgdal and it has done exactly what I wanted it to do. The code is as following:

    library(RPostgreSQL)
    library(rgeos)
    library(sp)
    
    # Load data from the PostGIS server
    conn = dbConnect(
      dbDriver("PostgreSQL"), dbname="dbNAME", host="localhost", port=5432, 
      user="username", password="pw"
    )
    
    strSQL = "SELECT osm_id, name, area, highway, railway, place, ST_AsText(way) AS wkt_geometry FROM table"
    df = dbGetQuery(conn, strSQL)
    
    #Geomtery column as R list
    geo_col = df$wkt_geometry
    
    polygon_list = suppressWarnings(lapply(geo_col, function(x){
    x <- gsub("POLYGON\\(\\(", "", x)
    x <- gsub("\\)", "", x)
    x <- strsplit(x, ",")[[1]]
    
    #Now each polygon has been parsed by removing POLYGON(( from the start and )) from the end
    #Now for each POLYGON its xValues and yValues are to be extracted to for Polygon object
    xy <- strsplit(x, " ")
    
    v_xy = suppressWarnings(sapply(xy, function(p){        
      xValue = p[1]
      yValue = p[2]
      vec = c(xValue, yValue)
    }))
    
    #Now we have all x values in first column of v_xy and all y values in second column of v_xy
    #Let us make the Polygon object now
    p_xvalues = as.numeric(v_xy[1, ])
    p_yvalues = as.numeric(v_xy[2, ])
    p_object <- Polygon(cbind(p_xvalues, p_yvalues))      
    }))
    
    #Now we have all of the polygons in polygon object format
    #Let us join it with main data frame, i.e. df
    df$object_polygon <- polygon_list
    #View(df)
    
    #Now Let us form SpatialPolygons() object out of it
    Ps_list = list()
    for (i in seq(nrow(df))) {
      Ps_list[[i]] <- Polygons(polygon_list[i], ID=df[i,][1])
    }
    SPs = SpatialPolygons(Ps_list)
    
    #Now FINALY its the time to form SpatialPolygonsDataFrame
    row.names(df) = df$osm_id
    SPDF = SpatialPolygonsDataFrame(Sr = SPs, data = df[, 1:6], match.ID = TRUE) 
    

    Therefore, essentially I had to write a parser to get the required data which readOGR() does it one line.