Duncan Golicher’s weblog

PostGIS does not support raster data as yet. Hopefully there should be raster support soon. Raster layers can be visualised from a local or remote server using geoserver. However having raster layers integrated into a PostGIS data base would be extremely useful for tasks such as overlaying herbarium records with continuous extrapolated climate data or elevation.

I have previously experimented with importing the points of a raster layer. This works quite well. It can be used for simple spatial overlays using proximity operations. However the visual impression produced by the points is not very satisfactory and a query using the contains operator would be more efficient.

QGIS has a plugin for graticule generation. This forms a grid of individual polygons. So if a grid is produced using the graticule plugin of QGIS and then imported into PostGIS using SPIT it is possible to partially mimic a raster layer, although with some computational overhead.

The document mnp-workshop3.pdf showed how worldclim data (http://www.worldclim.org/) can be imported into GRASS and R.

If R is run from GRASS the following lines can be used to import the appropriate raster layers at a chosen resolution into R.

Most of the lines of the script are concerned with calculations on the basic worldclim data. The dataset consists of 36 grids of monthly mean minimum temperature, maximum temperature and precipitation. Some mean annual values may be derived for simplicity.

a<-c(paste(”tmin”,1:12,sep=”"),paste(”tmax”,1:12,sep=”"),paste(”prec”,1:12,sep=”") ,”alt”)
d<-readRAST6(a)
d@data[,1:24]<-lapply(d@data[,1:24],function(x)x/10)
fullgrid(d)<-FALSE

f<-function(x)mean(x[1:12])
d@data[["meanmim"]]<-apply(d@data,1,f)

f<-function(x)mean(x[13:24])
d@data[["meanmax"]]<-apply(d@data,1,f)

f<-function(x)mean(x[1:24])
d@data[["meantemp"]]<-apply(d@data,1,f)

f<-function(x)sum(x[25:36])
d@data[["totprec"]]<-apply(d@data,1,f)

The trick for importing the layers to PostGIS from R is shown below. If the data is in R as a spatialPixelDataFrame the coordinates can be added to the data to form a table.

dd<-data.frame(coordinates(d),d@data)

This table can then be saved to a PostGIS database using ODBC.

library(RODBC)
con<-odbcConnect(”mydb”)
odbcQuery(con,”SET search_path =chiapas, pg_catalog;”)

sqlSave(con,dd,”worldclim”)

###########################################

Now either in pgadmin or running the queries directly through R by placing them inside an odbcQuery function.

ALTER TABLE chiapas.worldclim add column gid serial NOT NULL;
ALTER TABLE chiapas.worldclim ADD CONSTRAINT worldclim_pkey PRIMARY KEY(gid);
ALTER TABLE chiapas.worldclim ADD COLUMN the_geom geometry;

UPDATE chiapas.worldclim
SET the_geom = PointFromText(’POINT(’ || s1 || ‘ ‘ || s2 || ‘)’,4326);

CREATE INDEX
grat30arc_idx_geo on grat30arc USING GIST(the_geom GIST_GEOMETRY_OPS);

CREATE INDEX
worldclim_idx_geo on chiapas.worldclim USING GIST(the_geom GIST_GEOMETRY_OPS);

If a polygon graticule called grat30arc covering the same region with a 30 arch second grid has been produced using the graticule plugin in Qgis and imported into the public schema then the average values of the points within each grid square can then be calculated for the graticule and a new table produced using this as a model query.

CREATE TABLE chiapas.worldclim2 with oids as
SELECT a.*, g.totprec, g.meanmin, g.meanmax, g.elev from
grat30arc a,

(SELECT AVG(c.totprec) as totprec,
AVG(c.meanmim) as meanmin,
AVG(c.meanmax) as meanmax,
AVG(c.alt) as elev,d.gid
FROM chiapas.worldclim c, grat30arc d
WHERE c.the_geom && d.the_geom
AND contains(d.the_geom, c.the_geom) group by d.gid) g

WHERE a.gid=g.gid;

The result is very similar to the visualisation of an arcgrid in arcview when it is viewed in QGis over a shaded relief map. In this case the pixels are approximately 1 km x 1 km (30 arc seconds or 0.008333333 degrees) and the elevation lines up well with the contours derived from a 1:250,000 map.

Spatial indices on the herbarium table and the worldclim2 layer with the following

CREATE INDEX
herb_idx_geo on herb.herbario USING GIST(the_geom GIST_GEOMETRY_OPS);

CREATE INDEX
worldclim2_idx_geo on chiapas.worldclim2 USING GIST(the_geom GIST_GEOMETRY_OPS);
Then the following query runs in less than four seconds to provide 6,000 herbarium records with extrapolated climate data to the nearest km.

SELECT h.*, w.totprec, w.elev,w.meanmin,w.meanmax
FROM herb.herbario h, chiapas.worldclim2 w
WHERE h.the_geom && w.the_geom
AND contains(w.the_geom, h.the_geom);

Or

select h.*, w.totprec, w.elev,w.meanmin,w.meanmax
from herb.herbario h, chiapas.worldclim2 w
where ST_contains(w.the_geom, h.the_geom);

No creo que sea apto usar este sitio para hacer comentarios con respecto a características de la institución donde trabajo. Nada mas aprovecho una oportunidad para citar un icono cultural.

“Turning and turning in the widening gyre
The falcon cannot hear the falconer;
Things fall apart; the centre cannot hold;
Mere anarchy is loosed upon the world,
The blood-dimmed tide is loosed, and everywhere
The ceremony of innocence is drowned;
The best lack all conviction, while the worst
Are full of passionate intensity.
Surely some revelation is at hand;
Surely the Second Coming is at hand.”
- W. B. Yeats (‘The Second Coming’)

Mapservers have never taken centre stage in my thinking on how to make Ecosur’s geographical data more widely available.  When using a geographical data base the underlying information is almost by definition available for analysis. In contrast data served up by a mapserver is usually only suitable for visualisation. However mapservers have an important role in conjunction with PostGIS. The NASA JPL wms server is a an essential provider of satelite images that provide the backdrop that allows a visual intepretation of spatial information.

At a recent workshop Conabio provided me with a very nice shaded relief map for Mexico derived from a digital elevation model. I have been using it locally with Qgis. However I became concerned that it must also be potentially made available for use with POSTGIS at an institutional level. This clearly needs a mapserver.

I have previously wasted a fair ammount of time failing to configure mapservers correctly. So this morning I looked with some trepidation at geoserver. To my relief within an hour I had geoserver configured and the shaded image available to QGIS through a local host WMS.

The steps are extremely simple:

1. Download geoserver as a zip file

2. Extract the contents.

3. Set an environmental variable to point to your implementation of Java. I have sun java 6 installed using the medibuntu repository.

sudo gedit /etc/environment

Add this line

JAVA_HOME=”/usr/lib/jvm/java-6-sun”

4. Reboot and run mapserver using a command that goes something like

/home/duncan/geoserver/bin/startup.sh

5. Point the browser at

http://localhost:8080/geoserver/

Then if everything is working simply follow all the online instructions in order to add layers.

With geoserver started there is a new wms server that QGis can connect to.

It is also easy to include data from geoserver in Google Earth. Choose <add> <image overlay> from the menu. Then the “refresh tab”. Click on WMS parameters and add the geoserver details.

You then have the layer in Google Earth

Ecosur will soon be providing students and researchers with a networked data base for storing and visualising public domain spatial layers. This database should also store results from ongoing research. I am convinced that the advantages of providing a fully functional geographical data base at zero cost by using PostGIS/Postgresql are compelling. The release of a stable Postgresql 8.3 and Qgis 0.10 have made the task even simpler.

Postgresql 8.3 is now notably faster when running spatial queries than the version I had been using (8.1). It also has very advanced features.  These move PostgreSQL up to the power of the premium  proprietary enterprise data base, Oracle. They clearly go beyond those found in MySQL.

The only practical challenge I faced with PostgreSQL 8.1 involved a very simple issue of scalability when the database was viewed with Qgis 0.9. Despite some issues with the stability of Qgis (that did worry me, but appear to have been largely resolved) Qgis is a suitable visualisation tool for PostGIS. I have used it extensively since I started experimenting with PostGIS a few months ago and am generally satisfied. However a simple practical issue concerned me.

How could users easily be given access to the parts of a potentially very large database in such a way that they see only the part they need or are allowed to see?

The answer had clearly to be through the use of Postgresql schemas, but the Qgis interface didn’t show them in a user friendly way that matched the need.

This problem is now largely resolved by the new version of Qgis. Schemas now show up in an unfoldable tree structure. As groups of users and individual users can be given access to a subset of schemas it should now be feasible to design a database that hides all the complexity from users and allows them to quickly find the layers  they need.

Algunos estudiantes llevaron discos de Ubuntu Heron despues de mi presentación. Si han probado instalarlo estoy seguro que van a querer usarlo para ver videos. Para conseguir Flash para You Tube la mejor forma es abrir synaptics (sistema-administración-synaptics) y buscar “flashplugin-nonfree”. Luego los codecs se instalan cuando se necesitan.

Para tener Medibuntu (el repositorio para algunos codecs comerciales) la forma mas rapida es con comandos. Abre un terminal y pega este para abrir el archivo de los repositorios.
sudo gedit /etc/apt/sources.list

Luego pega este en el documento al fondo y guardarlo.

deb http://packages.medibuntu.org/ hardy free non-free

Para añadir la clave para activarlo pega este en un terminal.
wget -q http://packages.medibuntu.org/medibuntu-key.gpg -O- | sudo apt-key add - && sudo apt-get update

sudo apt-get update

Y luego nada mas añade los libs.

sudo apt-get install libdvdcss2 libdvdread3

No te asustas con el uso del terminal. Se puede conseguir lo mismo con el GUI, pero es mas facil comunicar la información como comandos.

I am impressed by Ubuntu Heron. The most important difference is  ease of installation. Both Feisty and Gutsy were simple enough to get running, given a reasonable degree of patience. However there were always a few problems to fix. These involved  screen resolution sound and video. Getting everything right needed  a search which produced multiple answers on the users forums.

A single issue that needs to be resolved by editing a config file is more than enough to get most Windows users, acustomed to  preinstalled operating systems, to give up and go strait back to what they know. Thus despite the clear merit of the underlying operating system, uptake remained much slower than it could have been among academic colleagues here.

In contrast my  experience installing Heron on a Toshiba Tecra was effortless. It all worked out of the box. What was more impressive was the short amount of time it took to get everything, including a large number of R packages, GRASS, QGis, Google Earth, MySQL, Apache, Postgresql etc all back together. I already knew what I wanted and where to get it, but even so the short time line is worth stressing.

3:15 Open Laptop and remove old hard disk (It was time to go from 80GB to 250GB and this is a great way to make the change without any fear of losing time)

3:20 New hard disk in place.

3: 25 Heron install CD is asking me for details of location, keyboard and partitioning.

3:41 Heron installed and running.

3:45 Reboot in order to use NVidia, this was recognised as necessary after activating advanced visual effects. No manual configuration needed at all.

5:03 R, GRASS, QGis, Google Earth, Flash plugin + codecs, Apache MySQL, Postgresql all installed.

So I will never worry about a hard disk crash again. Providing I have my personally generated data safe somewhere, all the rest is easy. Perhaps it is not so simple on all laptops, but it is just so encouraging to see that Ubuntu is getting within a hair of the grandma friendly Linux.

I occasionally do have need of Windows XP. So this morning I moved the .virtualbox directory over from my old home directory to the new one, after a quick install of virtual box and adding my usename to the virtualbox group. It worked exactly as it had before. The advantage of virtualising windows is so glaringly obvious I am constantly surprised by my colleagues reluctance to protect themselves against virus attack that way. It is understandable to want to continue using what you are used to (Windows XP). However it is very strange not to want to make sure that you can get back to work quickly after a disaster. Virtualising XP gives the confidence that you can be back up and running in less than an hour. As desktop PCs are more or less generic virtual machines can be backed up and moved freely between them even if an IT department hasn’t done what it should and provided a sensible networked option of virtual machines for Windows XP users.

It is worth being aware that Heron comes with Firefox 3 beta preinstalled, which is fine, as it will soon be the official version. However in the meantime you can install Firefox 2 in order to use  plugins that don’t yet work with the new version. Be aware that  you must  close down Firefox 3 in order to run firefox-2 (type firefox-2 on a command line or set up a launcher).

I also noticed that the Heron repositories has up to date versions of other programs that I didn’t realise had progressed. The most important development for me at the moment was that Postgresql8.3 is now the stable version and that QGis is 0.10.0 instead of 0.9.1. Despite still not getting to version 1.0, QGis  does look as if it has improved greatly in terms of both stability and appearance. This is encouraging after my recent presentations. QGis 0.9 did have a few uncomfortable issues that I didn’t mention. Also Open Office is better. The most notable improvement is when you try to change the language setting. Anyone who knows the frustrating illogicality of being presented with “thesaurus” when you want to check the spelling of a Spanish document will know what I mean. The options are now in much more sensible places.

Despite recent posts, I admit to being a long term Windows user myself. I therefore share many of the characteristics of Windows users. The most important of these is that I am naturally extremely conservative and very reluctant to change a laptop when I have things working the way I am used to. I have been using Ubuntu Feisty since late July 2007 without ever  bothering to upgrade to Gutsy (although I installed it for my son).

However I decided that the release of Heron really deserved my attention. I was becoming increasingly  frustrated with the limited disk space on my duel boot Toshiba Tecra with its original 80 GB HD (40 GB still in an XP partition that I now never use). So today I treated myself to a 250GB external Simpletech ATA hardrive. On getting it home I  invalidated the warranty by opening the case, taking out the drive and swopping it for my 80GB ATA drive. I mention the fact that it is ATA as older (my Toshiba was bought in October 2006) or cheaper laptops are likely to have IDE drives. Always check first and be aware that this simple trick will not work if your machine uses IDE. But if you need a new hard drive, and can afford one, simply switching the drives saves all worries about upgrading. If it goes wrong just put the old drive back in.

In went the Heron CD (I decided not to bother with any native windows install this time) and 15 minutes later all the basics were installed without any tweaking at all. I turned on the advanced visual effects. Ubuntu realised I needed an Nvidia driver and installed it. After rebooting I then had all the  wobbly windows and spinning cubes if might ever want, if I ever feel like using them. A completely painless, risk free upgrade in no time at all. Now all I need to do is copy the parts of my home director I want to keep back over from the external drive and it is all back where it was, with extras. Two hours maximum for the whole process, including putting back R, GRASS, Qgis MYSQL etc. PostGis is just a little bit more difficult in my case due to the fact that I previously used Postgreql 8.1 and it is time for 8.2. But that detail is fairly tecnical.  Contrast the basic experience with a complete re- install of Vista. It would be a nightmare that could last the best part of a  week.

Last night I read this very useful and sympathetic account of the install experience from a novices’ perspective http://www.bspcn.com/2008/04/28/the-great-ubuntu-girlfriend-experiment/

As I had handed out several disks of the Heron after the Ecosur seminar I thought I should be aware of these sort of problems.  I am quite likely to have a few students coming into my office in the next few days asking for help. The girlfriend gotchas really are all very minor. My  experience with the new Heron install was simple to the point of being disappointingly trivial. I missed all the thrills of fixing screen resolution and sound. The only important detail worth mentioning is to make sure that the flash plugin is installed from symantics before going to watch a you tube video. There are then no problems. You simply get a polite message asking for permission to install codecs and there is no need to try to install a gzipped file from adobe. Done.

1. ¿Necesito Linux para usarlo?

No. El modelo es lo que se llama “cliente-servidor”. Es decir que PostGIS es el servidor. Proporciona los datos, pero la visualización de los datos o el procesamiento de las consultas es al nivel del cliente. Tu puedes escoger el cliente que mas te gusta. En este momento hay varios programas clientes de acceso abierto (gratis) que pueden leer, procesar y visualizar información de un servidor de PostGIS. Lo mas fácil de usar es probablemente QGIS, aunque Udig también tiene buenas características. Open Jump es otra alternativa. Todos estos programas se puede instalar en Windows. Los clientes están en un proceso de desarrollo todavía y van a incrementar su utilidad con el tiempo. PostGIS ya tiene toda la funcionalidad de servir datos a clientes con mas capacidades.
Al mismo tiempo programas como Excel y Access pueden vincularse con el servidor para procesar datos no espaciales.

2. ¿Es difícil instalar Post GIS? Necesitamos apoyo técnico externo para hacerlo.

No. Se instala PostGIS con un servidor de Linux en alrededor de una media hora. Se requiere un poco de planificación y esfuerzo adicional para la configuración óptima de las opciones de seguridad. Pero es importante recordar que los usuarios de la información no van a instalar PostGIS. Ellos pueden bajar un archivo qgis.exe del internet y instalarlo bajo Windows en minutos. Si usan Linux también van a instalar Qgis, o Udig, no PostGIS. El uso de Qgis es sencillo e intuitivo. Se puede dar capacitación en su uso practico en el curso de una mañana.

3. Hemos intentado bajar capas de un mapserver, pero parece complicado.

La relación entre PostGIS y un mapserver es potencialmente un elemento potente, pero no se necesita correr un mapserver para usar PostGIS. Un mapserver recibe capas de PostGIS, Asi que es otra capa entre el usuario y servidor. La configuración de un mapserver para correr con PostGIS no es muy difícil. Hay dos mapservers de acceso abierto (mapserver y geoserver) que pueden usarse. Sin embargo un mapserver es realmente una forma de dar información al publico. El primer paso hacia el uso de un mapserver sería asegurar la velocidad de la conectividad al Internet de Ecosur.

4.Estoy acostumbrado de usar ArcView, tengo licencia para usarlo y muchos shapefiles disponibles. ¿Por qué me interesa esta iniciativa?

La ventaja de una base de datos espacial es que toda la información esta guardada en el mismo lugar, bajo el mismo sistema de coordenadas y actualizado. Las capas pueden comunicarse entre ellas. Si ves por ejemplo que las coordenadas del censo nacional no coinciden con los del conteo de población sería posible tratar de resolver el problema. Una vez resuelto todos los usuarios de la base no van a repetir los errores. Usuarios de ArcView pueden bajar la información como shapefiles con un clic y trabajarla al nivel local de sus propios Pcs. Es muy importante enfatizar que la iniciativa no fuerza usuarios de cambiar su forma de trabajar. Cada persona desarrolla sus propias preferencias por software al nivel de sus propia “desktop” y hay que respetar sus decisiones. Pero al nivel del servidor hay que asegurar que todos los datos están mantenidos bajo estándares rigurosos, uniformes y convencionales.

5.¿Cuando podemos contar con un servicio así?

Un prototipo ya esta funcional y ya puede ser muy útil. Se puede proporcionar una herramienta de uso interno en el espacio de una semana. Pero la verdadera potencia de un sistema así viene con mas esfuerzo. Se puede programar portales Web atractivos para el publico. Se puede organizar sistemas amigables para que estudiantes entren sus datos en el sistema sin tener que saber detalles sobre el uso de SQL. El control de calidad y mejoramiento de la información se consigue a través de esfuerzos de revisión meticulosos y programación con datos. Para contar con un sistema del cual podemos estar orgullosos como institución debemos emplear por lo menos una persona con capacidad en programación de bases de datos tiempo completo.
Sería vital la coordinación entre investigadores que producen nuevos datos de campo, las colecciones, LAIGE, informática y la biblioteca.

6. ¿Cuales son la ventajas de Linux para mi como individuo?

Mis propias razones estan expresados aqui.

http://duncanjg.wordpress.com/2008/02/08/por-que-linux/

Abajo he incluido las diapositivas de una presentación que voy a dar hoy en Ecosur sobre la potencia de PostGIS como un base de datos para la investigación institucional en Ecosur. Incluyo una introducción sobre la filosofía del movimiento de aceso libre y la potencia del nuevo modelo “profesionalizado” del software de aceso libre.

bases-de-datos-de-aceso-libre

Los ejemplos del uso de POSTGIS se ven en los videos de You Tube abajo.

Desafortunadamente la calidad del imagen no es muy alto. De todos modos dado la velocidad de la conexión institucional no creo que se pueden ver dentro del institución.

Until comparatively recently visualisation of geographical information tended to mean one thing. ArcView.

When Google Earth (and the perhaps under rated and certainly under used NASA Worldwind) appeared on the scene visualisation of free geographical information became popular and ubiquitous. Three dimensional spinning globes revolutionized the way none GIS professionals could work with geographical data. At around the same time some comparatively modest open source projects turned into mature easily used products that could be applied to two dimensional desktop map composition.

At the time of writing the best known of these are QGis (Quantum GIS) and Udig. Both are cross platform. Udig uses Java and runs well under Windows. Qgis still performs rather better under Linux, at least in my experience.

Neither of these programs are true geographical information systems, in the sense that they do not of themselves provide very much in the way of analytical tools. They simply allow geographical information to be visualised in a convenient way. Visualisation is a fundamental element of geographical analysis so it should be done well. These open source programs are acceptable alternatives to the traditional ESRI benchmark.

Users of Arcview commonly activate a range of plugins to give the program more features. Qgis provides a connection with GRASS to provide analytical capability, although I would assume that few GRASS users routinely run GRASS through Qgis. Both Qgis and Udig have wms (world map server) connectivity, which means that they can be used to customize the sort of maps that are found on web sites showing wms data.

An interesting recent development in GIS has been the rapid growth in the use of spatially aware data bases. Most major databases (MySQL, Oracle, Postgresql) can now be given a spatial component. However in the open source context “spatial database” has become almost synonymous with the postgis component of postgresql. It is the built in connectivity with posttgis feature that gives qgis and udig their rapidly expanding potential value. Spatial databases allow sophisticated extraction and combination of geographical data. The results usually usually have to be shown as maps. Command line sql statements are not visual of themselves. QGis and Udig can thus potentially provide the full functionality of an Arc* type environment through levering the capabilities of postgis.

At present this very powerful framework is still in a rather incipient phase, in the sense that there are no GUIs that carry out analyses comparable with Arcview/ArcGis. Most of the functionality of the expensive commercial software is now available through the open source alternatives. However while the documentation is comprehensive it does usually assume prior knowledge of data base programming. Most of the current postgis users are knowledgeable data base programmers with considerable experience using spatial sql queries in the context of online mapserver based web sites. However this is changing fast. The WKB (Well Known Binary) representation of geographical information used in postgis databases looks like a good standard and will undoubtedly be more widely used. This should result in the very rapid development of higher level forms of manipulating information in this format.

One advantage of using an online data base as opposed to separate shapefiles is that the coordinate system of all data has to be consistent, or made so through reprojection when the data base is set up.This is especially important in Mexico where data is often provided without details of the projection and datum used (NAD27 and WGS84 are both used and projection details are rarely provided). Also details of how to connect to the data base using Qgis or Udig can be provided in an email, while it can be difficult to send large shapefiles. Both Qgis and Udig allow POSTGIS data to be saved locally as shapefiles.

R is of course a statistical environment, so it really shouldn’t be expected to be doing the job of GRASS (Nviz), Google Earth, WorldWind, Qgis or Udig as well. However some nice visualizations can be done with rgl, including rgb composites from satellite images. They can even be more useful for quick regional overviews than the alternatives mentioned above.

The code to produce these images be found here bluemarble3d.doc. The first line will download the imported spatial object from this site. R does need a lot of memory for images of this size, so I guess you need a minimum of 1GB RAM (I am using a Toshiba Tecra laptop with Nvidia running Ubuntu Feisty 7.04 with 2GB). Given that, you should get very own zoomable 3d image of Mexico and Central America to play with in the same time it takes the screenshot below to run. It shows the steps in real time. (For R beginners, just open the file, paste all the code into the R console and wait a little for the download). The resolution is quite coarse (2 minute) so this is only suitable for a regional overview. You will need the rgl and sp packages installed first of course.

bluemarble3d.doc

Exactly he same can be done with Landsat imagery. landsat.doc

If you want to get your own coverages for this for your own region the starting point is to run R from GRASS in a lat-long location. For example

g.region -p

projection: 3 (Latitude-Longitude)
datum: wgs84
ellipsoid: a=6378137 es=0.006694379990141317
north: 32:51:47.86586N
south: 5:28:57.975592N
west: 118:30:20.609692W
east: 76:26:43.294846W
nsres: 0:02:00.060768
ewres: 0:01:59.982024
rows: 821
cols: 1262
cells: 1036102

Then get some of the Blue Marble (Modis) imagery that can be downloaded from the NASA wms strait into GRASS using r.in.onearth. From within R you can use system to run GRASS. This does work with GRASS under Cygwin in Windows with the shell command.

system(”r.in.onearth -b file=/tmp/test month=Apr time=2005-3-24 ‘wgetopt=-c -t 5 –user-agent=MSIE5.5′ “)

The layers can then be imported into R with spgrass6.

library(spgrass6)

d<-readRAST6(c(”BMNG_Apr.red”,”BMNG_Apr.blue”,”BMNG_Apr.green”,”alt”))

Notice that this line imported a digital elevation model of my own as well, as one is of course needed for the 3d terrain effect.

The nobel prize winning economist Joseph Stiglitz is very fond of quoting Thomas Jefferson. “Knowledge is like a candle; as it lights another candle, its original light is not diminished”.

Economists like Stiglitz refer to the opposite of this property as exclusivity. In other words if I eat a box of chocolates or turn a tree into timber, you can’t have the chocolates or enjoy the shade and scenic beauty provided by the tree. Only one of us can consume a consumable. On the other hand if I have an idea I can share it with you, if you are interested. If my idea happens to be a good one then denying access to it, or the right to use it, gives rise to an inefficiency. Traditional intellectual property is based on such exclusion. Knowledge is a global public good that is of potential benefit to anyone in the world. There is a global social cost in depriving anyone in the world the right to use available knowledge.

Most of our individual ideas are “memes”. In other words we have all been lit from someone else’s candle. Originality in science is extremely difficult and is only partially achieved even by the small handful of creative individuals that manage to provide a new spark. We all worry that we are wrong or being too derivative when we think or write about scientific issues but I have (quite recently) come to the conclusion that we really shouldn’t be worried that our ideas are not original providing we are casting some light, somewhere. The more candles that are out there, the more light is cast and the collective effect is to everyone’s benefit.

In the open source model of software development the light cast has an even more positive effect, because software is itself useful as a tool for achieving other goals. Software developers in fact make some of the candles. I have neither the time nor the ability to contribute code to open source projects, but I appreciate the tremendous generosity of all those who do and try to support and use open source alternatives to commercial software as a matter of principle.

Mickey wants to work with dangerous animals, particularly snakes, when he grows up. He has started already.

If you haven’t realised, Gimp is the open source equivalent of photoshop.

Compiling data for species distribution modelling has led me to look in some depth at the IPCC SRES scenarios. Developing plausible storylines for the future development of the Global economy over the next hundred years is an extremely challenging task, However the IPCC have drawn on their considerable expertise to develop a sophisticated, consensual and inclusive set of scenarios. Without going into great detail, a very simple message struck me as important. The more optimistic scenarios always involve greater convergence between regions and a change towards a service and information economy. In other words “Globalization”.

This suggests that the internet itself could be a major factor in combating global warming. However the section of the World’s population that will suffer most of the negative consequences of global warming, the rural poor, remain largely excluded from access to the internet. I am very interested in the way developments in Open Source software can help to redress this. One of the most positive examples is the Edubuntu project that allows high quality thin client - fat server networks to be built from recycled PCs. I will continue this theme in later posts.

These are the four groups of “storylines” used in the scenario building. Note that the worst case scenario in terms of emissions and consequent climate change arises from A2 or B2.
The A1 storyline and scenario family describes a future world of very rapid economic growth, low population growth, and the rapid introduction of new and more efficient technologies. Major underlying themes are convergence among regions, capacity building, and increased cultural and social interactions, with a substantial reduction in regional differences in per capita income. The A1 scenario family develops into four groups that describe alternative directions of technological change in the energy system. Two of the fossil-intensive groups were merged in the SPM.

The A2 storyline and scenario family describes a very heterogeneous world. The underlying theme is self-reliance and preservation of local identities. Fertility patterns across regions converge very slowly, which results in high population growth. Economic development is primarily regionally oriented and per capita economic growth and technological change are more fragmented and slower than in other storylines.

The B1 storyline and scenario family describes a convergent world with the same low population growth as in the A1 storyline, but with rapid changes in economic structures toward a service and information economy, with reductions in material intensity, and the introduction of clean and resource-efficient technologies. The emphasis is on global solutions to economic, social, and environmental sustainability, including improved equity, but without additional climate initiatives

The B2 storyline and scenario family describes a world in which the emphasis is on local solutions to economic, social, and environmental sustainability. It is a world with moderate population growth, intermediate levels of economic development, and less rapid and more diverse technological change than in the B1 and A1 storylines. While the scenario is also oriented toward environmental protection and social equity, it focuses on local and regional levels.

After setting up cdo and netcdf you need to download some scenarios from cera.

The distribution modelling project we are working on is looking at a wide range of simulations in order to look at the effects of credible patterns of change in the global climate on tree species diversity in MesoAmerica. There are still relatively few simulations that have specifically represented the new emisions scenarios in the IPCC fourth assessment report. Many of the available model data sets were derived using the Special Report on Emission Scenarios (IPCC 2000, SRES).

The SRES data sets were classified into four different scenario families (A1, A2, B1, B2).
SRES. The A2 storyline describes a very heterogeneous world with the underlying theme of self-reliance and preservation of local identities. It results in this scenario a continous increasing population together with a slower economic growth and technological change. The Hadley Centre Circulation Model is a 3-dim AOGCM described by (Gordon et al., 2000 and Pope et al., 2000).
The atmospheric component has a 19 levels horizontal resolution, comparable with spectral resolution of T42, while the ocean component has a 20 levels resolution. HADCM3(http://www.meto.gov.uk/research/hadleycentre/models/HadCM3.html )

Here I will show a simple use way to look at the HADCM3_A2 scenario in R.

For simple interpretation and use of climate model output three variables are of most interest. Maximum temperature, minimum temperature and precipitation. A great deal of processing can be done in cdo,but here we move the precipitation data set to R using cdo to rewrite the file in netcdf format.

cdo -f nc copy HADCM3_A2_prec_1-1800.grb HADCM3_A2_prec_1-1800.nc

Now in R we can open the file using ncdf.

library(ncdf)
library(maptools)
library(maps)
library(RColorBrewer)

prec.nc <- open.ncdf(”HADCM3_A2_prec_1-1800.nc”)

x <-get.var.ncdf( prec.nc, “lon”)
y <-get.var.ncdf( prec.nc, “lat”)
z<- get.var.ncdf( prec.nc, “var61″)
time <-get.var.ncdf( prec.nc, “time”)

This reads all the data in. Notice that this global grid is quite a low resolution so the whole simulation can be read into memory in one go. Subsets of netcdf files can be taken if they are larger.

The order of the x and y vector has to be changed to be strictly ascending and the matrix rearranged accordingly.

z<-z[,73:1,]
y<-sort(y)
x[x>180]<- (-(360-x[x>180]))
xorder<-order(x)
z<-z[xorder,,]
x<-sort(x)

Now we can average the first fifty years monthly precipitation and subtract this from each month’s simulated precipitation in order to show the simulated pattern of rainfall anomalies over the next fifty years.

for (i in 600:1200){
j<-i%%12;j[j<1]<-12
a<-seq(j,600,by=12)
aa<-apply(z[,,a],c(1,2),mean)
a<-z[,,i] -aa
image(x,y,a,col=brewer.pal(8,”RdBu”),zlim=c(-10,10))
map(”world”,add=T)
title(main=paste(substr(time[i],1,4),substr(time[i],5,6)))
Sys.sleep(0.2)
}

It is not easy to draw conclusions from this visualization alone. However what seems clear is that the largest anomalies and fluctuations in rainfall patterns (areas of dark red or dark blue) are concentrated in the tropics. This is in part due to the fact that the shading represents absolute variation. The tropics receive more rainfall, thus variability is naturally higher. What would be useful would be look in more detail of the pattern for different areas of the global simulation and contrast them after taking into account natural interannual variability.

In order to look at the trend for any part of the world we can make a small function in R that will select a pixel from the map using the locator function and then use stl to produce a seasonally adjusted analysis similar to the one used here.

plot.trend<-function(){
image(x,y,z[,,1],col=brewer.pal(8,”RdBu”)[8:1],zlim=c(-10,10))
map(”world”,add=T)
l<-locator(1)
a<-z[which(abs(l$x-x)==min(abs(l$x-x))),which(abs(l$y-y)==min(abs(l$y-y))),]
a<-ts(a,start=c(1950,1),frequency=12)
plot(stl(a,s.window=12))
}

Select Chiapas.

The units are in mm per day. Thus the third graph represents the annual trend after seasonal effects (second panel) are extracted. At present averaged over the Chiapas grid square annual rainfall averages just over 3 mm per day ( 1000 - 1200 mm per year) . There is around 40% variability in this between years. The simulation suggests that in general terms the pattern of variability will be maintained for the next fifty years. However the end of the simulation shows a rather sudden, sharp drop in rainfall by the end of the century.

There is a great deal of uncertainty in climate model simulations, not only regarding the amount of greenhouse gas forcing but also the complex patterns of feedbacks in the global systems. An alternative scenario from HADCM3 (the B2 scenario) shows a slightly different pattern for Southern Mexico.

Installing CDO and netcdf

Global circulation models that iterate calculations on a global grid over many time steps produce potentially vast amounts of data. Output from climate models and large historical climatic data are often held in either GRIB (WMO) or netCDF (Unidata) format. These formats have been specifically designed for efficient, compact storage of huge volumes of information. Research groups that work on climate modelling and monitoring have developed their own specialized tools for handling data in these formats under Unix. A Linux laptop can carry out many powerful operations on these large data sets using the same tools that super computer users run under Unix.

One of the differences between Linux and Windows users is that Linux users naturally tend to have a naturally clearer idea about how information is being held on their computers. While the Ubuntu GUI can do almost anything the Windows GUI can do Linux users tend to use the command line at least occasionally and are more comfortable editing configuration files.

Thus when work arises that involves moving around and reformatting large amounts of data Linux is a natural platform.

The tools for working with GRIB and netCDF have not been designed with ”user friendly” GUIs. However they fit easily into a scripting framework. At the time of writing the tool that appears to have the best documentation, ease of to installation and general simplicity of use is CDO (Climate Data Operators). The program can be downloaded as source code from

http://www.mpimet.mpg.de/fileadmin/software/cdo/.

This software offers a very strait-forward, powerful command line interface that can be used for a large number of pre-processing and processing tasks. Because both grib and netcdf file formats are in common use it is important to compile the source with netcdf support.

In Ubuntu or Debian install netcdf-bin and the libraries either using Synaptic or by typing

sudo apt-get install netcdf-bin

sudo apt-get install netcdfg-dev

sudo apt-get install libnetcdf3

This installs the programs ncdump and ncgen which convert NetCDF files to ASCII and back, respectively. NetCDF (network Common Data Form) is an interface for scientific data access and a freely-distributed software library that provides an implementation of the interface. The netCDF library also defines a machine-independent format for representing scientific data. Together, the interface, library, and format support the creation, access, and sharing of scientific data.

Now cdo can be built with netcdf support. Check the location of the netcdf header with

locate netcdf.h

In Ubuntu the header is found in /usr/include. So compile the cdo source code and install it with the following option.

sudo ./configure –with-netcdf=/usr/include

sudo make

sudo make install

Now if you get output after typing “cdo” in the console the program should be correctly installed. Full documentation and a reference card cdo_refcard.pdf can be found in the doc directory.

Next install the ncdf package for R. This can apparently be installed under Windows for R2.2.+ and includes all that is needed to run. However Windows users would have to use Cygwin to compile and run cdo, thus their lives would be much more difficult. Note that under Linux you must install the netcdf binary and headers before installing the R package. More details can be found here.

In the next post on the topic I will show how to begin analysing output from a major GCM in R.

Hay razones ideológicas de apoyar el movimiento de software libre. Tengo simpatía con ellos. Pero uso Linux (Ubuntu) en mi laptop por razones puramente prácticas.
Un sistema operativo no es más que una plataforma para correr programas. Como tal, mientras todos los programas corren bien, uno no debe pensar sobre el asunto. Este es especialmente importante cuando se trata de la maquina que usas todo el tiempo. El problema con Windows es que simplemente no hace lo que prometa hacer.
1. No quiero arriesgar la perdida de tiempo asociado con infecciones de virus.
2. No quiero spyware.
3. No quiero estar constantemente checando todos los archivos que recibo por virus
4. No quiero tener que reiniciar la maquina varios veces al día si abro mas de tres programas para recuperar la memoria (Windows XP).
5. No quiero tener que reinstalar el sistema operativa cada seis meses para mantenerlo funcionando a la velocidad inicial (Windows XP)
6. No quiero tener un sistema instalada que no puedo fácilmente reinstalar (Windows Vista)
7. No quiero que el sistema operativa decide cuales programas puedo instalar (Windows Vista)
8. No quiero comprar mas de 1GB de memoria para correr el sistema operativa (Windows Vista)
9. No quiero esperar constantemente para abrir programas.
10. No quiero pensar que mi laptop ha costado mil pesos mas que vale por haber sido forzado a comprar el sistema operativa al mismo tiempo (gracias a Dell este no es siempre cierto)
11. Quiero un diseño con interface moderno , atractivo y fácil de usar (Ubuntu)
12. Quiero tener la confianza que siempre puedo reinstalar el sistema rápidamente si pasa algo.
13. Quiero correr programas gratis y abiertos como Open Office, GRASS, Qgis, Gimp en una forma eficiente. (Nota: Los versiones de los mismos para Windows son disponibles pero no corren tan rápidamente)
14. Quiero instalar programas en línea en vez de buscar discos de instalación que se puede perder.
15. Quiero tener disponible herramientas practicas para procesar datos y imágenes usando “scripts” (Image Magick, Awk, Python, wget etc)
16. Quiero configurar mi laptop como servidor (Apache2, MySQL etc)
17. Quiero usar bases de datos potentes y geográficamente explicitas (Postgresql)
18. Quiero usar la memoria eficientemente para que puedo correr otros sistemas operativas como maquinas virtuales si un programa que tengo que usar no es disponible para el sistema nativa (e.g Word, Excel etc).
19. Quiero tener varios espacios de trabajo (Vista si tiene, XP no)
20. Quiero divertirme de vez en cuando haciendo el sistema visualmente atractiva (Beryl o Compiz )
21. Quiero encontrar información para arreglar mis propios problemas si el sistema no funciona.

Quiero que la pantalla que esta en frente de mis ojos tantas horas por dia se ve asi …….