Demonstrating R and PostgreSQL Inter-Operability: Rice farm Site Selection

INTRODUCTION

In this project we are to find the best locations for rice farming in Ghana using a series of SQL statements.We would demonstrate by operating PostgreSQL in R.

In this projected we are provided with a districts table (spatial data) in a form of table and also a rivers table (spatial data)

Originally this data was a district.shp and a rivers.shp, but we should have it at the back of our minds that,as soon as we load such data into Postgres they now regarded as tables.

NB:All data used in this work-flow are in refrence to the yaer 2010

Objectives

1. To familiarize oneself with R inter-operability with Spatial Databases(eg.PostgreSQL) and to be able to set up the necessary environment.
2. Learn to use some basic Structured Query language(SQL) and Spatial Query language(‘SQL’).

NB:The operations we would be performing can equally be done using the sf package in R with some assistance with the dplyr package. It can also be performed in GIS applications like QGIS ,ArcGIS etc.

Loading of libraries

library(DBI)
library(odbc)
library(sf)
library(dplyr)
library(htmltools)
library(mapview)

Setting up connection

Lets first set up the connection to allow the inter-operability between R and PostgreSQL.

#Already set up the named DSN in windows.
mycon<- DBI::dbConnect(odbc::odbc(),dsn='benDSN')

You can equally set up your connection also in R. In terms of security ,I would advise you setup the connection in windows(preferable).

mycon <- DBI::dbConnect(odbc::odbc(),
                      Driver ="{PostgreSQL ANSI}",
                      Server ='localhost',
                      Port =5432,
                      Database ='Lesson4db',
                      UID ='postgres',
                      PWD ='password' 
                      ) 

Let’s view the first 6 features of the districts table.

SELECT *
FROM ghana.districts
LIMIT 6;

Let’s view the fisrt 6 features of the rivers table

SELECT *
FROM ghana.rivers
LIMIT 6;

Let’s examine the SRID of the districts table column geom(this hold the geometry makeup).

SELECT 
    Find_SRID('ghana','districts','geom');
    

SRID 25000 is a Transverse Mercartor projection with units in metres.

We also examine the SRID of the rivers table.

SELECT 
  Find_SRID('ghana','rivers','geom');

Suitable Location for Rice farming

1. Population of district capital less than 10000m.
2. Within at least 10km (10000m) of river.
3. Should be within the Northern Region or the Volta Region.

SELECT d.region,d.district, d.capital,d.geom
FROM ghana.districts d
INNER JOIN ghana.rivers r
ON ST_overlaps(d.geom, ST_buffer(r.geom, 10000))
WHERE d.popcap < 10000
AND (d.region = 'Northern' OR d.region = 'Volta')
GROUP BY (d.capital, d.district, d.region, d.geom)
ORDER BY (d.region) ;

Create View

Let’s create a virtual table(view), this allows us to easily pipe the results of our analysis into any GIS application. A typical example is Qgis. We would need the column gid as primary Key when piping our view into Qgis.

NB:Th view created ricefarm_locations would be stored in the schema= ghana. Also, tables and views are written directly into our database in PostgreSQL.

CREATE VIEW  ghana.ricefarm_bestlocation AS(
SELECT d.gid, d.region,d.district, d.capital, d.geom
FROM ghana.districts d
INNER JOIN ghana.rivers r
ON ST_overlaps(d.geom, ST_buffer(r.geom, 10000))
WHERE d.popcap < 10000
AND (d.region = 'Northern' OR d.region = 'Volta')
GROUP BY (d.gid,d.capital, d.district, d.region, d.geom)
ORDER BY (d.region)
);

Create Table

We would love to plot the ouput of our geo-processing analysis in R. Hence the is a need to save the output as a table and an output variable in R.

Note: This would output a dataframe.

SELECT region,district,capital, ST_ASEWKT(geom) AS geom
FROM ghana.ricefarm_bestlocation ;

Let’s convert ricefarm_bestlocation from a dataframe to an sf object using st_as_sf(). We would just re-assign to same variable.

ricefarm_bestlocation <- st_as_sf(ricefarm_bestlocation,
                              wkt = 'geom')

Let’s ascertain the class and crs(epsg) of the ricefarm_selected object.

cat("Class ricefarm_bestlocation \n"); class(ricefarm_bestlocation) ; 
Class ricefarm_bestlocation 
[1] "sf"         "data.frame"

cat("CRS ricefarm_bestlocation \n"); st_crs(ricefarm_bestlocation)$epsg 
CRS ricefarm_bestlocation 
[1] 25000

Interactive Map with Mapview

Let’s create a basic interactive map of selected suitable location ricefarm_bestlocation in mapview.

addLabel <- function(data) {                          
  data$label <- paste0(
    '<table style="width:120px;">
      <tr>  
        <td><b>Region</b>: ',ricefarm_bestlocation$`region`,'</td>
      </tr>
      <tr>  
        <td><b>District</b>: ',ricefarm_bestlocation$`district`,'</td>
      </tr>
    </table>'
  )
  data$label <- lapply(data$label, HTML)

  return(data$label)
}

# We now create the interactive map
mapView(ricefarm_bestlocation,
        label= addLabel(ricefarm_bestlocation),
        layer.name= "Suitable Locations",
        map.types="OpenStreetMap",
        col.regions='#31a354')

DBI::dbDisconnect(mycon)