Archive for the ‘ PostGIS ’ Category

Postgis 2.0.2 hits Arch repository

Arch has posted the PostGIS 2.0.2 release that came out December 3rd to the Community repository. After Pacman updates the package Postgres needs to be made aware of the upgrade.  Use:

ALTER EXTENSION postgis UPDATE TO "2.0.2"; 
ALTER EXTENSION postgis_topology TO "2.0.2";

Chernoff Faces of Bike Rides

Chernoff faces of 4 bike rides

Chernoff faces from 4 bike rides

I plugged some of the fit file data that I stored in a Postgres database into R to generate the image above.  Two R libraries, RpostgreSQL to connect to Postgres, and aplpack, to generate the faces were used.  Once the libraries were loaded the following R command pulls the data out of Postgres:

rs <- dbSendQuery(con, "
SELECT ride_date
, ST_Length(ST_Transform(ST_Makeline(the_geom ORDER BY ride_time),26916)) As length
, date_part('hour', max(ride_time))-date_part('hour', min(ride_time)) As ridetime
, avg(temperature) As temp, avg(speed) As avg_speed
, avg(altitude) As alt, max(altitude)-min(altitude) As alt_diff
, avg(cadence) As rpm 
FROM ride GROUP BY ride_date ORDER BY ride_date
")

Then fetch the data from the record set with:

rides<-fetch(rs, n=-1)

And finally plot the faces with:

faces(rides[,2:8], face.type=0, labels=rides$ride_date)

The faces function draws the faces based on the order of the variables.  The features are modified in this order:

height of face 
width of face 
structure of face 
height of mouth 
width of mouth 
smiling 
height of eyes 
width of eyes 
height of hair 
width of hair 
style of hair 
height of nose 
width of nose 
width of ear 
height of ear

If you don’t have enough variables to fill this list it will wrap around and start again from the top.  For more configuration options use ?faces after you load the aplpack.

There are two small problems with the data, and they both caused by stopping the timer on your bike computer, then restarting at a later time and place.   The time ridden is calculated here by simply subtracting the minimum time from the maximum.  If you rode for 3 hours and took an hour break the query would return a time of 4 hours.  Similarly, if you ride for a mile, hop in a truck and drive 10, then ride for another mile you’d get a result of 12 miles.  There’s a way to work around this, I just haven’t figured it out yet.  I suppose the best solution is to harden up and not take breaks.

The data is stored as points in the database, to get the length ST_Makeline “connects the dots” from point to point.  An awesome feature of Postgres 9.0+ is the Order By in the middle an aggregate function.  It helped in this case because the first ride graphed doubled up and crossed over itself several times.  This lead to the query planner making some interesting decisions on where the line should go.  Forcing the process to follow the points in order the line followed my route perfectly.

Ranking the Indiana Breweries by Population

In my last post I used PostGIS to calculate the number of Hoosiers that live within 5 miles of an Indiana microbrewery.  Now I want to rank each brewery by population.  The total is going to be a little more than the 2.5 million-ish I came up in the last post due to the lucky duckies who live with 5 miles of more than one brewery.

Flat 12 of Indianapolis, with over 243,000 people within 5 miles of their location is our winner.  No big surprise, but the top 11, and 12 of the top 13 are all Indy breweries.

A few breweries, most notably Three Floyds in Munster are lower than they should be because I only used Indiana census data.  Three Floyds is only a half mile from the Indiana / Illinois border.  If I included Illinois data the nearby Chicago population would have pushed them up the list.

The full list and the SQL used to create it are below the jump.

Continue reading

How many Hoosiers live with 5 miles of an Indiana brewery?

2,539,798.  Approximately.  I think.  Maybe.  Here’s what I did to get there.

First, the Census department has put together some nice pre-digested block layers, they include the total population for the blocks so you don’t have to compile that data yourself.

Second, create a 5 mile (acutally I used 8046.72 meters because my data in a UTM projection) buffer around my brewery points.  Shp2pgsql, use it to shove the buffers into a PostGIS database and call the layer buffer_5mile.  Fun!

Third, load the block data into QGIS and add a numeric field (I used totarea) to store the area of the blocks.  Use shp2pgsql to reproject the census data to match your buffers and load the results into PostGIS.  I called my layer block_pop.

Fourth, update that new area field with

UPDATE block_pop SET totarea = st_area(geom);

Fifth, to get an estimate of the population within the buffers use:

SELECT 
round(sum(c.pop10*(ST_Area(ST_Intersection(c.geom, b.the_geom))) / c.totarea)) 
FROM
(SELECT ST_Union(geom) as the_geom FROM buffer_5mile) as b, block_pop c
WHERE ST_Intersects(the_geom, c.geom);

Bam.  You done.

Updated Indiana Breweries Map

A few months ago I followed Darren Cope’s Geography of Tim Horton’s blog post to create a similar map that showed breweries in Indiana.  A new brewery opened up down south so I decided to update the map.

To get the location I looked up the location on Google Maps then used the QGIS plugin OpenLayers to find the brewery on Google Maps’ Hybrid Map.  I added the point to my breweries PostGIS layer and re-imported the layer into GRASS along with a shapefile of the state.  From there I used v.to.rast.constant to convert those layers to rasters.

With the vectors converted to rasters I used r.mask to use the Indiana raster as a mask for the next step which is to use r.grow.distance to create the final map.  The result was exported with r.out.png and brought into GIMP to make the no data areas transparent.

To create the transparencies use Layer -> Transparency -> Add Alpha Channel then use the Select Region by Color tool to select the no data areas, then press the delete key to zap them away.  The image below is the result.

 

Location of breweries in Indiana

Indiana breweries

Views, PostGIS 2, and QGIS

If you wanted to create a view with spatial capabilities with previous versions of PostGIS all you needed to do was manually register the view in the geometry_columns table.  This is no longer possible since from PostGIS 2.0 forward geometry_columns is a view, not a table.  However, if you use typemods to define the geometry column your view the geometry_columns view will pick it up.  For example, I wanted to take some of the data I loaded from my .fit files and convert the metric data to imperial units, and reproject the data from WGS84 to UTM.   Here’s the view that does that for me:

CREATE OR REPLACE VIEW rides.ride_utm AS
SELECT id, no
, utility.units_from_to('meter', 'feet', altitude) AS edge_elev
, 32+(temperature*(1.8)) AS temp_f, speed*2.2369 AS mph, cadence
, ride_date, ride_time, elev_dem
, ST_Transform(the_geom, 26916)::geometry(Point, 26916) As the_geom FROM rides.ride;

The units_from_to function is something I picked up from PostGIS in Action.  The geometry column uses ST_Transform to reproject the data, and the typemod after defines the column as having point geometry with my UTM projection.

With the geometry column defined in that manner the geometry_columns view is properly populated.  When you want to use this view in QGIS it is visible as a loadable layer.  However, before you can use it you must tell QGIS which field is the unique row identifier.  To the left of the layer’s SRID number there will be a drop box on the PostGIS layer load screen.  This dropbox will give you a list of possible fields to use as your ID column.  Just select the column that contains unique values and you’re good to go.

PostGIS in Action 2nd edition is underway!

A few weeks ago Regina Obe posted that work has begun on the 2nd edition of PostGIS in Action.  I can’t wait until it is released, the first edition is an amazing resource and I’m sure the new edition will be better at covering topics covered in the first edition, plus be the best resource for all the new stuff that PostGIS 2.0 brought us such as topology and raster support.

Understandably it’s going to take some time to put all this information together, so until the new edition drops keep an eye on the BostonGIS page, the BostonGIS blog, and everything else Regina Obe and Leo Hsu do in the world of PostGIS and Postgresql.  It’s really too much to keep track of!