Deforestation in Southern Mexico: Displaying PostGIS raster queries using Google maps

The recent publication of Hansen and colleagues’ global deforestation map in a user friendly Google maps earth explorer format (including street view for Mexico) has received a great deal of publicity.
So I placed the Landsat classification for Southern Mexico that our group published in Google maps using more or less the same format.
deformeasure1

Click here to explore the map.

We are using a desktop PC as a provisional server from Ecosur. The link will not be 100% reliable as we do have fairly frequent power failures. However the map layers use GeoWebCache which works well to speed up the visualisation even with the very slow internet connections we have

I have been struggling to get my head around the logic of the Google maps API. The map needs the addition of a transparency slider and a legend to be really usable, which I have not managed to add yet. The colours follow Hansen’s scheme with red for deforestation between 1990 and 2000 and orange for deforestation between 2000 and 2005. This should be shown on the map in some way.

The fun bit was implementing some code to make simple calculations around a buffer using PostGIS and R. This part was very easy within PostGIS. However as I have no experience with the Google maps API, nor knowledge of Java, I had more difficulty showing the results in the Google interface. For some reason whenever I added a click listener I lost the toggle control (and vice versa). I adapted the Google map code from examples so I wasn’t quite sure what I was doing with this part. So for the moment there are two versions of the map that need unifying.

Click here to explore the second map. Drop markers and click on them to get the calculations.

The second map produces results like the screenshot below. The Goggle maps source code for both can be read in the browser. So if anyone with experience with Google maps could help me out by adding a transparency slider and unifying the two maps into one I would be very thankful!

deformeasure

The second map sends the result of a click to a php file. This runs a query in PostGIS that includes some PLR code to produce the figure. Scripting this bit was the easy part. The approach could be adapted to do many useful things. I used php to link the two as it is very simple to understand.

<?php

include('/usr/share/secret/dbaccess.php');
// Connecting, selecting database
$dbconn = pg_connect("host=localhost dbname=reddeam user=$user password=$pass")
    or die('Could not connect: ' . pg_last_error());

// Performing SQL query

$Lat = $_GET["Lat"];
$Long = $_GET["Long"];



echo "<h1 Deforestation>";
$query = "select (ret).id, (ret).prop*100 percent from
(select forcov((st_dumpvalues(st_union(st_clip(rast,1,geom)))).valarray) ret from 
(select st_buffer(st_setsrid(st_makepoint('$Long','$Lat'),4326)::geography,10000)::geometry geom) o, 
defor  where st_intersects(rast,geom)) s";
$result = pg_query($query) or die('Query failed: ' . pg_last_error());
// Printing results in HTML

$src = '/var/www/apps/figs/temp/pie.png';
$new = "/var/www/apps/figs/temp/pie$Lat.png";
rename($src , $new);

echo "<h2>Cubertura forestal a un radio de 10 km</h2>";

echo "<table>\n";
while ($line = pg_fetch_array($result, null, PGSQL_ASSOC)) {
    echo "\t<tr>\n";
    foreach ($line as $col_value) {
        echo "\t\t<td>$col_value</td>\n";
    }
    echo "\t</tr>\n";
}
echo "</table>\n";
echo "<center>\n";
echo "<img src='figs/temp/pie$Lat.png'?forceref==$Lat height='100' width='100' alt=''>";
echo "</center>";

// Free resultset
pg_free_result($result);

// Closing connection
pg_close($dbconn);
?>

The only real difficulty I had was forcing the browser to refresh the little pie chart. As the file had the same name each time the browser used a cached version. I eventually added these lines to change the file name.

$src = ‘/var/www/apps/figs/temp/pie.png’;
$new = “/var/www/apps/figs/temp/pie$Lat.png”;
rename($src , $new);

Then used a system command from R to clean out the figures before making the next one.

The PLR function is below.

create type ret2 as (id text, prop float);

CREATE OR REPLACE FUNCTION forcov (float[]) RETURNS setof ret2 AS '
rm(x)
system("rm /var/www/apps/figs/temp/pie*.png") # get rid of old images
x<-as.vector(unlist(arg1)) # Not sure how much of all this cleaning up is needed!
x<-gsub("\\{","",x)
x<-gsub("\\}","",x)
x<-gsub("NULL","NA",x)
x<-as.numeric(as.character(x))
x<-na.omit(x)
xx<-x
#xx<-x[x%in%c(111,222,122,112,211,221)]
## This bit merges the results into a dataframe with the right colours for the pie.
x[x==111]<-"Bosque 1990-2006"
x[x==222]<-"No bosque antes de 1990"
x[x==122]<-"Deforestado entre 1990 y 2000"
x[x==112]<-"Deforestado entre 2000 y 2006"
x[x==211]<-"Regen entre 1990 y 2006"
x[x==221]<-"Regen entre 1990 y 2006"
x[x==212]<-NA
x[x==444]<-NA
x<-na.omit(x)
d<-data.frame(xx=c(111,222,122,112,211,221,212),colors=c("green","black","red","orange","blue","blue","pink"))
png(width=400,height=400,file="/var/www/apps/figs/temp/pie.png")
dd<-na.omit(merge(d,data.frame(table(xx))))
pie(dd$Freq,col=as.character(dd$colors))
dev.off()
## File permissions need to be changed.
system("chmod 635 /var/www/apps/figs/temp/pie.png")

data.frame(round(prop.table((table(x))),4))

' 
LANGUAGE 'plr' STRICT;

Reference

http://www.plosone.org/article/info:doi/10.1371/journal.pone.0042309

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s