Archive for the ‘ PostGIS ’ Category

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:

round(sum(c.pop10*(ST_Area(ST_Intersection(c.geom, b.the_geom))) / c.totarea)) 
(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 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:

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!

Arch updates Postgresql & PostGIS

Postgresql 9.2.1 and PostGIS 2.0.1 just hit Arch’s repositories.  Going from Postgres 9.1 to 9.2 is a major upgrade and will require your 9.1 databases to be migrated.  For steps on how to upgrade your data see this post in the Arch Wiki.  If you want to postpone updating your databases add postgresql and postgresql-libs to your IgnorePkg line in your /etc/pacman.conf file.

Getting .fit file data into PostGIS

GPSBabel released an update (1.4.4) on Labor Day.  This update fixed GPSBabel’s FIT file compatibility.  Now FIT files created on Garmin Edge and other devices can be converted into formats that can be imported into PostGIS.  Follow along after the jump to see how.

Continue reading

Temporary PostGIS-2.0.1 PKGBUILD

A few months ago I put out a PKGBUILD for PostGIS 2.0 to use while we waited for the official Arch repositories to catch up.  I just updated it to 2.0.1.  This time I changed it so it will build a package called Postgis-2.0.1-0, so when the official repositories add 2.0.1 it will be replaced automatically.  If you currently have 2.0.0 installed you can use this to upgrade but you’ll still have to upgrade your Postgresql extensions by using:

ALTER EXTENSION postgis_topology UPDATE TO "2.0.1";

As always, I only kind of know what I’m doing, use at your own risk.

PostGIS 2.0 in Arch’s official repository!

A  little over a week ago PostGIS 2.0 hit the official Arch repository.  Hooray!


Trying to improve the TIGER Geocoder query for QGIS

In an earlier post I showed how you can use the Fast SQL Layers plugin (available from the QGIS Plugin Repository) to execute a TIGER Geocoder query.  I’ve simplified the query slightly, instead of using a UUID for the unique ID field you can use a window query and the row_number() function instead.

SELECT ROW_NUMBER() OVER (ORDER BY g.rating) as id, 
 g.rating, ST_X(g.geomout) As lon, ST_Y(g.geomout) As lat, 
 (addy).address As stno, (addy).streetname As street, 
 (addy).streettypeabbrev As styp, (addy).location As city, (addy).stateabbrev As st,(addy).zip, g.geomout As the_geom
 FROM geocode('5200 Duesner Rd, Evansville, IN 47712') As g;

While simpler, it doesn’t quite work as intended.  If you run the query in Pgadmin or psql it works just fine.  However, running it in Fast SQL Layers the window doesn’t work, the ID for each feature is zero.  This causes a crash when you try to open the attribute table.  The identify tool and exporting to shapefile both work, but as mentioned the ID field will show 0 for all records returned.

How far away are you from beer on a Sunday in Indiana?

Answer:  At most 50 miles.

Inspired by GIS blogger Darren Cope’s analysis of Tim Horton’s locations in Canada I decided to figure out how far you are at any given time from a microbrewery in Indiana.   Instead of Canadian coffee and donuts I mapped brewery locations.  I pulled the Brewer’s Guild of Indiana’s membership list (and took out Bee Creek, I’m pretty sure they’re closed).  This is a vital map because Indiana’s blue laws prohibit alcohol sales on Sundays except at microbreweries and vineyards.   Two things I don’t like on Sundays, wine and crossing state boundaries so this was a vital project for me.  I found that I really don’t want to be in northeast Dubois County on a Sunday. Seriously, 50 miles.  Dang.  At least people in Vincennes and Terre Haute can jump over to Illinois for a 6 pack.

Distances to Indiana breweries

I followed Cope’s steps except to get the lat / long for the breweries I had to geocode them against the TIGER street centerlines and instead of using GRASS directly I used QGIS’s GRASS plugin.  This article was a big help on getting a GRASS workspace set up to create the raster.