Smathermather's Weblog

Remote Sensing, GIS, Ecology, and Oddball Techniques

Posts Tagged ‘Contours’

Proper (ab)use of a database, contour interpolation using #postgresql #postgis #arcgis

Posted by smathermather on August 6, 2012

Anyone who has been following along at home knows I don’t think much like a DBA.  Sometimes that’s good; mostly it’s probably bad.  In this post, I hope it will be interesting.

The problem of the day is how to take engineering contours derived from breaklines, a lidar point cloud, and all the lot, and do a good job interpolating that to a DEM.  This is an interesting problem, as it’s another one of those problems with so many sub-optimal solutions.  In part, it’s an impossible problem– going from a surface model to contours is inherently a loss of information, so going the other way is difficult to do believably, both with respect to believable smoothness and the requisite feature retention.  This is a little like the pugilistic pie puppets problem: ‘flaky vs. tender’ of making good pie crust.  Doing one well seems to preclude doing the other well.

Now, Chris Hormann, of PovRay fame has a solution, there are solutions in ArcGIS, and in GRASS, and plenty of other attempts.  Here’s a decent list, which I haven’t checked for timeliness:  http://vterrain.org/Elevation/contour.html.  Since we have a ArcGIS license, it handles the dataset I have OK, the results checked out OK, this is what we chose to use for now.  I tried GRASS, but it was too slow for the amount of data being processed.  The command in ArcGIS ends up looking something like this:


TopoToRaster_sa "N2260635.shp elevation Contour;N2260640.shp elevation Contour;N2260645.shp elevation Contour;N2265633.shp elevation Contour;N2265640.shp elevation Contour;N2265645.shp elevation Contour;N2270635.shp elevation Contour;N2270640.shp elevation Contour;N2270645.shp elevation Contour;N2275650.shp elevation Contour" "C:\Temp\TopoToR_N2269.tif" 2 "2265031 645070 2270065 639945" 20 # # ENFORCE CONTOUR 40 # 1 0 # # # # # #

And thus we go from this:

To this:

Now to do it 623 more times… .  And in walks PostgreSQL for some abuse as a spatial engine with string manipulation to boot.

Each tile of contours I run through this approach needs the adjacent 8 tiles of contours loaded in order to avoid dreaded edge effects:

Except that sometimes, there aren’t 8 tiles available, like on the edge of large water bodies:

So, we can’t just use a heuristic that loads the 8 adjacent tiles– we need to enumerate these tiles.  An ST_Touches should do the trick here:

SELECT p.tile, a.tile AS tiles
FROM cuy_contour_tiles AS p, cuy_contour_tiles As a
WHERE ST_Touches(p.geom, a.geom);

But, that does not arrange our table very nicely:

I want something that looks like this:

N2225655, N2225660, N2220650 …

N2225660, N2225665, N220665 … etc.

Thank goodness for Postgres’ String_agg, and moreover for Leo Hsu and Regina Obe’s blog.

Adapting their entry, we can reorganize those entries like this:

SELECT
STRING_AGG(a.tile, ',' ORDER BY a.tile)
As tiles
FROM cuy_contour_tiles AS p, cuy_contour_tiles As a
WHERE ST_Touches(p.geom, a.geom)
GROUP BY p.tile
ORDER BY p.tile;

Ah, now that’s better:

And if we want, we can use an arbitrary delimiter, and build our entire command in one swell foop:

SELECT p.tile,
'TopoToRaster_sa "'
|| 'W:\contours\cuy_contours\shapefiles\' || p.tile || '.shp elevation Contour;'
|| 'W:\contours\cuy_contours\shapefiles\'
|| STRING_AGG(a.tile, '.shp elevation Contour;W:\contours\cuy_contours\shapefiles\' ORDER BY a.tile)
|| '.shp elevation Contour;" '
|| '"v:\svm\dem\cuyahoga\' || p.tile || '.tif" 2'
|| ' "' || replace(replace(replace(ST_Extent(p.geom)::text, 'BOX(', ''), ')', ''), ',', ' ') || '"'
|| ' 20 # # ENFORCE CONTOUR 40 # 1 0 # # # # '
As tiles
FROM cuy_contour_tiles AS p, cuy_contour_tiles As a
WHERE ST_Touches(p.geom, a.geom)
GROUP BY p.tile
ORDER BY p.tile

Ok.  I have 623 more tiles to do.  Until next time… .

Posted in Analysis, ArcGIS, Database, PostGIS, PostgreSQL, SQL | Tagged: , , , , , , , , , , | Leave a Comment »

Contours– Symbolized from a single table

Posted by smathermather on July 15, 2011

In a previous post, I restructured the contour data for display in GeoServer, e.g.:


UPDATE base.contours_2
	SET 	div_10 = CAST( contours_2.elevation % 10 AS BOOLEAN ),
		div_20 = CAST( contours_2.elevation % 20 AS BOOLEAN ),
		div_50 = CAST( contours_2.elevation % 50 AS BOOLEAN ),
		div_100 = CAST( contours_2.elevation % 100 AS BOOLEAN ),
		div_250 = CAST( contours_2.elevation % 250 AS BOOLEAN );

The SLD styling for the contours based on the new data structure will be forthcoming with my intern’s upcoming guest post, but in the meantime, I have a teaser of a map snapshot, served up through a GeoExt interface (note the correct contour intervals showing in the legend simply and elegantly because of the data structure and a single SLD):

Contour Map, 10ft and 50ft contours

Over-zoom showing 2ft and 10ft Contours

Posted in Database, Database Optimization, GeoServer, PostGIS | Tagged: , , , , , , , | Leave a Comment »

Contours– Structuring PostGIS data for viewing with GeoServer

Posted by smathermather on May 25, 2011

Naively structured data is my bane– the desire (and need) to get stuff done so often overtakes the time needed to do things the better way. So, we bootstrap.

A long time ago, we managed to load in a few tens of gigs of contour data into PostGIS, partitioned it into 2ft, 10ft, 20ft, 50ft, 100ft and 250ft tables using select queries with a modulus operator, e.g.


CREATE TABLE base.cuy_contours_10
	AS
	SELECT elevation, the_geom
		FROM base.cuy_contours_2
		WHERE base.cuy_contours_2.elevation % 10 = 0;

And then we built SLD‘s to display the data in GeoServer and formed the data into layer groups, and away we went… .

However, layer groups can often be replaced by properly structured PostGIS tables. For large (and homogenous) datasets like this, it’s the only way to go. In addition, properly structuring this allows us to take advantage of GeoServer’s ability to modify the legend based on zoom extent, which for representing contours at a range of scales is an almost “must have”.

To structure the table, we could be disciplined and build this out as a proper normalized relational dataset where our gid is used to determine if a contour is divisible by 10, 20, 50 etc., and while “I don’t wanna” isn’t a good enough reason not to do this, I think the computational overhead of a database view piecing these data back into a single table each time we need to access this would not be justified in the light of the static nature of the table. So database normalization be darned, disk space is cheap, full speed ahead. Let’s add some boolean fields for flagging whether a contour is divisible by our numbers and calculate that out:


UPDATE base.contours_2
	SET div_10 = CAST( contours_2.elevation % 10 AS BOOLEAN );


UPDATE base.contours_2
	SET div_250 = CAST( contours_2.elevation % 250 AS BOOLEAN );

yields the following (with highlights added):

Sort of the opposite of what I intended, the “falses” maybe should be “trues” and vice versa, but hey, it’ll work anyway. BTW, we did test doing this calculation a few different ways, and while I can’t remember all the details, doing this as a calculation instead of doing an update query with a while statement testing the modulus was much faster (3x faster).

Ok, so now we style an sld for it, and sit back and enjoy (pics later…).

Posted in Database, Database Optimization, GeoServer, PostGIS, SQL | Tagged: , , , , , , , , , | Leave a Comment »

SLD for contour data

Posted by smathermather on May 25, 2011

See other post for explanation:

<?xml version="1.0" encoding="ISO-8859-1"?>
<StyledLayerDescriptor version="1.0.0" xmlns="http://www.opengis.net/sld"

xmlns:ogc="http://www.opengis.net/ogc"
xmlns:xlink="http://www.w3.org/1999/xlink" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance"
xsi:schemaLocation="http://www.opengis.net/sld

http://schemas.opengis.net/sld/1.0.0/StyledLayerDescriptor.xsd">
<NamedLayer>
<Name>contours</Name>
<UserStyle>
<Title>contours</Title>
<Abstract>Contour lines with index</Abstract>
<FeatureTypeStyle>

<Rule>

<Name>rule01</Name>
<Title>2 ft contours</Title>
<Abstract>Abstract</Abstract>

<ogc:Filter>
<ogc:PropertyIsEqualTo>
<ogc:PropertyName>div_10</ogc:PropertyName>
<ogc:Literal>1</ogc:Literal>
</ogc:PropertyIsEqualTo>
</ogc:Filter>

<MinScaleDenominator>1</MinScaleDenominator>
<MaxScaleDenominator>4799</MaxScaleDenominator>

<LineSymbolizer>
<Stroke>
<CssParameter name="stroke">
<ogc:Literal>#ffd700</ogc:Literal>
</CssParameter>
<CssParameter name="stroke-linecap">
<ogc:Literal>butt</ogc:Literal>
</CssParameter>
<CssParameter name="stroke-linejoin">
<ogc:Literal>miter</ogc:Literal>
</CssParameter>
<CssParameter name="stroke-width">
<ogc:Literal>1</ogc:Literal>
</CssParameter>

</Stroke>
</LineSymbolizer>
</Rule>

<Rule>

<Name>rule02</Name>
<Title>10 ft contours index</Title>
<Abstract>Abstract</Abstract>

<ogc:Filter>
<ogc:PropertyIsEqualTo>
<ogc:PropertyName>div_10</ogc:PropertyName>
<ogc:Literal>0</ogc:Literal>
</ogc:PropertyIsEqualTo>
</ogc:Filter>

<MinScaleDenominator>1</MinScaleDenominator>
<MaxScaleDenominator>4799</MaxScaleDenominator>

<LineSymbolizer>
<Stroke>
<CssParameter name="stroke">
<ogc:Literal>#D95F02</ogc:Literal>
</CssParameter>
<CssParameter name="stroke-linecap">
<ogc:Literal>butt</ogc:Literal>
</CssParameter>
<CssParameter name="stroke-linejoin">
<ogc:Literal>miter</ogc:Literal>
</CssParameter>
<CssParameter name="stroke-width">
<ogc:Literal>1</ogc:Literal>
</CssParameter>

</Stroke>
</LineSymbolizer>
</Rule>

<Rule>

<Name>rule03</Name>
<Title>10 ft contours</Title>
<Abstract>Abstract</Abstract>

<ogc:Filter>
<ogc:PropertyIsEqualTo>
<ogc:PropertyName>div_10</ogc:PropertyName>
<ogc:Literal>0</ogc:Literal>
</ogc:PropertyIsEqualTo>
</ogc:Filter>

<MinScaleDenominator>4800</MinScaleDenominator>
<MaxScaleDenominator>7199</MaxScaleDenominator>

<LineSymbolizer>
<Stroke>
<CssParameter name="stroke">
<ogc:Literal>#ffd700</ogc:Literal>
</CssParameter>
<CssParameter name="stroke-linecap">
<ogc:Literal>butt</ogc:Literal>
</CssParameter>
<CssParameter name="stroke-linejoin">
<ogc:Literal>miter</ogc:Literal>
</CssParameter>
<CssParameter name="stroke-width">
<ogc:Literal>1</ogc:Literal>
</CssParameter>

</Stroke>
</LineSymbolizer>
</Rule>

<Rule>

<Name>rule04</Name>
<Title>50 foot contours index</Title>
<Abstract>Abstract</Abstract>

<ogc:Filter>
<ogc:PropertyIsEqualTo>
<ogc:PropertyName>div_50</ogc:PropertyName>
<ogc:Literal>0</ogc:Literal>
</ogc:PropertyIsEqualTo>
</ogc:Filter>

<MinScaleDenominator>4800</MinScaleDenominator>
<MaxScaleDenominator>7199</MaxScaleDenominator>

<LineSymbolizer>
<Stroke>
<CssParameter name="stroke">
<ogc:Literal>#D95F02</ogc:Literal>
</CssParameter>
<CssParameter name="stroke-linecap">
<ogc:Literal>butt</ogc:Literal>
</CssParameter>
<CssParameter name="stroke-linejoin">
<ogc:Literal>miter</ogc:Literal>
</CssParameter>
<CssParameter name="stroke-width">
<ogc:Literal>1</ogc:Literal>
</CssParameter>

</Stroke>
</LineSymbolizer>
</Rule>

<Rule>

<Name>rule05</Name>
<Title>20 ft contours</Title>
<Abstract>Abstract</Abstract>

<ogc:Filter>
<ogc:PropertyIsEqualTo>
<ogc:PropertyName>div_20</ogc:PropertyName>
<ogc:Literal>0</ogc:Literal>
</ogc:PropertyIsEqualTo>
</ogc:Filter>

<MinScaleDenominator>7200</MinScaleDenominator>
<MaxScaleDenominator>20999</MaxScaleDenominator>

<LineSymbolizer>
<Stroke>
<CssParameter name="stroke">
<ogc:Literal>#ffd700</ogc:Literal>
</CssParameter>
<CssParameter name="stroke-linecap">
<ogc:Literal>butt</ogc:Literal>
</CssParameter>
<CssParameter name="stroke-linejoin">
<ogc:Literal>miter</ogc:Literal>
</CssParameter>
<CssParameter name="stroke-width">
<ogc:Literal>1</ogc:Literal>
</CssParameter>

</Stroke>
</LineSymbolizer>
</Rule>

<Rule>

<Name>rule06</Name>
<Title>100 foot contours index</Title>
<Abstract>Abstract</Abstract>

<ogc:Filter>
<ogc:PropertyIsEqualTo>
<ogc:PropertyName>div_100</ogc:PropertyName>
<ogc:Literal>0</ogc:Literal>
</ogc:PropertyIsEqualTo>
</ogc:Filter>

<MinScaleDenominator>7200</MinScaleDenominator>
<MaxScaleDenominator>20999</MaxScaleDenominator>

<LineSymbolizer>
<Stroke>
<CssParameter name="stroke">
<ogc:Literal>#D95F02</ogc:Literal>
</CssParameter>
<CssParameter name="stroke-linecap">
<ogc:Literal>butt</ogc:Literal>
</CssParameter>
<CssParameter name="stroke-linejoin">
<ogc:Literal>miter</ogc:Literal>
</CssParameter>
<CssParameter name="stroke-width">
<ogc:Literal>1</ogc:Literal>
</CssParameter>

</Stroke>
</LineSymbolizer>
</Rule>

<Rule>

<Name>rule07</Name>
<Title>50 ft contours</Title>
<Abstract>Abstract</Abstract>

<ogc:Filter>
<ogc:PropertyIsEqualTo>
<ogc:PropertyName>div_50</ogc:PropertyName>
<ogc:Literal>0</ogc:Literal>
</ogc:PropertyIsEqualTo>
</ogc:Filter>

<MinScaleDenominator>21000</MinScaleDenominator>
<MaxScaleDenominator>100000</MaxScaleDenominator>

<LineSymbolizer>
<Stroke>
<CssParameter name="stroke">
<ogc:Literal>#ffd700</ogc:Literal>
</CssParameter>
<CssParameter name="stroke-linecap">
<ogc:Literal>butt</ogc:Literal>
</CssParameter>
<CssParameter name="stroke-linejoin">
<ogc:Literal>miter</ogc:Literal>
</CssParameter>
<CssParameter name="stroke-width">
<ogc:Literal>1</ogc:Literal>
</CssParameter>

</Stroke>
</LineSymbolizer>
</Rule>

<Rule>

<Name>rule08</Name>
<Title>250 foot contours index</Title>
<Abstract>Abstract</Abstract>

<ogc:Filter>
<ogc:PropertyIsEqualTo>
<ogc:PropertyName>div_250</ogc:PropertyName>
<ogc:Literal>0</ogc:Literal>
</ogc:PropertyIsEqualTo>
</ogc:Filter>

<MinScaleDenominator>21000</MinScaleDenominator>
<MaxScaleDenominator>100000</MaxScaleDenominator>

<LineSymbolizer>
<Stroke>
<CssParameter name="stroke">
<ogc:Literal>#D95F02</ogc:Literal>
</CssParameter>
<CssParameter name="stroke-linecap">
<ogc:Literal>butt</ogc:Literal>
</CssParameter>
<CssParameter name="stroke-linejoin">
<ogc:Literal>miter</ogc:Literal>
</CssParameter>
<CssParameter name="stroke-width">
<ogc:Literal>1</ogc:Literal>
</CssParameter>

</Stroke>
</LineSymbolizer>
</Rule>

</FeatureTypeStyle>
</UserStyle>
</NamedLayer>
</StyledLayerDescriptor>

Posted in GeoServer | Tagged: , , , | 1 Comment »

Geometry Collections and Small Headaches (pics)

Posted by smathermather on February 10, 2011

Posted in Database, GeoServer, PostGIS, SQL | Tagged: , , , , | Leave a Comment »

Geometry Collections and Small Headaches

Posted by smathermather on February 10, 2011

The disk that held the personal geodatabases of our contour datasets died a while back, but not before I loaded the contours into PostGIS and started serving them up. Our new intern is working on putting together some shapefiles layer groups in ArcGIS for map production, and asked for a missing one… . I don’t have it, but in principle can always extract it from PostGIS and serve it back as a shapefile.

So I performed an intersection:


CREATE TABLE public.contour_2_cut AS
	SELECT contour.elevation, (ST_Intersection(contour.the_geom), cutline) AS the_geom, gid FROM
		public.contour_2 AS contour,  (ST_GeomFromText('POLYGON((2000000 860000, 2500000 860000, 2500000 460000, 2000000 460000,2000000 860000))', 3734)) AS cutline;

After adding a primary key:


ALTER TABLE public.contour_2_cut ADD PRIMARY KEY (gid);

I couldn’t display it in uDig, which complained that it was a geometry collection. Ah, of course. In the ESRI world (whence I hark) an intersection takes the form of the least of the geometries involved, so the intersection of a polyline and a polygon is a polyline. Not so in PostGIS. PostGIS supports geometry collections, so naturally the intersection includes both types of geometry. (Side note– the viewing in uDig didn’t fail until it got to the polygon– it first loaded all the contour polylines, then discovering the polygon, threw an error).

What to do? Well, I hoped pgsql2shp would bail me out:


pgsql2shp -f contours_block_08_01 -h cmac-srv-gis -u postgres CM public.contour_2_cut
Initializing... type 'GEOMETRYCOLLECTION' is not Supported at this time.
The DBF file will be created but not the shx or shp files.
You've found a bug! (pgsql2shp.c:2864)

But it doesn’t support GEOMETRYCOLLECTIONs either.

Time to dig deeper. So how do we tell what a geometry is?


SELECT ST_GeometryType(st_geom) FROM contour_2_cut;

and I see ST_LineString is the descriptor for our lines. So, let’s create a copy of the table with just linestrings (we could remove records from the existing table as well):


CREATE TABLE contour_2_cut_line AS
	SELECT * FROM contour_2_cut WHERE
		ST_GeometryType(the_geom) ='ST_LineString';

Posted in Database, PostGIS, SQL | Tagged: , , , , | 5 Comments »

Optimizing PostGIS Geometry Functions– Mmm. Donut.

Posted by smathermather on March 6, 2010

Well, how can I PostGIS without referring back to BostonGIS.  One day I’ll know… .  But first, the Simpsons:  “Mmm… … Donut.”

A while back I read Regina Obe’s BostonGIS post on “Map Dicing,” essentially using a grid to dice up an existing dataset to a small granularity in order to take full advantage of spatial indexing, etc..  The post really intrigued me, and is often at the back of my mind as I’m working on doing heavy duty PostGIS queries (and waiting hours/days for them to complete).

There was a time when I worked a job (not long ago) and we would start data processes days before we expected to see them complete– sometimes 6-7 days.  So, generally, I don’t mind waiting (especially if I know something will likely complete).  Imagine my surprise at my own lack of patience then when I ran a couple of analyses in PostGIS, and they were still running 2 days later.  Well, I don’t want to be impatient, but I also don’t want to be inefficient!

A few posts back, I deleted some bad data from my contour dataset by chopping it out willy-nilly with an ST_Difference command.  It ran pretty quickly.  As in, overnight, which by most measures is fast enough.

So now, I’m trying to paste together the best of what I have for contours.  For our most urban county, I have 2ft contours derived from the combo of LiDAR and 3D breaklines– a beautiful (but PitA) dataset.  For the remaining 6 counties we want to map, I have 2 ft contours derived from (what I think is) a TIN interpolation of LiDAR points, that isn’t quite so detailed, or as effective at representing roads, cliffs, and waterways, but good enough for most purposes.  Sadly, the urban county sits surrounded by the other counties, and the bounds of the two datasets do not mesh perfectly.

In other words, I’ve got one cookie cutter and with it I want to cut a hole in one data set (our future donut), and trim another dataset to the same shape (our future donut hole), and then stick the two datasets together.

So using my two favorite cookie cutters, ST_Difference and ST_Intersection, I thought I’d use the county boundary of the urban area to remove all data within that county from the 7 county dataset (ST_Difference), and “clip” the urban county to the same strict boundary, stick the two datasets together, and wallah!

Problem is, it’s not as fast this time.  There’s two reasons for this.  On the dataset I am making into a donut, it is a dataset about 3 times as big.  So, naturally it should run more slowly.  On the dataset I’m making into the donut hole, however, the problem is with the inefficiency of meaningless bounding boxes– e.g. the bounding box of the cookie cutter is nearly as big as the full extent of the contour dataset, therefore, nearly all the tests against the bounding box as a filter to see whether there’s any likelihood of an intersection return a positive result, whether the contours in question are near the edge or not.

Enter “Map Dicing”.  I need smaller chunks of data.

Since I already have a grid I like to use for the county, I can skip Regina’s first step of creating a grid, and jump right into dicing up my cookie cutter:

CREATE TABLE county_cutter
(
 gid serial NOT NULL PRIMARY KEY,
 name character varying(15)
 );

SELECT AddGeometryColumn('public', 'county_cutter', 'the_geom', 9102722, 'MULTIPOLYGON', 2);

INSERT INTO county_cutter(name,the_geom)
SELECT name,
 ST_multi(ST_Intersection(c.the_geom, cg.the_geom))
 FROM county c INNER JOIN county_grid cg ON ST_Intersects(c.the_geom, cg.the_geom);

CREATE INDEX county_cutter_the_geom_idx ON county_cutter USING gist(the_geom);

Now, I can get back to geo-processing.

CREATE TABLE public.cuy_contours_2_donut_hole
 AS
 SELECT elevation,
       ST_Intersection(a.the_geom, b.the_geom) AS the_geom
    FROM base.cuy_contours_2_1 as a,
         public.county_cutter as b;

CREATE TABLE public.contours_2_donut
 AS
 SELECT elevation,
       ST_Difference(a.the_geom, b.the_geom) AS the_geom
    FROM base.contours_2 as a,
         public.county_cutter as b;

Well, I should do this with a union statement in order to combine it into one table in one go.  I’ll revise before I run it, and report back whether it works.

Posted in Database, Database Optimization, PostGIS, SQL | Tagged: , , , , , | 1 Comment »

Contour data and table management in PostGIS

Posted by smathermather on February 23, 2010

My contour dataset is almost ready to go live.  Just a few more tweaks.  I was manipulating it in the public schema.  I really want it in my base map schema “base”.

Want to change schemas on a table in Postgre? No problem:

ALTER TABLE public.cuy_contours_2_clip SET SCHEMA base;

Rename a table? No problem:

ALTER TABLE base.cuy_contours_2_clip RENAME TO base.cuy_contours_2;

Save space by not storing higher precision data than you need to (e.g. NUMERIC data for what should be an integer elevation):

ALTER TABLE base.cuy_contours_2
 ALTER COLUMN elevation TYPE smallint;

Now I want pre-built queries that store the two foot contour data as only the 10 foot, 20 foot, 50 foot, etc..  These different versions of the dataset will then have zoom factors applied to them when I serve them up in GeoServer.  To achieve this, we use Postgre’s modulus operator with WHERE statement, e.g.:

CREATE TABLE base.cuy_contours_250
 AS
 SELECT elevation, the_geom
 FROM base.cuy_contours_2
 WHERE base.cuy_contours_2.elevation % 250 = 0;

And the rest:

CREATE TABLE base.cuy_contours_10
 AS
 SELECT elevation, the_geom
 FROM base.cuy_contours_2
 WHERE base.cuy_contours_2.elevation % 10 = 0;
CREATE TABLE base.cuy_contours_20
 AS
 SELECT elevation, the_geom
 FROM base.cuy_contours_2
 WHERE base.cuy_contours_2.elevation % 20 = 0;
CREATE TABLE base.cuy_contours_50
 AS
 SELECT elevation, the_geom
 FROM base.cuy_contours_2
 WHERE base.cuy_contours_2.elevation % 50 = 0;
CREATE TABLE base.cuy_contours_100
 AS
 SELECT elevation, the_geom
 FROM base.cuy_contours_2
 WHERE base.cuy_contours_2.elevation % 100 = 0;

Posted in Database, Database Optimization, PostGIS, SQL | Tagged: , , , , | 1 Comment »

Mapping places unknown– free global datasets and FOSS GIS are a great combo

Posted by smathermather on January 31, 2010

I wanted to put together a quick and dirty map of a biological reserve in Ecuador, sort of a laptop exploration of a place quite distant.  At first, I thought I’d use Shuttle Radar Topography Misssion data to get the elevation information.  Then I discovered the ASTER Global DEM which is 30m resolution for the whole world.  Wow.  Cool cool data. (I used the Japanese site to download, as I had trouble getting the US site to load).

So, first I downloaded and extracted the zipped geotiffs that come from the ASTER Global DEM.  I extracted the DEMs into the same directory, and mosaicked them into a single DEM:

gdalwarp AST*.tif ecuador_aster_dem.tif

Then I created a hillshade version of the DEM for mapping:

gdaldem hillshade ecuador_aster_dem.tif ecuador_aster_hillshade.tif

and then 100 and 50 meter contours for the site:

gdal_contour -a elev -3d -i 100 ecuador_aster_dem.tif ecuador_aster_contour_100.shp
gdal_contour -a elev -3d -i 50 ecuador_aster_dem.tif ecuador_aster_contour_50.shp

I loaded it all into Quantum GIS, symbolized it, and added some OpenStreetMap data.  To do this, I navigated to the location in OSM, chose the download tab, and chose OpenStreetMap XML Data as my output.  Using the OSM plug-in in QGIS, I imported the file, and walla– a nearly functional map.

Posted in GDAL | Tagged: , , , | Leave a Comment »

Clipping Data w/ PostGIS continued

Posted by smathermather on January 29, 2010

It finished, yay!  A little Friday 5PM exuberance, I’m afraid.  See previous post for explanation.

Here is the clipped version:

Posted in Database, PostGIS, SQL | Tagged: , , , , , | Leave a Comment »