Welcome to Step Five in our exploration of open source GIS using Ubuntu. In the previous installment, we uploaded some spatial data into a PostGIS. database. The next logical step is to install QGIS and add start visually analyzing this data, which we’ll get to in Step Six. But before that, I wanted to take a short time out and give a quick demonstration of some of the spatial analysis functions PostGIS can perform on its own through SQL database queries.
Warning: If you haven’t completed the previous steps, installing and configuring PostgreSQL and PostGIS, and uploading Portland streets data, none of the following will work at all – you’ll get all kinds of crazy errors. So unless you have these programs and datasets installed and working on your machine already, you’ll need to go back and do that before proceeding.
1. Connect to the Database (using the terminal):
psql -d portland
In the previous step, we added Portland’s custom coordinate system definition (SRID 6856) to the spatial_ref_sys database. However, that definition was originally missing a piece of information – the proj4 text. Proj4 is a reprojection system used by PostGIS to transform coordinates. If we add the proj4 text, we’ll be able to project coordinates in the terminal. I went back and updated the previous step to include the proj4 text, but if you missed that, here’s the instructions to do it now:
UPDATE spatial_ref_sys SET proj4text=’+proj=lcc +lat_1=44.33333333333334 +lat_2=46.0 +lat_0=43.66666666666666 +lon_0=-120.5 +x_0=2500000 +y_0=0 +ellps=GRS80 +units=ft +no_defs’ WHERE srid=6856;
UPDATE 1
2. Simple Selection
The first thing to try is a basic selection. Here’s the code to select SE 4th Street:
SELECT prefix, streetname, ftype FROM public.streets WHERE streetname=’4TH’ AND prefix=’SE’ AND ftype=’ST’;
prefix | streetname | ftype
--------+------------+-------
SE | 4TH | ST
SE | 4TH | ST
SE | 4TH | ST
SE | 4TH | ST
SE | 4TH | ST
SE | 4TH | ST
SE | 4TH | ST
SE | 4TH | ST
SE | 4TH | ST
SE | 4TH | ST
SE | 4TH | ST
SE | 4TH | ST
SE | 4TH | ST
SE | 4TH | ST
SE | 4TH | ST
SE | 4TH | ST
SE | 4TH | ST
SE | 4TH | ST
(18 rows)
There is only one SE 4th Street, but it is stored in the database as 18 segments, extending from intersection to intersection. This makes it easier to do network analysis and geocoding.
3. Finding an Intersection
Since each street segment has its own database entry, every two street cross intersection has four intersection points, all of which have the same coordinates. And it’s likely that each street will include ten or more segments. This requires a complicated SQL statement that will Select each street as a subtable, and then Select the intersection from the subtables. For this example, we’ll be selecting the intersection of NW Lovejoy Street and NW 23rd Avenue.
SELECT distinct(AsText(intersection(b.geom, a.geom))) AS the_intersection FROM (SELECT geom FROM public.streets WHERE streetname=’LOVEJOY’ AND prefix=’NW’ AND ftype=’ST’) a, (SELECT geom FROM public.streets WHERE streetname=’23RD’ AND prefix=’NW’ AND ftype=’AVE’) b WHERE a.geom && b.geom AND intersects (b.geom, a.geom);
the_intersection
------------------------------------------
POINT(7638738.59481627 686914.222769022)
(1 row)
The results are in projection units – feet. That might be useful in coordinate space, but perhaps latitude and longitude might be more useful. So here’s the same query again, with an added transformation function that reprojects to SRID 4326, which is unprojected WGS84 coordinates.
SELECT distinct(AsText(ST_Transform(intersection(b.geom, a.geom),4326))) AS the_intersection FROM (SELECT geom FROM public.streets WHERE streetname=’LOVEJOY’ AND prefix=’NW’ AND ftype=’ST’) a, (SELECT geom FROM public.streets WHERE streetname=’23RD’ AND prefix=’NW’ AND ftype=’AVE’) b WHERE a.geom && b.geom AND intersects (b.geom, a.geom);
the_intersection
-------------------------------------------
POINT(-122.698608851154 45.5297891204646)
(1 row)
You can check this result by transposing the numbers (so that the latitude is first) and then pasting them into the search bar in Google Maps:
45.5297891204646, -122.698608851154
4. Geocoding
Each street segment in the database contains 4 numbers that describe the address ranges used by the segment. Leftadd1 is the lowest number on the left, and rgtadd1 is the lowest number on the right; while leftadd2 and rgtadd2 are the highest numbers. To geocode a street address, the first step is to find the street segment containing the address. As an example, let’s search for the House of Ramen restaurant in downtown Portland, which is located at 223 SW Columbia Street.
SELECT gid, prefix, streetname, ftype, leftadd1, leftadd2, rgtadd1, rgtadd2 FROM public.streets WHERE streetname=’COLUMBIA’ AND prefix=’SW’ AND ((leftadd1 < 223) OR (rgtadd1 < 223)) AND ((leftadd2 > 223) OR (rgtadd2 > 223));
gid | prefix | streetname | ftype | leftadd1 | leftadd2 | rgtadd1 | rgtadd2
-------+--------+------------+-------+----------+----------+---------+---------
40183 | SW | COLUMBIA | ST | 200 | 298 | 201 | 299
(1 row)
Based on this info, it looks like our odd-numbered address will be on the right side of the street. Take note of the gid, the unique number we can use later to isolate this segment. But first, we have to determine how far down the street 223 is likely to be. We could estimate that it will be about 22% of the way down, and that’s probably close enough. But there’s a mathematical way to get a more precise value.
proportion = (number – rgtadd1)/(rgtadd2 – rgtadd1) = (223-201)/(299-201) = 22/98 = 0.2245
Therefore our query will be searching for the coordinates of the point that is 22.45% of the way down street segment 40183:
SELECT ST_AsText(ST_Line_Interpolate_Point(geom, 0.2245)) FROM public.streets WHERE gid=40183;
st_astext
------------------------------------------
POINT(7644036.17590994 680787.382718017)
(1 row)
Again, this is in projected units, so let’s use an on-the-fly transformation to get an answer in decimal degrees:
SELECT ST_AsText(ST_Transform(ST_Line_Interpolate_Point(geom, 0.2245),4326)) FROM public.streets WHERE gid=40183;
st_astext
-------------------------------------------
POINT(-122.677299994613 45.5133847845873)
(1 row)
If you transpose these results and search for them in Google Maps, you will see some interesting results. There will be a green arrow on the map showing a point that is exactly 22.45% of the way down the street. However, it shows House of Ramen being much further down the street. Google Maps also assumes that you’re looking for something else, perhaps in the building on the other side of the street, because its search algorithm is trying to anticipate your needs by choosing the nearest business location. What we’ve done is mathematically accurate, but it doesn’t actually reflect the real-life addressing of SW Columbia Street.
Well, that’s just a quick demonstration of how SQL can be used to do spatial analysis on a PostGIS table. This may not be helpful for desktop analysis, but it is an important part of developing spatially enabled web applications. In Step Six, we’ll be installing QGIS, and that’s where the real fun begins.
Comments are closed.
