Mapping with R and Postgresql

+ Mapping with R and Postgresql + A somewhat more extended example, will use the RPostgreSQL library and several other libaries to produce a map from the coordinates stored in the database. We will define a few helper functions to provide a layer of abstraction, which will make the R code more reusable. + #load some dependent libraries + library(maps) + library(maptools) + library(ColorBrewer) + library(ClassInt) + library(RPostgreSQL) + +#Define some helper functions + +#Returns a dataframe from the connection for a valid statement +dfFromSQL<-function (con,sql){ + rs<-dbSendQuery(con,sql) + result<-fetch(rs,n=-1) + return(result) +} +#Returns a list of latitudes and + longitudes from the orgunit table +dhisGetFacilityCoordinates<- function(con,levelLimit=4) { +sqlCoords<-paste("SELECT ou.organisationunitid, ou.name, +substring(ou.coordinates from E'(?=,?)-[0-9]+\\.[0-9]+')::double precision as latitude, +substring(ou.coordinates from E'[0-9\\.]+')::double precision as + longitude FROM organisationunit ou where ou.organisationunitid + in (SELECT DISTINCT idlevel",levelLimit, " from _orgunitstructure) + and ou.featuretype = 'Point' + ;",sep="") + result<-dfFromSQL(con,sqlCoords) + return(result) + } + +#Gets a dataframe of IndicatorValues, +# provided the name of the indicator, +# startdate, periodtype and level +dhisGetAggregatedIndicatorValues<-function(con, +indicatorName, +startdate, +periodtype="Yearly", +level=4) +{ + sql<-paste("SELECT organisationunitid,dv.value FROM aggregatedindicatorvalue dv +where dv.indicatorid = +(SELECT indicatorid from indicator where name = \'",indicatorName,"\') and dv.level + =", level,"and + dv.periodid = +(SELECT periodid from period where +startdate = \'",startdate,"\' +and periodtypeid = +(SELECT periodtypeid from periodtype + where name = \'",periodtype,"\'));",sep="") + result<-dfFromSQL(con,sql) + return(result) + } + +#Main function which handles the plotting. +#con is the database connection +#IndicatorName is the name of the Indicator +#StartDate is the startdate +#baselayer is the baselayer +plotIndicator<-function(con, +IndicatorName, +StartDate, +periodtype="Yearly", +level=4,baselayer) +{ +myDF<-dhisGetAggregatedIndicatorValues(con, +IndicatorName,StartDate,periodtype,level) +coords<-dhisGetFacilityCoordinates(con,level) +myDF<-merge(myDF,coords) +myDF<-SpatialPointsDataFrame(myDF[, +c("longitude","latitude")],myDF) +IndColors<-c("firebrick4","firebrick1","gold" +,"darkolivegreen1","darkgreen") +class<-classIntervals(myDF$value,n=6,style="quantile" +,pal=IndColors) +colCode<-findColours(class,IndColors) +myPlot<-plot.new() +plot(baselayer) +points(myDF,col=colCode,pch=19) +title(main=IndicatorName,sub=StartDate) +return(myPlot) } + + + + Up until this point, we have defined a few functions to help us make a map. We need to get the coordinates stored in the database and merge these with the indicator which we plan to map. We then retrieve the data from the aggregated indicator table, create a special type of data frame (SpatialPointsDataFrame), apply some styling to this, and then create the plot. + +#Now we define the actual thing to do +con <- dbConnect(PostgreSQL(), user= "dhis", password="SomethingSecure", dbname="dhis") +#Define the name of the indicator to plot +MyIndicatorName<-"Total OPD Attendance" +MyPeriodType<-"Yearly" +#This should match the level where coordinates are stored +MyLevel<-4 +#Given the startdate and period type, it is enough +#to determine the period +MyStartDate<-"2010-01-01" +#Get some Some Zambia district data from GADM +con <- url("http://www.filefactory.com/file/c2a3898/n/ZMB_adm2_RData") +print(load(con))#saved as gadm object +#Make the map +plotIndicator(con,MyIndicatorName,MyStartDate,MyPeriodType,MyLevel,gadm) + The results of the plotIndicator function are shown below. + + + + + + In this example, we showed how to use the RPostgreSQL library and other helper libraries(Maptools, ColorBrewer) to create a simple map from the DHIS2 data mart. +