Happy PostGIS day!
Posted on 17 November 2011 | No responses
Happy PostGIS day everyone!
Until about year ago, most of my work involved using proprietary GIS software, most of it ArcGIS. Today, almost 100 % of the work i do involve open-source software. Contributing much to this shift has without a doubt been PostGIS.
One of my main projects the last year have been the development of the PRIO-GRID project for the peace and conflict research arena. PRIO-GRID is a spatio-temporal grid structure constructed to aid the compilation, management, and analysis of spatial data within a time-consistent framework. It consists of quadratic grid cells that jointly cover all terrestrial areas of the world. Although PRIO-GRID was designed with peace and conflict research in mind, its potential applicability extends well beyond the study of civil war. What it really is about is to collect spatial data from diverse sources and combine these into a common unit of observation.
The first version of the PRIO-GRID was developed using ArcGIS. We quickly realized that for scientific purposes, ArcGIS did not cover the requirements we had for data replication, easy documentation of all processes we did and software capabilities. In mid 2010 we started discovering PostGIS, and quickly realized that we could create far more advanced queries that covered all of our needs. One example of this was when we wanted to calculate the distance from a set of centroids in the grid to the nearest neighbouring polygon belonging to a different country but still “touching” the departure point’s respective polygon.
In ArcGIS we would have to do this one polygon pair at the time.
In PostGIS we simply wrote this as:
SELECT a.gid, a.gwcode, a.gridyear, MIN(ST_Distance(a.centroid, b.geom))/1000 INTO borddist FROM gridcells a, countrypolygons b, countrypolygons c WHERE a.gwcode != b.gwcode AND a.gwcode = c.gwcode and ST_Intersects(b.geom, c.geom) GROUP BY a.gid, a.gwcode, a.gridyear;
All in all. What i really would like to say is that PostGIS have helped us save time, money and some potential grey hair.
Happy PostGIS day!
PostGIS How-to: NetCDF to PostGIS
Posted on 15 November 2011 | 3 responses
Disclaimer: I am not sure if this is the quickest way of importing NetCDF files to PostGIS, but it works, and for me that is the most important. If you have any other suggestions for improvements, please post them in the comment field!
In my work we need to combine a lot of spatial data from diverse sources. This includes socio-economic, physical and climatic data. Our latest requirement is to be able to import NetCDF files containing climate model data into our PostGIS database. NetCDF (Network Common Data Format) is a multidimensional spatio-temporal data structure with the usual x, y spatial coordinates and a z coordinate or value in addition to the t variable which represents time since some start date.
Our approach is to read the NetCDF in Python using the Scipy.io NetCDF module. This makes it possible to read the NetCDF file into Python and then output it to PostGIS using the psycopg2 module.
Read more
PostGIS How to: Calculate the average raster pixel value inside a polygon
Posted on 2 November 2011 | No responses
In my daily work I use both vector and raster data and the combination of both. Often, I want to include raster data into our common vector grid. Hence, I often want averaged raster pixel values inside each of the grid cells in our vector grid.
To do this, we need to use both the vector and raster functionality in PostGIS.
In this example i want to calculate the average infant mortality rate pixel value inside each vector grid cell.
One query we use for this is:
SELECT gid, CAST(AVG(((foo.geomval).val)) AS decimal(9,3)) as avgimr FROM (SELECT p.gid, ST_Intersection(r.rast, p.cell) AS geomval FROM imrgrid2p5m r, priogrid_land p WHERE ST_Intersects(p.cell, r.rast)) AS foo WHERE (foo.geomval).val >=0 GROUP BY gid ORDER BY gid;
To explain, we select the vector grid cell’s gid, and the average of (foo.geomval).val FROM the Intersection of the infant mortality raster and the vector grid. Hence, the geomval is a list of each pixels intersection with the vector grid cell. This geomval contain both the actual geometry of this intersection and the value of the pixel. The we apply a where clause of where the cell and raster intersects since the raster is tiled. We also do not want to include negative values (missings) since infant mortality rate is not negative. Then we group the query by the grid cell’s gid since we are aggregating the average pixel value.
PostGIS raster values become integer when Polygonizing
Posted on 1 November 2011 | No responses
Lately I have been trying to find out why my raster decimal values are truncated to integer when polygonizing.
After some research I found out that this is a known bug (http://trac.osgeo.org/postgis/ticket/650), and the problem lies in the GDALPolygonize function. This function is used in PostGIS to convert raster to vector, and since this function do not support decimal values, the values are truncated to integers. Particularly problematic if you are working with a raster with pixel values between min 0 and max 1.
The solution is to use GDAL 1.9dev Subversion from trunk. Hence, you have to compile GDAL 1.9dev yourself, and then recompile PostGIS 2.0.0SVN to include GDAL 1.9dev.
Then, when querying raster to vector with values, you get the correct result.
The below query tests this by querying the average pixel values within a quadratic polygon of 0.5 x 0.5 decimal degree’s.
SELECT gid, AVG((foo.geomval).val)
FROM (SELECT p.gid, ST_Intersection(ST_SetSRID(p.cell,4326), ST_SetSRID(r.rast,4326), 1) AS geomval FROM mountain r, priogrid_land p
WHERE ST_Intersects(ST_SetSRID(p.cell,4326), ST_SetSRID(r.rast,4326),1)) AS foo
WHERE gid = 183230 AND (foo.geomval).val >= 0
GROUP BY gid
ORDER BY gid;
Result
gid;avg
183230;0.404999998294645
PostGIS Back-up and restore
Posted on 4 October 2011 | No responses
Ever wondered on what is the most appropriate way of backing up and restoring your precious PostGIS database? Or have you experienced errors when restoring your old database backup file?
The solution is, put simple: Do not store your data in the public schema in PostGreSQL!
Paul Ramsey explains more in detail on his blog:
http://blog.cleverelephant.ca/2010/09/postgis-back-up-restore.html
Soon back on track
Posted on 18 May 2010 | No responses
Hey all visitors. Sorry for not updating the blog for a while. I have been very busy working on my thesis.
Updates are coming in a week!
Bounding Containers
Posted on 4 March 2010 | No responses
I had to make some minimum bounding circles from a set of points in ArcGIS. After some searching i came across a very useful script by Dan Patterson. Actually it comes as a set of scripts in a toolbox. Ready to be imported into ArcGIS 9.3.
The toolbox has scripts for generating convex hull, minimum area bounding rectangle, minimum area circle and extent polygons.
You find it here.
Natural Earth – Free Raster and Vector Data
Posted on 17 February 2010 | 1 response
Somewhat unfamiliar to me is the Natural Earth project, providing free raster and vector data of administrative borders, disputed areas, populated places, urban polygons and water boundaries in addition to much more.
Great for updated maps. For historical shapefiles. Check my post on the cShapes dataset.
From their website:
“We developed Natural Earth as a convenient resource for making custom maps. Unlike other map data intended for scientific analysis or military mapping, Natural Earth is designed to meet the needs of production cartographers using a variety of software applications. Maximum flexibility is a goal.”
Go to Natural Earth Website
Reference: http://hansvandermaarel.wordpress.com/
Population Estimate Service
Posted on 16 February 2010 | No responses
Many within the GIS community working on global data are familiar with the Gridded Population of the World (GPW). This global gridded population estimate has now been more available through a web-service. This lets the user zoom around to investigate the data without having to download the rasterdata.
Additionally and more interesting are the population estimation service, that can be accessed through OGC, REST and SOAP interface to estimate the population within a drawn polygon.
One of the new mapping tools also released, based on the technology used by Google Maps, demonstrates the Population Estimation Service. It lets users select an area of interest by drawing a polygon on the map and submit the request to the service, and it displays the results. (CIESIN News)
The google map based Population Estimation Service can be accessed here.
FME – Raster to polygon
Posted on 15 February 2010 | No responses
The last couple of weeks i have been working on a raster to vector conversion in FME, and thought it would be nice to share my experiences here at my blog. This might be the start to a FME special series here at the blog with short explanations of different topics within FME. To read this raster to vector special click the read more link below.
What i wanted to do was to create polygons representing each of the rastercells, and extract the value from the raster and copy this value to the attribute table of the polygon.
I started off importing the raster to my workspace. One could also do this as a batch operation, importing multiple rasters.
To convert the raster cell to polygon i used the the ‘RasterCellCoercer’. This transformer decomposes each of the raster features into individual points or polygons. However, the transformer does not extract the value, even when enabling the “preserve attributes” option.
To be able to extract the value from each of the raster cells into an attribute i used the transformer ‘ElevationExtractor’. This tool extracts the z-value from each of the raster cells into a defined attribute. The default is “_elevation”.
That is basically it. You can now add your output writer and add the user attribute you want to store the raster cell value from the “_elevation” attribute. Remember that your output feature must be equal to either the point or polygon option that you selected in the ‘RasterCellCoercer’.
If you have selected multiple rasters, you could write each out to its original filename using fanout by “fme_basename”.
Here is my workspace:
