tag:blogger.com,1999:blog-90415776667262269472024-03-13T14:27:38.242-06:00Moving SpatialA personal gallery of geospatial musings.Anonymoushttp://www.blogger.com/profile/16259376360116803538noreply@blogger.comBlogger24125tag:blogger.com,1999:blog-9041577666726226947.post-82170068159023406072015-12-06T08:02:00.001-07:002015-12-06T08:03:45.026-07:00Federal to State Land Transfers - Colorado MapThere have been numerous attempts across the West to transfer Federal land ownership over to the states (<a href="http://www.denverpost.com/willoughby/ci_28056540/willoughby-federal-land-transfer-idea-fails-but-wont">Denver Post</a>). Luckily, they have been defeated, however, this issue is not likely to go away.<br />
<br />
The idea behind transferring land comes from this idea that states can manage land within their borders better than the federal government. Although most would not argue that federal management is perfect, it would be foolish to think that states could and would manage them better. States have a history of selling off their public lands for economic gain. I encourage you to read <a href="http://www.trcp.org/assets/pdf/Locked-Out-Report-Sportsmens-Access.pdf">Locked Out</a>, a report put out by the <a href="http://www.trcp.org/">Theodore Roosevelt Conservation Partnership</a> (TRCP). The report describes - in plain English - the implications and consequences of any land transfer.<br />
<br />
As for me, I decided to look into State lands in Colorado to see how they are managed. In Colorado, I think most outdoors folk are aware, state land is by default <b>NOT open to the public (</b><a href="http://trustlands.state.co.us/Projects/Pages/Recreation.aspx">Colorado State Land Board</a><b>)</b>. Of the 3 million acres of state land, only 500,000 are open to the public for wildlife related activities. This public access is made possible by a public access program between Colorado Parks and Wildlife (CPW) and the Colorado State Lands Board. The rest are leased to generate revenue for Colorado public schools and other public institutions (<a href="http://cpw.state.co.us/placestogo/Pages/StateTrustLands.aspx">CPW</a>). This means that only <b>~16% of state lands are open to the public</b>.<br />
<div>
<br /></div>
<div>
So if the state of Colorado were to take over all Forest Service lands, what would public access look like on them? We can only assume that unless specifically specified otherwise, the state would treat these lands in the same way they treat current lands under their management - 16% access to the public. So what does 16% of Forest Service land in Colorado look like? In the map below, the light green depicts current Forest Service ownership. The dark green represents 16% of the land area. A stark reduction in public access.<br />
<br />
Please, write your representatives and let them know that you are against federal to state land transfers. The TRCP has put together <a href="https://secure3.convio.net/trcp/site/Advocacy?cmd=display&page=UserAction&id=415">a petition</a> ready for you to sign.</div>
<div>
<br /></div>
<div>
<br /></div>
<iframe allowfullscreen="" frameborder="0" height="520" mozallowfullscreen="" msallowfullscreen="" oallowfullscreen="" src="https://alopez04.cartodb.com/viz/2935ea64-9b75-11e5-88a5-0e787de82d45/embed_map" webkitallowfullscreen="" width="100%"></iframe>Anonymoushttp://www.blogger.com/profile/16259376360116803538noreply@blogger.com0tag:blogger.com,1999:blog-9041577666726226947.post-31474928942812190242015-11-30T14:25:00.000-07:002015-12-01T10:06:15.823-07:00H.R. 1931 American Land Act - Map SimulationIn April of 2015, the American Land Act bill was introduced by Ted Poe (R-TX). As of the date on this blog post, the bill is still active (<a href="https://www.congress.gov/bill/114th-congress/house-bill/1931">congress.gov</a>). <b>The bill would force the sell of 8% of National Forest Service land and 8% of BLM land for each year of 2016-2021</b>. This amounts to about 1/3 of each BLM and Forest Service land being sold with proceeds going to federal transportation projects.<br />
<br />
Being an avid outdoorsman, I was alarmed at the sheer magnitude of land to be sold. And to be honest, appalled such a bill would even be introduced. In any regard, I am not as articulate as others who might make an impassioned argument against the sell of our Public lands as well as I am sure there are other folks out there that can argue the finer points of such legislation. These are not my intentions with this post. My intention is to allow the American public a chance to conceptualize the magnitude of land this bill would sell off.<br />
<br />
Below you'll see a map of Forest Service lands (green) and Privatized lands (dark blue). This represents a simulation of lands to be sold off by 2021 if the American Land Act were to come to fruition.<br />
<br />
The simulation iteratively reduced Forest Service lands by 8% for each year of 2016-2021. Each year, the total amount of land left was recalculated then reduced. <b>The parcels removed were chosen at random</b>. I did not make an attempt to eliminate lands from this simulation I thought might be off limits as I didn't want to make that judgement call. <br />
<br />
<br />
<iframe allowfullscreen="" frameborder="0" height="520" mozallowfullscreen="" msallowfullscreen="" oallowfullscreen="" src="https://alopez04.cartodb.com/viz/a2c1fbf8-9796-11e5-92a4-0e787de82d45/embed_map" webkitallowfullscreen="" width="100%"></iframe>
<br />
<br />
The random draw has the effect of spatially distributing the hypothetical privatized lands. Perhaps a more powerful example is to sum the area required to meet the bills 1/3 reduction state-by-state. In such case, it would roughly require <b>all</b> the Forest Service land in Colorado, Utah, Nevada, Wyoming, Montana, and Idaho.<br />
<br />
<iframe allowfullscreen="" frameborder="0" height="520" mozallowfullscreen="" msallowfullscreen="" oallowfullscreen="" src="https://alopez04.cartodb.com/viz/0a192686-97c5-11e5-8cd5-0e674067d321/embed_map" webkitallowfullscreen="" width="100%"></iframe>
<br />
Keep in mind this is only looking at Forest Service lands and does not include the 1/3 of BLM lands to be sold off. Luckily, the bill <a href="http://poe.house.gov/2015/4/congressman-poe-introduces-the-american-land-act">does not authorize</a> the sell of National Parks or U.S. Fish and Wildlife lands....<br />
<br />
<br />
I urge anyone who loves the outdoors and our public lands to share this post, write your congressman/woman and make your voice heard. Our public lands should not be for sell!Anonymoushttp://www.blogger.com/profile/16259376360116803538noreply@blogger.com0tag:blogger.com,1999:blog-9041577666726226947.post-58344734741069124302013-06-14T12:36:00.000-06:002013-06-14T12:38:15.551-06:00PostgreSQL - Copy Table to Another Server AND Different SchemaHere I present a method for dumping a table from one server/database and loading into another server/database with a different schema.<br />
<br />
This is a spinoff of a previous post (<a href="http://movingspatial.blogspot.com/2012/11/postgresql-copy-table-to-another-server.html" target="_blank">PostgreSQL - Copy Table to Another Server</a>) which shows you how to dump a table from one server to another. The major downfall of that was the inability to specify a different schema - the schema of the table being dumped had to exist on the new server and database.<br />
<br />
<ul class="ul1">
<li class="li1"><i>schema.table</i>: this is your table to be transferred</li>
<li class="li1">-O: do not save ownership</li>
<li class="li1"><i>from_db</i>: from database</li>
<li class="li1"><span class="s1"><b>from_schema</b></span>: from schema</li>
<li class="li1"><span class="s2"><b>to_schema</b></span>: to schema</li>
<li class="li1"><i>to_db</i>: to database</li>
</ul>
<br />
<div class="p1">
<b>pg_dump</b> -h serverA -t <i>schema.table</i> -O <i>from_db</i> | <b>sed</b> -e '/^SET search_path = /s/<span class="s1"><b>from_schema</b></span>/<b>to_schema</b>/g' | <b>psql</b> -h serverB<i> to_db</i></div>
<div class="p1">
<i><br /></i></div>
<div class="p1">
<i><br /></i></div>
Anonymoushttp://www.blogger.com/profile/16259376360116803538noreply@blogger.com3tag:blogger.com,1999:blog-9041577666726226947.post-45795848119960594282013-01-31T20:21:00.000-07:002013-01-31T20:48:02.438-07:00PostGIS 1.5 - Batch Populate Geometry Columns<br />
<div class="p1">
When creating new spatial tables – say as the result of an analysis – you should register and add constraints to your table. At least this is the proper thing to do… If you are like me, I get a little lazy and the next thing I know I have 50 tables without proper registration. This is fine if you do not plan on using the tables again or simply plan to dump out as a shapefile, but if you plan on keeping them around and using again, you should register and add constraints.</div>
<div class="p1">
<br /></div>
<div class="p1">
You can do this on one table simply by running: </div>
<div class="p1">
<code style="font-size: 12px;"><span style="color: blue;">SELECT </span>Populate_Geometry_Columns<span style="color: grey;">(</span><span style="color: red;">'schema.table'</span><span style="color: grey;">::</span>regclass<span style="color: grey;">)<span style="font-family: arial, verdana, sans-serif;">;</span></span></code></div>
<div class="p1">
<br /></div>
<div class="p1">
However, if you have many tables, you can run the following on an entire schema:</div>
<div class="p1">
<code style="font-size: 12px;"><span style="color: blue;">WITH</span></code><br />
<code style="font-size: 12px;">t <span style="color: blue;">AS </span><span style="color: grey;">(</span><span style="color: blue;">SELECT </span>table_schema<span style="color: grey;">, </span>table_name<br /> <span style="color: blue;">FROM </span>information_schema.tables<br /> <span style="color: blue;">WHERE </span>table_schema <span style="color: blue;">= </span><span style="color: red;">'foo' </span><span style="color: grey;">AND </span>table_type <span style="color: blue;">= </span><span style="color: red;">'BASE TABLE'</span><span style="color: grey;">)</span></code><br />
<code style="font-size: 12px;"><span style="color: blue;">SELECT </span>Populate_Geometry_Columns<span style="color: grey;">((</span>table_schema::text <span style="color: grey;">|| </span><span style="color: red;">'.' </span><span style="color: grey;">|| </span>table_name::text<span style="color: grey;">)::</span>regclass<span style="color: grey;">)</span><span style="color: blue;">FROM </span>t<span style="font-family: arial, verdana, sans-serif;">;</span></code></div>
Anonymoushttp://www.blogger.com/profile/16259376360116803538noreply@blogger.com0tag:blogger.com,1999:blog-9041577666726226947.post-88284385784550071952013-01-29T11:28:00.001-07:002013-01-29T11:30:02.559-07:00PostgreSQL - Determining Table and Index SizeThis is somewhat of a follow up to my previous post. If you are interested in knowing the size of your postgres table and associated indexes, sequences, etc., you can run the following SQL to find out. Enjoy.<br />
<br />
<code style="font-size: 12px;"><span style="color: blue;">SELECT </span><span style="color: black;">pg_size_pretty</span><span style="color: grey;">(</span><span style="color: black;">pg_relation_size</span><span style="color: grey;">(</span><span style="color: black;">pg_class.oid</span><span style="color: grey;">)), </span><span style="color: black;">pg_class.relname</span><span style="color: grey;">, </span><span style="color: black;">pg_namespace.nspname</span><span style="color: grey;">,<br /> </span><span style="color: magenta;">CASE </span><span style="color: black;">pg_class.relkind </span><span style="color: blue;">WHEN </span><span style="color: red;">'r' </span><span style="color: blue;">THEN </span><span style="color: red;">'table' </span><span style="color: blue;">WHEN </span><span style="color: red;">'i' </span><span style="color: blue;">THEN </span><span style="color: red;">'index' </span><span style="color: blue;">WHEN </span><span style="color: red;">'S' </span><span style="color: blue;">THEN </span><span style="color: red;">'sequence' </span><span style="color: blue;">WHEN </span><span style="color: red;">'v' </span><span style="color: blue;">THEN </span><span style="color: red;">'view' </span><span style="color: blue;">WHEN </span><span style="color: red;">'t' </span><span style="color: blue;">THEN </span><span style="color: red;">'TOAST' <br /> </span><span style="color: blue;">ELSE </span><span style="color: black;">pg_class.relkind::text </span><span style="color: blue;">END<br /> FROM </span><span style="color: black;">pg_class<br /> </span><span style="color: magenta;">LEFT </span><span style="color: grey;">OUTER </span><span style="color: blue;">JOIN </span><span style="color: black;">pg_namespace </span><span style="color: blue;">ON </span><span style="color: grey;">(</span><span style="color: black;">pg_namespace.oid </span><span style="color: blue;">= </span><span style="color: black;">pg_class.relnamespace</span><span style="color: grey;">)<br /> </span><span style="color: blue;">ORDER BY </span><span style="color: black;">pg_relation_size</span><span style="color: grey;">(</span><span style="color: black;">pg_class.oid</span><span style="color: grey;">) </span><span style="color: blue;">DESC</span><span style="color: grey;">;</span></code>Anonymoushttp://www.blogger.com/profile/16259376360116803538noreply@blogger.com0tag:blogger.com,1999:blog-9041577666726226947.post-77938685328202920472012-11-21T12:55:00.002-07:002012-11-21T18:08:02.121-07:00PostGIS - Raster Tabulate Area This post compares - in a very basic sense - tabulating raster area in ESRI and PostGIS.<br />
<br />
<b><span style="color: #073763;">1. Data</span></b><br />
- Raster: 90 meter resolution covering most of the contiguous United States. Represents available land for solar development.<br />
- Vector: Native American lands from the Bureau of Indian Affairs. Roughly 300 distinct boundaries.<br />
<br />
<b><span style="color: #073763;">2. Storage</span></b><br />
-ESRI: 113 MB<br />
-PostGIS: Loaded in using 100x100 tiles. Table = 137 MB, spatial index = 19 MB<br />
<br />
<span style="color: #073763;"><b>3. Tabulate Area - Time Trial</b></span><br />
-ESRI: tabulate area tool finished in roughly 1 minute.<br />
-PostGIS: Using the SQL below, finished in roughly 13 minutes. Confirmed spatial index was used.<br />
<br />
<code style="font-size: 12px;"><span style="color: blue;">CREATE TABLE </span><span style="color: black;">tabulate_area </span><span style="color: blue;">AS<br />SELECT </span><span style="color: black;">a.tribe_name</span><span style="color: grey;"> </span></code><br />
<code style="font-size: 12px;"><span style="color: grey;"> ,</span><span style="color: magenta;">SUM</span><span style="color: grey;">(</span><span style="color: black;">ST_Count</span><span style="color: grey;">(</span><span style="color: black;">ST_Clip</span><span style="color: grey;">(</span><span style="color: black;">b.rast</span><span style="color: grey;">, </span><span style="color: black;">1</span><span style="color: grey;">, </span><span style="color: black;">a.the_geom</span><span style="color: grey;">)) </span></code><span style="color: grey; font-family: monospace; font-size: 12px;">* </span><span style="color: black; font-family: monospace; font-size: 12px;">pow</span><span style="color: grey; font-family: monospace; font-size: 12px;">(</span><span style="color: black; font-family: monospace; font-size: 12px;">90</span><span style="color: grey; font-family: monospace; font-size: 12px;">,</span><span style="color: black; font-family: monospace; font-size: 12px;">2)</span><span style="color: grey; font-size: 12px;">)</span><span style="color: grey; font-size: 12px;"> </span><span style="color: blue; font-size: 12px;">AS </span><span style="color: black; font-size: 12px;">square_meters</span><br />
<code style="font-size: 12px;"><span style="color: blue;">FROM </span><span style="color: black;">tribal_boundaries a</span></code><br />
<code style="font-size: 12px;"><span style="color: blue;">INNER JOIN </span><span style="color: black;">available_land b </span><span style="color: blue;">ON </span><span style="color: black;">ST_Intersects</span><span style="color: grey;">(</span><span style="color: black;">a.the_geom</span><span style="color: grey;">, </span><span style="color: black;">b.rast</span><span style="color: grey;">)</span></code>
<br />
<code style="font-size: 12px;"><span style="color: grey;"><br /></span></code>
<b><span style="color: #073763;">4. Compare Results</span></b><br />
I found that PostGIS likes to automatically snap to the raster whereas ESRI defaults to not snapping. If not snapping in ESRI, results between ESRI and PostGIS are different (as would be expected). If snapping in ESRI, then results are exactly the same. Loaded in the ESRI results and compared the two with the following SQL:<br />
<br />
<code style="font-size: 12px;"><span style="color: blue;">SELECT<br /> </span><span style="color: black;">post.tribe_name<br /> </span><span style="color: grey;">,</span><span style="color: black;">post.sqm </span><span style="color: blue;">AS </span><span style="color: black;">post_sqm<br /> </span><span style="color: grey;">,</span><span style="color: black;">esri.sqm </span><span style="color: blue;">AS </span><span style="color: black;">esri_sqm<br /> </span><span style="color: grey;">,</span><span style="color: black;">post.sqm </span><span style="color: grey;">- </span><span style="color: black;">esri.sqm </span><span style="color: blue;">AS </span><span style="color: black;">tabulate_diff</span></code><br />
<code style="font-size: 12px;"><span style="color: blue;">FROM </span><span style="color: black;">tabulate_area </span><span style="color: blue;">AS </span><span style="color: black;">post</span></code><br />
<code style="font-size: 12px;"><span style="color: magenta;">LEFT </span><span style="color: grey;">OUTER </span><span style="color: blue;">JOIN </span><span style="color: black;">esri </span><span style="color: blue;">ON </span><span style="color: black;">post.tribe_name </span><span style="color: blue;">= </span><span style="color: black;">esri.tribe_name</span><span style="color: grey;">;</span></code>
<br />
<br />
<b><span style="color: #073763;">5. Why is PostGIS so slow? It's not! Read on...</span></b><br />
Well...12 minutes slower is simply not acceptable...so what is the reason? I decided to look at the tile size on the raster. My initial thought was, the tighter the tiles, the better the spatial index, so I set it at 100x100 or roughly 9 kilometers ((90cell sizex100tile size)/1000meters). Thinking on it, sure my ST_Intersects() was fast, but my clip was choking on all the individual tiles. So I loaded in the raster again, this time at 1000x1000tile size. Ran the query again and que bueno! <b>New runtime for PostGIS is 1 1/2 minutes</b>.<br />
<br />
Lesson learned: think about tile size and the analysis you'll be performing before you load your raster.<br />
<br />Anonymoushttp://www.blogger.com/profile/16259376360116803538noreply@blogger.com0tag:blogger.com,1999:blog-9041577666726226947.post-37040392716965823412012-11-01T22:24:00.000-06:002012-11-02T09:03:57.647-06:00PostGIS - Majority Point Value in Polygon<div class="separator" style="clear: both; text-align: left;">
Here I present a simple method for determining the majority point value within each polygon in PostGIS. This method utilizes the postgres Window function to rank order the count of point values within each polygon. </div>
<div class="separator" style="clear: both; text-align: center;">
<br /></div>
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjjPEEFdhMH4cOjQxr50xktmP-uflSIobURVp0MyHT6q0MYc57T-OVOr8mBI9bPEJKQFTNOWTSXTZTzthNFdJ6tAvvlN79kt7HW5GrEs-4ToybgPsH0TBwZ_ZdXfq0lY6bd0JyoMgKwLEH0/s1600/Untitled.jpg" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" height="267" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjjPEEFdhMH4cOjQxr50xktmP-uflSIobURVp0MyHT6q0MYc57T-OVOr8mBI9bPEJKQFTNOWTSXTZTzthNFdJ6tAvvlN79kt7HW5GrEs-4ToybgPsH0TBwZ_ZdXfq0lY6bd0JyoMgKwLEH0/s320/Untitled.jpg" width="320" /></a></div>
<br />
<br />
<br />
<code style="font-size: 12px;"><span style="color: green;">--Majority Count of Point Values</span></code><br />
<code style="font-size: 12px;"><span style="color: blue;">SELECT<br /> </span><span style="color: black;">i.</span><span style="color: grey;">*</span></code><br />
<code style="font-size: 12px;"><span style="color: blue;">FROM</span><span style="color: grey;">(</span></code><br />
<code style="font-size: 12px;"><span style="color: blue;">SELECT<br /> </span><span style="color: black;">polygons.gid<br /> </span><span style="color: grey;">,</span><span style="color: black;">points.value<br /> </span><span style="color: grey;">,</span><span style="color: magenta;">COUNT</span><span style="color: grey;">(*) </span><span style="color: blue;">AS </span><span style="color: black;">point_count<br /> </span><span style="color: grey;">,</span><span style="color: black;">RANK</span><span style="color: grey;">() </span><span style="color: blue;">OVER</span><span style="color: grey;">(</span><span style="color: black;">PARTITION </span><span style="color: blue;">BY </span><span style="color: black;">polygons.gid </span><span style="color: blue;">ORDER BY </span><span style="color: magenta;">COUNT</span><span style="color: grey;">(*) </span><span style="color: blue;">DESC</span><span style="color: grey;">) </span><span style="color: blue;">AS </span><span style="color: black;">point_value_rank</span></code><br />
<code style="font-size: 12px;"><span style="color: blue;">FROM </span><span style="color: black;">polygons</span><span style="color: grey;">, </span><span style="color: black;">points</span></code><br />
<code style="font-size: 12px;"><span style="color: blue;">WHERE </span><span style="color: black;">ST_Intersects</span><span style="color: grey;">(</span><span style="color: black;">polygons.geom</span><span style="color: grey;">, </span><span style="color: black;">points.geom</span><span style="color: grey;">)</span></code><br />
<code style="font-size: 12px;"><span style="color: blue;">GROUP BY </span><span style="color: black;">polygons.gid</span><span style="color: grey;">, </span><span style="color: black;">points.value </span><span style="color: grey;">) </span><span style="color: black;">i</span></code><br />
<code style="font-size: 12px;"><span style="color: blue;">WHERE </span><span style="color: black;">i.point_value_rank </span><span style="color: blue;">= </span><span style="color: black;">1</span><span style="color: grey;">; </span></code><br />
<br />
Oh, and by the way. My SQL has always look terrible because I usually get lazy halfway through edited the color etc (apologies). But now I have come across <a href="http://extras.sqlservercentral.com/prettifier/prettifier.aspx" target="_blank">SQL "Prettify"</a>. Neat.Anonymoushttp://www.blogger.com/profile/16259376360116803538noreply@blogger.com0tag:blogger.com,1999:blog-9041577666726226947.post-88744523390078786832012-11-01T21:55:00.001-06:002012-11-01T21:56:08.205-06:00PostgreSQL - Copy Table to Another Server<br />
<div class="p1">
My previous workflow for copying tables from one server to another was (1) pg_dump on server "A" then (2) SCP to server "B" then (3) restore via psql. Not a terrible workflow, but my new workflow is much more eloquent. </div>
<div class="p1">
<br /></div>
<div class="p1">
First, I login to server B (my new home for the data), then pg_dump with connection params set to serverA (current location of data). Then pipe directly to psql with connection params set to serverB a.k.a. localhost. </div>
<div class="p1">
<br /></div>
<div class="p1">
<b>pg_dump</b> -t <i>schema.table</i> -h serverA -U <i>username</i> databaseA | <b>psql</b> -h localhost -U <i>username</i> databaseB;</div>
<div class="p1">
<br /></div>
<div class="p1">
An important side note:</div>
<div class="p1">
<br /></div>
<div class="p1">
If your postgres databases are of different vintages, you should always pg_dump using the latest vintage i.e. If serverB is 9.1 and serverA is 8.4, dump using serverB. pg_* are usually backwards compatible, but not always forwards compatible. </div>
Anonymoushttp://www.blogger.com/profile/16259376360116803538noreply@blogger.com1tag:blogger.com,1999:blog-9041577666726226947.post-19905635757384810842012-10-27T14:30:00.000-06:002012-10-28T08:43:22.179-06:00Split CSV into multiple filesReceiving files in a format that is...less than helpful for your project needs is almost an everyday occurrence. On this occurrence I received roughly 30,000 files of hourly solar irradiance data. The files had multiple years all crammed together i.e. file 1 had hourly data from years 2002 - 2008. In my project I needed to separate the individual years of hourly data - Bash to the rescue.<br />
<br />
#Get a list of all CSV files<br />
FILES=$(find -name *.csv)<br />
#loop through the files<br />
<b>for</b> file <b>in</b> $FILES<br />
<b>do</b><br />
#Get the basename of the current file<br />
base=${file##*/}<br />
#get the path of the current file<br />
path=${file%/*}<br />
<br />
#Split file into yearly data files.<br />
<b>split</b> -l $file 8760<br />
<br />
#Rename files to appropriate years. Split automatically names your output xaa, xab..and so<br />
# forth depending on how many files result from your split. If you know how many files<br />
#will result from your split, its easy to rename appropriately.<br />
<b>mv</b> xaa $path"/"$base"_2002.csv"<br />
<b>mv</b> xab $path"/"$base"_2003.csv"<br />
<b>mv</b> xac $path"/"$base"_2004.csv"<br />
<b>mv</b> xad $path"/"$base"_2005.csv"<br />
<b>mv</b> xae $path"/"$base"_2006.csv"<br />
<b>mv</b> xaf $path"/"$base"_2007.csv"<br />
<b>mv</b> xag $path"/"$base"_2008.csv"<br />
<b>echo</b> wrote $file<br />
<b>done</b><br />
<br />
<br />
<br />Anonymoushttp://www.blogger.com/profile/16259376360116803538noreply@blogger.com0tag:blogger.com,1999:blog-9041577666726226947.post-25765581673381437762012-10-27T13:53:00.000-06:002012-10-27T13:54:33.148-06:00Add columns of data to an existing CSV - a case for Bash over PythonI do a lot of scripting in Python, I love it, but sometimes things are just <i>more</i> easy in Bash. Take for instance appending new columns of data into a CSV file. To go about this in Python you would have to do something like the following: <a href="http://stackoverflow.com/questions/3376236/how-do-i-add-a-new-column-of-data-to-a-csv-file">Stackoverflow Example</a>. Back yet? Ok, and you would still have to add some logic to open the other file and read through its lines, etc...BUT in Bash, this command is so, so simple:<br />
<br />
#paste is the utility and "-d ," specifies the delimiter type <br />
<b>paste</b> -d , file_1.csv file_2.csv > combined_file.csv<br />
<br />
<br />Anonymoushttp://www.blogger.com/profile/16259376360116803538noreply@blogger.com0tag:blogger.com,1999:blog-9041577666726226947.post-6785699119214774682012-07-19T20:16:00.002-06:002012-07-19T20:17:06.453-06:00U.S. Renewable Energy Technical Potentials ReportFinally, finally, finally published! I've been working on this for some time and now, its finished. If your interested in using GIS to estimate renewable energy electric power potential, this report is for you.<br />
<br />
See it on the NREL GIS webpage: <a href="http://www.nrel.gov/gis/re_potential.html">http://www.nrel.gov/gis/re_potential.html</a>Anonymoushttp://www.blogger.com/profile/16259376360116803538noreply@blogger.com0tag:blogger.com,1999:blog-9041577666726226947.post-20892043789184211052012-07-15T19:35:00.001-06:002012-07-21T09:15:54.022-06:00PostGIS 2.0 - Intersect Raster and PolygonI've just taken up the task to experiment with PostGIS 2.0 raster support. At this point in time, it seems that there is little documentation (by other folks) regarding raster analysis using PostGIS, so I figured I'd write down some of what I'm experiencing. First, raster support in PostGIS is HUGE! I'm eternally thankful that ESRI is will no longer be my sole raster analysis tool (I know, there's more out there than just ESRI...).<br />
<br />
<b>The Data:</b><br />
-One 10x10km grid cell from the NREL/Perez Direct Normal Irradiance grid<br />
-90m Landscan Raster detailing population for greater Denver area from Oak Ridge National Lab<br />
<br />
<b>The Analysis Intent:</b><br />
-Summarize population within grid cell, first using ESRI and then PostGIS and compare. Easy, no? <br />
<br />
<b style="color: #073763;">1. Load Some Data</b><br />
First, loading the raster into PostGIS. This is pretty straight forward with one exception, tile size. What tile size do you use? I would love to see theoretical or empirical analysis detailing the optimal tile size per given cell size and extent. Anyone? Until then, I'll guess.<br />
<br />
-C : applies contraints, -M : vacuum analyze, -Y : use copy rather than insert, -t : tile size, -s : SRID, -U : username, -d : database<br />
<div style="color: #274e13;">
raster2pgsql -s 4326 -C -M -Y landscanden -t 100x100 ls_den | psql -U postgres -d postgis2</div>
<br />
<b style="color: #073763;">2. Get a Base Estimate </b><br />
Get a base estimate using ESRI. Using Zonal Statistics I get the following:<br />
<div style="text-align: center;">
<b>Population: 207,578</b></div>
<div style="text-align: center;">
<b>Cell Count: 14,400</b></div>
<br />
<b style="color: #073763;">3. How Much Can We Explain?</b><br />
If my results are different between ESRI and PostGIS, how much could be explained? Assume one program summarizes all border cells while the other does not. This will give us a percent difference we can easily explain. So, if I have a 10x10km grid cell as the zone, then I could fit 111 90m cells on each border for a total of 444 border cells. Square 111.9 to get 12,321 90m cells that could fit inside the 10x10km grid cell. Divide the border cells by the total 444/12,321*100 = 3.6%. We can explain<b> 3.6%</b> of the difference (if there is any) by just the border cells.<br />
<br />
<div style="color: #073763;">
<b>4. First Attempts at Raster Analysis in PostGIS</b></div>
I thought I was thinking logically when I wrote this...should have just looked at the documentation first!<br />
<div>
<span style="color: #274e13;">--Create a vectorized result of the geometric intersection between my raster layer and vector polygon </span></div>
<div>
<span style="color: #3d85c6;">CREATE TABLE</span> sum_pop <span style="color: #3d85c6;">AS</span></div>
<div>
<span style="color: #3d85c6;">SELECT</span><span style="color: #3d85c6; white-space: pre-wrap;"> </span>gid</div>
<div>
<span style="white-space: pre-wrap;"> </span>,(gv).val AS population</div>
<div>
<span style="white-space: pre-wrap;"> </span>,(gv).geom AS geom</div>
<div>
<span style="color: #3d85c6;">FROM </span>(</div>
<div>
<span style="color: #3d85c6;">SELECT </span>a.gid<span style="white-space: pre-wrap;"></span><span style="white-space: pre-wrap;"></span></div>
<div>
<span style="white-space: pre-wrap;"> </span>,ST_Intersection(b.rast, a.geom) gv<span style="white-space: pre-wrap;"> </span></div>
<div>
<span style="color: #3d85c6;">FROM </span><span style="white-space: pre-wrap;"></span>perez_grid a, ls_den b</div>
<div>
<span style="color: #3d85c6;">WHERE </span><span style="color: #3d85c6; white-space: pre-wrap;"> </span>ST_Intersects(b.rast, a.geom) ) foo;</div>
<br />
Spatial results looked good. Just needed to summarize<br />
<span style="color: #3d85c6;">SELECT </span>SUM(population), COUNT(*) <span style="color: #3d85c6;">AS </span>cell_count <span style="color: #3d85c6;">FROM </span>sum_pop;<br />
<br />
Results:<br />
<div style="text-align: center;">
<b>Population: 193,492</b></div>
<div style="text-align: center;">
<b>Cell Count: 9,811</b></div>
<div style="text-align: center;">
<b>Population Percent Difference: ~7% </b></div>
<br />
Well...that doesn't look right...how do I explain the additional 3.4%? First things first, what is the deal with the count? Well, that one is pretty easy, when the raster is vectorized, it groups cells of like values. So count doesn't really matter. <br />
<br />
<div style="color: #073763;">
<b>5. A Hunch or Swing in the Dark?</b></div>
Ok, not sure why I decided to try this, but I loaded in the exact same raster dataset except for this time I used 200x200 tile size. Shouldn't make a difference, right? WRONG!!! This still vexes me, if anyone can explain I would be grateful. <br />
<br />
Results:<br />
<div style="text-align: center;">
<b>Population: 193,446</b></div>
<div style="text-align: center;">
<b>Count: 9,787 </b></div>
<div style="color: #073763;">
<b><br /></b></div>
<b style="color: #073763;">6. The Correct Method</b><br />
After reading the documentation, I have a plan of attack. I'll follow their example and see what I get.<br />
<div>
</div>
<div>
<span style="color: #3d85c6;">CREATE TABLE</span> sum_pop2 <span style="color: #3d85c6;">AS </span></div>
<div style="color: #3d85c6;">
WITH </div>
<div>
feat <span style="color: #3d85c6;">AS </span>(<span style="color: #3d85c6;">SELECT </span>gid, geom <span style="color: #3d85c6;">FROM </span>perez_grid <span style="color: #3d85c6;">AS </span>b ),</div>
<div>
b_stats AS</div>
<div>
<span style="white-space: pre-wrap;"> </span>(<span style="color: #3d85c6;">SELECT </span>gid, (stats).*</div>
<div>
<span style="white-space: pre-wrap;"> </span> <span style="color: #3d85c6;">FROM </span>(</div>
<div>
<span style="white-space: pre-wrap;"> </span><span style="color: #3d85c6;">SELECT </span>gid, ST_SummaryStats(ST_Clip(rast,<wbr></wbr>1,geom)) <span style="color: #3d85c6;">AS</span> stats</div>
<div>
<span style="white-space: pre-wrap;"> </span><span style="color: #3d85c6;">FROM </span>ls_den</div>
<div>
<span style="white-space: pre-wrap;"> </span><span style="color: #3d85c6;">INNER JOIN</span> feat</div>
<div>
<span style="white-space: pre-wrap;"> </span><span style="color: #3d85c6;">ON </span>ST_Intersects(feat.geom,rast) ) <span style="color: #3d85c6;">AS</span> foo )</div>
<div>
<span style="color: #3d85c6;">SELECT </span>gid, SUM(count) <span style="color: #3d85c6;">AS</span> cell_count</div>
<div>
,SUM(sum) AS population</div>
<div>
<span style="white-space: pre-wrap;"> </span><span style="color: #3d85c6;">FROM </span>b_stats</div>
<div>
<span style="color: #3d85c6;">WHERE </span>count > 0</div>
<div>
<span style="white-space: pre-wrap;"> </span><span style="color: #3d85c6;">GROUP BY</span> gid</div>
<div>
<span style="white-space: pre-wrap;"> </span><span style="color: #3d85c6;">ORDER BY</span> gid;</div>
<br />
Results:<br />
<div style="text-align: center;">
<b>Population: 207,578</b></div>
<div style="text-align: center;">
<b>Cell Count: 14,400 </b></div>
<div style="text-align: center;">
<b>Population Percent Difference: 0%</b></div>
<br />
SUCCESS, however, I'm not sure why examples must be so verbose! It is a beautiful query don't get me wrong...<br />
<div style="color: #073763;">
<b><br /></b></div>
<b><span style="color: #073763;">7. Discussion</span></b><br />
So whats difference between this amd my first attempt? <br />
The WITH clause simply allows you to join queries, you could right this as a query with a subquery and get the same results. Using ST_Clip(), the raster is never converted to vector, the vector is used to clip the raster than a simple SummaryStats() is performed. That's the bulk of it. The GROUP BY clause aggregates all the rids (aka tiles).<br />
The BIG question is, why does the first attempt produce the wrong answer? The intersection produces geomval pairs for only those cells that intersect, so why the difference? I would have guessed that it would have resulted in more population because it grabs even the slightest sliver that intersects (its vector on vector) while ESRI grabs only cells whose centroid falls within the intersect, but PostGIS resulted in less! <br />
<br />
<b><span style="color: #073763;">8. A Simple Alternative to the Correct Method</span></b><br />
For those that find the above query verbose and hard to understand, I'll provide an alternative query that produces the same result in similar return time.<br />
<br />
<span style="color: #3d85c6;">CREATE TABLE</span> sum_pop3 <span style="color: #3d85c6;">AS</span><span style="color: #3d85c6;"> </span><br />
<span style="color: #3d85c6;">SELECT </span>gid, SUM((ST_SummaryStats(ST_Clip(rast,1,geom))).sum)<br />
<div>
<span style="color: #3d85c6;">FROM </span>perez_grid, ls_den </div>
<div>
<span style="color: #3d85c6;">WHERE </span>ST_Intersects(geom,rast) </div>
<div>
<span style="color: #3d85c6;">GROUP BY</span> gid</div>Anonymoushttp://www.blogger.com/profile/16259376360116803538noreply@blogger.com7tag:blogger.com,1999:blog-9041577666726226947.post-77920555665552605862012-06-16T12:21:00.000-06:002012-07-21T12:31:25.584-06:00PostgreSQL TidbitsCouple things here.<br />
<br />
<br />
<div class="p1">
1. Counting records in a postgres table is slow!</div>
<div class="p1">
<b><span style="color: #0b5394;">SELECT</span> count(*) <span style="color: #0b5394;">FROM</span> foo;</b> is slow as molasses…a good alternative is:<br />
<b><span style="color: #0b5394;">SELECT</span> reltuples <span style="color: #0b5394;">FROM</span> pg_class <span style="color: #0b5394;">WHERE</span> relname = 'foo';</b></div>
<div class="p1">
Be careful through, this method relies on the statistics built when running ANALYZE. So if your table is static, good, else make sure you keep up with analyze if using this method.</div>
<div class="p2">
<br /></div>
<div class="p1">
2. Adding a serial column in a postgres table is slow!</div>
<div class="p1">
If you need a unique id, adding a serial column will do the job and will update if you insert new records. But this is an extremely slow process. If you need a unique id and know your table will be static, you can use:<br />
<b>ROW_NUMBER() <span style="color: #0b5394;">OVER</span>() AS unique_id </b>when generating your table, and it is relatively fast. <br />
This method can also be used when generating a view to be visualized in QGIS. If you have ever tried to add a view without a unique id in QGIS, you'll get an error. To get around this, use the method above to create the unique id on the fly.</div>Anonymoushttp://www.blogger.com/profile/16259376360116803538noreply@blogger.com0tag:blogger.com,1999:blog-9041577666726226947.post-84688349912047308242012-02-12T18:02:00.000-07:002012-06-16T12:22:31.351-06:00PostgreSQL - Generate SeriesGenerate Series (<a href="http://www.postgresql.org/docs/9.1/static/functions-srf.html">PostgreSQL Doc</a>) is a powerful uncomplicated function to have in your tool-belt. Straight from the PostgreSQL help docs, generate series, generates a series from a start and stop value.<br />
<div>
<br /></div>
<div>
i.e.</div>
<div>
<span style="color: blue;">SELECT</span> * <span style="color: blue;">FROM</span> generate_series(1,3) <span style="color: blue;">AS</span> series</div>
<div>
<span style="color: #38761d;">series</span></div>
<div>
<span style="color: #38761d;">_ _ _ _ _ </span></div>
<div>
<span style="color: #38761d;">1</span></div>
<div>
<span style="color: #38761d;">2</span></div>
<div>
<span style="color: #38761d;">3</span><br />
<span style="color: #38761d;"><br /></span><br />
So, how you ask can you use this? The other week I was processing weather files in an hourly PV simulation model via Python and passing results to a PostgreSQL table. I had ~2.4million weather files (100,000 unique spatial locations * 8 years * 3 different system types) to process but was sure some would fail, but I did not want to have to restart the process every time a file failed so I threw everything in a try, except block. With the processing finished I needed to check what finished and what failed, so I needed a table with all possible outcomes to join my results to - in came generate_series. <br />
<br />
<span style="color: #38761d;">--Create table of all possible outcomes</span><br />
<span style="color: blue;">CREATE TABLE</span> unique_outcomes <span style="color: blue;">AS</span><br />
<span style="color: blue;">SELECT</span><br />
gid <span style="color: blue;">AS</span> spatial_location_id,<br />
generate_series(1998,2005) <span style="color: blue;">AS</span> years,<br />
generate_series(0,2) <span style="color: blue;">AS</span> system<br />
<span style="color: blue;">FROM</span><br />
weather_files<br />
<br />
<span style="color: #38761d;">spatial_location_id, years, system</span><br />
<span style="color: #38761d;">_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _</span><br />
<span style="color: #38761d;">1, 1998, 0</span><br />
<span style="color: #38761d;">1, 1998, 1</span><br />
<span style="color: #38761d;">1, 1998, 2</span><br />
<span style="color: #38761d;">1, 1999, 0</span><br />
<span style="color: #38761d;">1, 2999, 1</span><br />
<span style="color: #38761d;">...</span><br />
<span style="color: #38761d;"><br /></span><br />
Next step was to simply join my results to the unique_outcomes table and select null joined values. A simple and powerful function; and yes, I could have just saved my failed ids in my try, except block, but then I wouldn't have anything to post, ha.<br />
<br /></div>Anonymoushttp://www.blogger.com/profile/16259376360116803538noreply@blogger.com0tag:blogger.com,1999:blog-9041577666726226947.post-78114679005960245032011-10-22T20:00:00.001-06:002011-10-22T20:28:05.742-06:00PostGIS Nearest Neighbor Speed Tests<br />
<div class="p1">
I've used several different techniques for finding the nearest neighbor in PostGIS but never really knew which was the fastest. Below are some speed tests using 3 different approaches, all limited the spatial analysis to within 200 miles. Test 1 was the slowest, but returns the entire dataset with Nulls for rows that did not meet the maximum distance requirement. Test 2 was consistently the fastest and used a PostgreSQL Windows Function to rank by distance, it can also be used to return X many nearest neighbors, but requires a little more logic than the others. Test 3 uses a combination of Distinct and Order By to find the nearest neighbor, the advantage here is that if you would like to return the distance, you only have to calculate it once as it is required in your primary Select statement. Test results are in bold green and are in milliseconds. Obviously you'll have different reasons for using each, and I still need to test Boston GISs pg<i>fn</i>nn function, but hopefully this helps give you some guidance</div>
<br />
<ul>
</ul>
The tests were run against 43,000 point features in the source dataset and 20 polyline features in the neighboring dataset.<br />
<br />
<span class="Apple-style-span" style="color: #6aa84f;"><b>--Test 1 1882, 1729, 1820, 1778, 1597, 1595</b></span><br />
<span class="Apple-style-span" style="font-size: x-small;">SELECT</span><br />
<span class="Apple-style-span" style="font-size: x-small;"><span class="Apple-tab-span" style="white-space: pre;"> </span>a.gid,</span><br />
<span class="Apple-style-span" style="font-size: x-small;"><span class="Apple-tab-span" style="white-space: pre;"> </span>(SELECT b.gid</span><br />
<span class="Apple-style-span" style="font-size: x-small;"><span class="Apple-tab-span" style="white-space: pre;"> </span>FROM epa_iag.wecc_proposed b</span><br />
<span class="Apple-style-span" style="font-size: x-small;"><span class="Apple-tab-span" style="white-space: pre;"> </span>WHERE ST_DWithin(a.the_geom, b.the_geom, 200*1609.344) AND a.fallow = 3</span><br />
<span class="Apple-style-span" style="font-size: x-small;"><span class="Apple-tab-span" style="white-space: pre;"> </span>ORDER BY ST_Distance(a.the_geom, b.the_geom) LIMIT 1) AS b_gid</span><br />
<span class="Apple-style-span" style="font-size: x-small;">FROM</span><br />
<span class="Apple-style-span" style="font-size: x-small;"><span class="Apple-tab-span" style="white-space: pre;"> </span>epa_iag.fallow_lands_cnty a</span><br />
<span class="Apple-style-span" style="font-size: x-small;">WHERE</span><br />
<span class="Apple-style-span" style="font-size: x-small;"><span class="Apple-tab-span" style="white-space: pre;"> </span>a.fallow = 3</span><br />
<br />
<b><span class="Apple-style-span" style="color: #6aa84f;">--Test 2 1667, 857, 1612, 694, 755, 702</span></b><br />
<span class="Apple-style-span" style="font-size: x-small;">SELECT</span><br />
<span class="Apple-style-span" style="font-size: x-small;"><span class="Apple-tab-span" style="white-space: pre;"> </span>i.gid,</span><br />
<span class="Apple-style-span" style="font-size: x-small;"><span class="Apple-tab-span" style="white-space: pre;"> </span>i.b_gid</span><br />
<span class="Apple-style-span" style="font-size: x-small;"><span class="Apple-tab-span" style="white-space: pre;"> </span>--ST_Distance(i.the_geom, i.b_geom) / 1609.344 AS dist</span><br />
<span class="Apple-style-span" style="font-size: x-small;">FROM(</span><br />
<span class="Apple-style-span" style="font-size: x-small;">SELECT</span><br />
<span class="Apple-style-span" style="font-size: x-small;"><span class="Apple-tab-span" style="white-space: pre;"> </span>a.gid,</span><br />
<span class="Apple-style-span" style="font-size: x-small;"><span class="Apple-tab-span" style="white-space: pre;"> </span>b.gid AS b_gid,</span><br />
<span class="Apple-style-span" style="font-size: x-small;"><span class="Apple-tab-span" style="white-space: pre;"> </span>a.the_geom, b.the_geom AS b_geom,</span><br />
<span class="Apple-style-span" style="font-size: x-small;"><span class="Apple-tab-span" style="white-space: pre;"> </span>rank() OVER (PARTITION BY a.gid ORDER BY</span><br />
<span class="Apple-style-span" style="font-size: x-small;"><span class="Apple-tab-span" style="white-space: pre;"> </span>ST_Distance(a.centroid_geom, b.the_geom) / 1609.344 ASC) AS pos</span><br />
<span class="Apple-style-span" style="font-size: x-small;">FROM</span><br />
<span class="Apple-style-span" style="font-size: x-small;"><span class="Apple-tab-span" style="white-space: pre;"> </span>epa_iag.wecc_proposed b,</span><br />
<span class="Apple-style-span" style="font-size: x-small;"><span class="Apple-tab-span" style="white-space: pre;"> </span>epa_iag.fallow_lands_cnty a</span><br />
<span class="Apple-style-span" style="font-size: x-small;">WHERE</span><br />
<span class="Apple-style-span" style="font-size: x-small;"><span class="Apple-tab-span" style="white-space: pre;"> </span>ST_DWithin(a.the_geom, b.the_geom, 200*1609.344) AND a.fallow = 3 ) i</span><br />
<span class="Apple-style-span" style="font-size: x-small;">WHERE</span><br />
<span class="Apple-style-span" style="font-size: x-small;"><span class="Apple-tab-span" style="white-space: pre;"> </span>pos = 1</span><br />
<span class="Apple-tab-span" style="white-space: pre;"> </span><br />
<b><span class="Apple-style-span" style="color: #6aa84f;">--Test 3 2409, 1057, 1192, 1098, 1089, 2167</span></b><br />
<span class="Apple-style-span" style="font-size: x-small;">SELECT<span class="Apple-tab-span" style="white-space: pre;"> </span>DISTINCT ON (a.gid)</span><br />
<span class="Apple-style-span" style="font-size: x-small;"><span class="Apple-tab-span" style="white-space: pre;"> </span>a.gid,</span><br />
<span class="Apple-style-span" style="font-size: x-small;"><span class="Apple-tab-span" style="white-space: pre;"> </span>b.gid AS b_gid,</span><br />
<span class="Apple-style-span" style="font-size: x-small;"><span class="Apple-tab-span" style="white-space: pre;"> </span>ST_Distance(a.centroid_geom, b.the_geom) / 1609.344 AS dist</span><br />
<span class="Apple-style-span" style="font-size: x-small;">FROM</span><br />
<span class="Apple-style-span" style="font-size: x-small;"><span class="Apple-tab-span" style="white-space: pre;"> </span>epa_iag.wecc_proposed b,</span><br />
<span class="Apple-style-span" style="font-size: x-small;"><span class="Apple-tab-span" style="white-space: pre;"> </span>epa_iag.fallow_lands_cnty a</span><br />
<span class="Apple-style-span" style="font-size: x-small;">WHERE</span><br />
<span class="Apple-style-span" style="font-size: x-small;"><span class="Apple-tab-span" style="white-space: pre;"> </span>ST_DWithin(a.the_geom, b.the_geom, 200*1609.344) AND a.fallow = 3</span><br />
<span class="Apple-style-span" style="font-size: x-small;">ORDER BY</span><br />
<span class="Apple-style-span" style="font-size: x-small;"><span class="Apple-tab-span" style="white-space: pre;"> </span>a.gid, dist ASC;</span>Anonymoushttp://www.blogger.com/profile/16259376360116803538noreply@blogger.com0tag:blogger.com,1999:blog-9041577666726226947.post-76190862692170154572011-07-30T08:59:00.004-06:002011-08-02T20:46:06.322-06:00Connect GAMS to PostgreSQL-PostGIS (w/ sql2gms)First things first, if you are running a 32bit PostgreSQL/PostGIS on a 64bit Windows box, you will need to read my previous blog post (<a href="http://movingspatial.blogspot.com/2011/07/configure-32bit-psqlodbc-on-64bit.html">how to configure psqlODBC</a>) before moving on. Back? Ok, once you have your DSN setup you are ready to move onto the next step, testing your connection with sql2gms. I found the easiest way to test my connection was through the sql2gms GUI, located in the root directory of your GAMS install (i.e. c:\Program Files\GAMS\23.1\...). Once there...<br />
<br />
Click the Options button to set up your connection<br />
<div class="separator" style="clear: both; text-align: center;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEiVi3soSXvojSv_0Io9dHlif4sPjm9-vd3Zm6_H85FFX-lIQ5jiDaZolIjEZ5LYmpdRLOe8L4wzsoMlAwdzNQuuegST-6ZKhcorsDxDAvDZGUah348nuIpDdirroPmYT2nmrjLfUGC0Tm9K/s1600/Screen+shot+2011-07-30+at+8.04.45+AM.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" height="280" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEiVi3soSXvojSv_0Io9dHlif4sPjm9-vd3Zm6_H85FFX-lIQ5jiDaZolIjEZ5LYmpdRLOe8L4wzsoMlAwdzNQuuegST-6ZKhcorsDxDAvDZGUah348nuIpDdirroPmYT2nmrjLfUGC0Tm9K/s400/Screen+shot+2011-07-30+at+8.04.45+AM.png" width="400" /></a></div><br />
Enter your user name and password. If the connections string is not present, Select the ODBC Data Source/Driver to match mine below. The Connection String should auto populate, if not, replicate the input below. Test the connection, it should connect then disconnect (as seen below). Success! Click ok<br />
<div class="separator" style="clear: both; text-align: center;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjLzLMnBh4RFT9BMzheRKEFkFrQiOT7h0THtnKIEhuJi6pvUZaKatY-f3hMTB0NlPQNs3u3xyp0Ev7nxAeSDNpWPeX0b4DHZBRt5vD3v6fDiJohLfIK9Ewn8JPTNMymQBE1sNSCo6OvVcvc/s1600/Screen+shot+2011-07-30+at+8.05.25+AM.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" height="292" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjLzLMnBh4RFT9BMzheRKEFkFrQiOT7h0THtnKIEhuJi6pvUZaKatY-f3hMTB0NlPQNs3u3xyp0Ev7nxAeSDNpWPeX0b4DHZBRt5vD3v6fDiJohLfIK9Ewn8JPTNMymQBE1sNSCo6OvVcvc/s400/Screen+shot+2011-07-30+at+8.05.25+AM.png" width="400" /></a></div><br />
Now do a test query. Success!<br />
<div class="separator" style="clear: both; text-align: center;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEiaWxx0GyArJPGDff4leJmmeR0B8gqZnZ0s8jEfMf1XWCTR3U1WOoIMfqrZ0FgAxew4G8it54ayWStNDu2UrUMAXqRA5YDEwADMQfAuzJgHeuvg6UiMEHSNwwgWHrEPoCJZ7NxuFYfrCsOZ/s1600/Screen+shot+2011-07-30+at+8.06.39+AM.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" height="280" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEiaWxx0GyArJPGDff4leJmmeR0B8gqZnZ0s8jEfMf1XWCTR3U1WOoIMfqrZ0FgAxew4G8it54ayWStNDu2UrUMAXqRA5YDEwADMQfAuzJgHeuvg6UiMEHSNwwgWHrEPoCJZ7NxuFYfrCsOZ/s400/Screen+shot+2011-07-30+at+8.06.39+AM.png" width="400" /></a></div>Now that you know your connection works you can program away. The sql2gms documentation via GAMS/help is quite good.<br />
<br />
Below you'll see the standard GAMS tutorial transportation problem, but with sql2gms modifications. Q1 represents the 'Spatially Enabled' part. Q1 is Parameter 'd', and ST_Distance() is used to calculate the distance between the cities. Q2 is Parameter i, etc. <br />
<br />
<b><span style="font-size: x-small;"><span style="color: #666666;">*Create cmd text to execute multiple queries</span><br />
$onecho > cmd.txt<br />
<span style="color: #666666;">C="DSN=PostgreSQL35W"</span><br style="color: #666666;" /><span style="color: #666666;">X="all_data.gdx"</span><br style="color: #666666;" /><br style="color: #666666;" /><span style="color: #666666;">Q1="SELECT b.city_name AS plants, a.city_name AS markets, CAST(ST_Distance(a.the_geom, b.the_geom) / 1609.344 * 0.001 AS numeric(8,2)) AS dist FROM (SELECT city_name, the_geom FROM gams.city_comp WHERE city_type = 'market') a, (SELECT city_name, the_geom FROM gams.city_comp WHERE city_type = 'plant' ) b WHERE a.city_name <> b.city_name"</span><br style="color: #666666;" /><span style="color: #666666;">A1=d</span><br style="color: #666666;" /><br style="color: #666666;" /><span style="color: #666666;">Q2="SELECT city_name FROM gams.city_comp WHERE city_type = 'plant';"</span><br style="color: #666666;" /><span style="color: #666666;">S2=i</span><br style="color: #666666;" /><br style="color: #666666;" /><span style="color: #666666;">Q3="SELECT city_name FROM gams.city_comp WHERE city_type = 'market';"</span><br style="color: #666666;" /><span style="color: #666666;">S3=j</span><br style="color: #666666;" /><br style="color: #666666;" /><span style="color: #666666;">Q4="SELECT city_name, capacity FROM gams.city_comp WHERE city_type = 'plant'"</span><br style="color: #666666;" /><span style="color: #666666;">A4=a</span><br style="color: #666666;" /><br style="color: #666666;" /><span style="color: #666666;">Q5="SELECT city_name, demand FROM gams.city_comp WHERE city_type = 'market'"</span><br style="color: #666666;" /><span style="color: #666666;">A5=b</span><br />
$offecho<br />
<br />
<span style="color: #666666;">*Execute query</span><br />
$call =sql2gms @cmd.txt<br />
$call =gdxviewer all_data.gdx<br />
<br />
<span style="color: #666666;">*Create sets and Parameters</span><br />
<span style="color: #3d85c6;">SETS</span><br />
i canning plants<br />
j markets ;<br />
<br />
<span style="color: #3d85c6;">Parameters</span><br />
a(i) capacity of plant i in cases<br />
b(j) demand at market j in cases<br />
d(i,j) distance between plants and markets ;<br />
<br />
<span style="color: #666666;">*Load all the data from the gdx</span><br />
$gdxin all_data.gdx<br />
$load i j d a b<br />
$gdxin<br />
<br />
<span style="color: #3d85c6;">Scalar </span>f freight in dollars per case per thousand miles <span style="color: #38761d;">/90/</span> ;<br />
<br />
<span style="color: #3d85c6;">Parameter </span>c(i,j) transport cost in thousands of dollars per case ;<br />
<br />
c(i,j) = f * d(i,j) / 1000 ;<br />
<br />
<span style="color: #3d85c6;">Variables</span><br />
x(i,j) shipment quantities in cases<br />
z total transportation costs in thousands of dollars ;<br />
<br />
<span style="color: #3d85c6;">Positive Variable </span>x ;<br />
<br />
<span style="color: #3d85c6;">Equations</span><br />
cost define objective function<br />
supply(i) observe supply limit at plant i<br />
demand(j) satisfy demand at market j ;<br />
<br />
<span style="color: #666666;">*Constraints</span><br />
supply(i) .. sum(j, x(i,j)) =l= a(i) ;<br />
demand(j) .. sum(i, x(i,j)) =g= b(j) ;<br />
<br />
<span style="color: #666666;">*objective</span><br />
cost .. z =e= sum((i,j), c(i,j) * x(i,j)) ;<br />
<br />
<span style="color: #3d85c6;">Model </span>transport <span style="color: #38761d;">/all/</span> ;<br />
<br />
<span style="color: #3d85c6;">Solve </span>transport using lp minimizing z ;<br />
<br />
<span style="color: #3d85c6;">Display </span>x.l, x.m ; </span></b>Anonymoushttp://www.blogger.com/profile/16259376360116803538noreply@blogger.com0tag:blogger.com,1999:blog-9041577666726226947.post-78169512096340950082011-07-29T21:13:00.001-06:002011-07-29T21:15:52.794-06:00Configure 32bit psqlODBC on 64bit Windows Server 2008Installing and configuring psqlODBC was not as easy as I had anticipated. I guess, I was amazed not many other folks were having the same issues I had. In the end, when all the pieces are in one place, its really not that difficult. As the title said, this post walks you through the steps for installing 32bit psqlODBC on a 64bit Windows Server 2008. Why you ask install 32bit on a 64bit? Simply because postGIS is 32bit only.<br />
<br />
<ul><li>Install psqlODBC </li>
<ul><li>Unfortunately I kept getting errors using the Stack Builder so I went to the <a href="http://www.postgresql.org/ftp/odbc/versions/">source</a></li>
</ul></ul><ul><li>Add ODBC to Data Source Administrator (tricky part)</li>
<ul><li> You cannot add a 32bit driver using the 64bit Data Source Admin</li>
<li>To activate 32bit Data Source Admin:</li>
<ul><li> Go to Start</li>
<li>Run</li>
<li>type: %WINDIR%\SysWOW64\odbcad32.exe</li>
</ul><li>Add System DSN, select PostgreSQL Unicode and enter the server parameters </li>
</ul></ul>Now, that wasn't so bad now was it?Anonymoushttp://www.blogger.com/profile/16259376360116803538noreply@blogger.com1tag:blogger.com,1999:blog-9041577666726226947.post-58328156472879708972011-06-21T08:20:00.001-06:002011-06-22T08:28:19.442-06:00Setup VirtualEnv & Psycopg2 on Linux - redhatWe recently setup a boss of a computer at work to handle large analysis tasks. Luckily I got out of the initial setup, but when it came time to setup my virtualenv with the usual, I ran into a roadblock I had never encountered; pg_config cannot be found. So, here is the low down on how to setup your own virtualenv with psycopg2 on redhat.<br />
<br />
<span style="font-size: x-small;">-Create python virtualenv in your home director ("pymain" is what I call my virtualenv, you can call yours whatever you like):</span><br />
<span style="font-size: x-small;"> <b><span style="color: #38761d;">virtualenv --no-site-packages pymain </span></b><br />
-To activate your virtualenv, first navigate to the "bin" folder you just created: <b style="color: #38761d;"> </b></span><br />
<span style="font-size: x-small;"><b style="color: #38761d;">cd pymain/bin/ </b><br />
-Then activate using (you'll notice in the terminal the name of your virtualenv appears before yours after doing this): </span><br />
<span style="font-size: x-small;"><b><span style="color: #38761d;">source activate</span></b><br />
- Setup symbolic link to pg_config (so psycopg2 knows where to find it. THIS WAS THE ROADBLOCK): </span><br />
<span style="font-size: x-small;"><b style="color: #38761d;">ln –s /usr/local/pgsql/bin/pg_config ~/pymain/bin/pg_config</b><br />
-Install psycopg2:</span><br />
<span style="font-size: x-small;"> <b><span style="color: #38761d;">easy_install psycopg2</span></b><br />
-To deactivate your virtualenv simply type:</span><br />
<span style="font-size: x-small;"><b style="color: #38761d;">deactivate</b></span>Anonymoushttp://www.blogger.com/profile/16259376360116803538noreply@blogger.com0tag:blogger.com,1999:blog-9041577666726226947.post-7295504946534150972011-05-24T19:02:00.001-06:002011-06-21T09:58:24.553-06:00Batch Import Shapefiles into PostgreSQLImporting shapefiles into PostgreSQL is simple using shp2pgsql, unless you want batch import multiple shps. Bash to the rescue. Below is a simple script that imports all shapefiles within a directory.<br />
<br />
To run:<br />
<ul><li>Update the schema in the script</li>
<li> Update the SRID (make sure all shps in the folder have the same SRID)</li>
<li>Update the database name and user </li>
<li>Navigate to the folder where your data is sitting</li>
<li>Decide whether or not you want to remove the shps afterwards (if not, remove the line above '<i>done</i>') </li>
<li>Run</li>
</ul>#-------------------------------------------------------------- <br />
<div style="color: #38761d;"><span style="font-size: small;">FILES=$(find *.shp);</span><br />
<span style="font-size: small;">for FILE in $FILES</span><br />
<span style="font-size: small;"> do</span><br />
<span style="font-size: small;"> basename=${FILE%.shp}</span><br />
<span style="font-size: small;"> shp2pgsql -D -s 4326 -I $FILE schema.$basename | psql -U postgres -d postgis</span><br />
<span style="font-size: small;"> rm $basename.*</span><br />
<span style="font-size: small;">done</span></div>#--------------------------------------------------------------Anonymoushttp://www.blogger.com/profile/16259376360116803538noreply@blogger.com0tag:blogger.com,1999:blog-9041577666726226947.post-66151523032634584252011-05-18T18:52:00.006-06:002012-10-27T17:20:22.474-06:00PostGIS Polygon to Point GridI've occasionally run into the need to convert my PostGIS polygon geometries to points. Usually this need arises as my datasets get large and larger and simple geometry helps speed up processing (treating the dataset as a raster). Hopefully, with the release of PostGIS 2.0 this function will lose its value (hopefully), since 2.0 will fully support rasters.<br />
<br />
To the function<br />
<br />
The below function converts polygons to point grids letting the user specify the resolution. It also takes on the projection of the polygonal dataset. Simple yet effective. Note, that I cannot take full credit for this function, I had seen the beginnings of it on some user forum. Unfortunately I simply cannot track it down...if I stumble into it again I will post the url. Note, that the resolution is in projected units i.e. if your projection is in meters, then inputting 10 = 10 meter resolution. <br />
<div style="color: #0b5394;">
<br />
<table align="center" cellpadding="0" cellspacing="0" class="tr-caption-container" style="margin-left: auto; margin-right: auto; text-align: center;"><tbody>
<tr><td style="text-align: center;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgxwun1KdYf1ahlMBpPS09l9UXJYq94ydxog7E_f-GHrqYmJU0EKe4I19HJFTWFXpok_BeN2KIKGS-vt3o-jUaRgLM7FT8ilgRi4GGXu5EXgg2c_vv_DvsCYfzol1tn4heABt20TDccQk6X/s1600/grid.jpg" imageanchor="1" style="margin-left: auto; margin-right: auto;"><img border="0" height="200" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgxwun1KdYf1ahlMBpPS09l9UXJYq94ydxog7E_f-GHrqYmJU0EKe4I19HJFTWFXpok_BeN2KIKGS-vt3o-jUaRgLM7FT8ilgRi4GGXu5EXgg2c_vv_DvsCYfzol1tn4heABt20TDccQk6X/s200/grid.jpg" width="200" /></a></td></tr>
<tr><td class="tr-caption" style="text-align: center;"><i>point grid inside polygon</i></td></tr>
</tbody></table>
<div style="color: black;">
<br /></div>
</div>
<div style="color: #666666;">
<span style="color: black;">-----------------------------------------------------------------------------------------</span><br />
<span style="font-size: x-small;"><b><span style="color: #0b5394;">CREATE OR REPLACE FUNCTION</span> st_polygrid(<span style="color: #0b5394;">geometry</span>, <span style="color: #0b5394;">integer</span>) <span style="color: #0b5394;">RETURNS geometry AS</span><br />
$$</b></span></div>
<div style="color: #666666;">
<span style="font-size: x-small;"><b><span style="color: #0b5394;">SELECT </span><br />
ST_SetSRID(ST_Collect(ST_POINT(x,y)), ST_SRID($1))<br />
<span style="color: #0b5394;">FROM </span><br />
generate_series(floor(ST_XMin($1))::</b></span><b style="font-size: small;">numeric</b><b style="font-size: small;">, ceiling(ST_xmax($1))::numeric, $2) as x,</b><br />
<span style="font-size: x-small;"><b>
generate_series(floor(ST_ymin($1))::</b></span><b style="font-size: small;">numeric</b><b style="font-size: small;">, ceiling(ST_ymax($1))::numeric,$2) as y </b><br />
<span style="font-size: x-small;"><b>
<span style="color: #0b5394;">WHERE</span><br />
ST_Intersects($1,ST_SetSRID(ST_POINT(x,y), ST_SRID($1)))</b></span></div>
<div style="color: #666666;">
<span style="font-size: x-small;"><b>$$<br />
<span style="color: #0b5394;">LANGUAGE</span> sql <span style="color: #0b5394;">VOLATILE<span style="color: #666666;">;</span></span></b></span><br />
<div style="color: #38761d;">
<span style="color: black;">-----------------------------------------------------------------------------------------</span><br />
<br /></div>
<div style="color: #38761d;">
<span style="font-size: x-small;"><b>--Example using the above function</b></span></div>
<span style="font-size: x-small;"><b><span style="color: #0b5394;"><span style="color: #666666;"><span style="color: #0b5394;">CREATE TABLE</span> some.other_table <span style="color: #0b5394;">AS</span></span></span></b></span><br />
<div style="color: #0b5394;">
<span style="font-size: x-small;"><b>SELECT</b></span></div>
<span style="font-size: x-small;"><b><span style="color: #0b5394;"><span style="color: #666666;"> a.gid,</span></span></b></span><br />
<span style="font-size: x-small;"><b><span style="color: #0b5394;"><span style="color: #666666;"> st_polygrid(a.the_geom, 1000) <span style="color: #0b5394;">AS </span>the_geom</span></span></b></span><br />
<div style="color: #0b5394;">
<span style="font-size: x-small;"><b>FROM</b></span></div>
<span style="font-size: x-small;"><b><span style="color: #0b5394;"><span style="color: #666666;"> some.table a</span></span></b></span></div>
Anonymoushttp://www.blogger.com/profile/16259376360116803538noreply@blogger.com0tag:blogger.com,1999:blog-9041577666726226947.post-5557448603040552362011-05-18T10:58:00.004-06:002011-05-19T07:15:26.701-06:00Google Geocoder in PostGIS (with reverse geocode!)Well, it's been a while since my last post. It seems there is a negative correlation between having twins and blog posts! <br />
<br />
Alright, onto the reason for this post, PostGIS geocoding using Google! OK, I realize that there are other methods for geocoding in PostGIS, but I had a couple of drivers motivating me: <br />
<b>1.</b> I wanted to create a Function in Postgres using PL-Python. A way of expanding my horizons we'll say.<br />
<b>2.</b> Geopy (get this python module <a href="http://code.google.com/p/geopy/">here</a>) makes it extremely easy to geocode using Google, Yahoo, GeoNames, MediaWiki and more. <br />
<b>3.</b> I wanted to be able to reverse geocode i.e. give a lat, long and return an address. <br />
<br />
You'll need a couple things first before these functions will work (don't worry, nothing drastic!). <br />
<b>1.</b> Make sure you have PL-Python installed. To check if you have it; if you have pgAdmin, go to Preferences, and navigate to the list of objects displayed in the tree and check Languages. If its not there, pgAdmin to the rescue. Right click languages, new language, plpythonu, done!<br />
<b>2.</b> Make sure you have Geopy installed, see link above. If you operate in a virtualenv, no worries, I've added a simple workaround for that in the function. So virtualenv away!<br />
<b>3.</b> Get a Google API key.<br />
<b>4.</b> If you are planning on using the reverse geocoder, you'll need to download and install the experimental (but stable) branch of geopy <a href="http://code.google.com/p/geopy/wiki/ReverseGeocoding">here</a>. <br />
<br />
Now for the code!<br />
-------------------------------------------------------------------------------<br />
<div style="color: #38761d;"><b><span style="font-size: x-small;">--Converts an address to lat, long</span></b></div><div style="color: #666666;"><b><span style="font-size: x-small;"><span style="color: #0b5394;">CREATE OR REPLACE FUNCTION</span> google_latlng(<span style="color: #0b5394;">IN </span>address <span style="color: #0b5394;">text</span>)<span style="color: #0b5394;"> RETURNS text ARRAY[2] AS</span></span></b></div><div style="color: #666666;"><b><span style="font-size: x-small;">$$</span></b></div><div style="color: #666666;"><b><span style="font-size: x-small;">#If using virtualenv set path to appropriate environment</span></b></div><div style="color: #666666;"><b><span style="font-size: x-small;">import site</span></b></div><div style="color: #666666;"><b><span style="font-size: x-small;">site.addsitedir('/your/virtualenv/lib/pythonX.Y/site-packages')</span></b></div><div style="color: #666666;"><b><span style="font-size: x-small;"><br />
</span></b></div><div style="color: #666666;"><b><span style="font-size: x-small;">from geopy import geocoders</span></b></div><div style="color: #666666;"><b><span style="font-size: x-small;">gc = geocoders.Google('Your key here')</span></b></div><div style="color: #666666;"><b><span style="font-size: x-small;">try:</span></b></div><div style="color: #666666;"><b><span style="font-size: x-small;"> geo = list(gc.geocode(address, exactly_one=True))</span></b></div><div style="color: #666666;"><b><span style="font-size: x-small;"> lat, lng = geo[1][0], geo[1][1]</span></b></div><div style="color: #666666;"><b><span style="font-size: x-small;">except:</span></b></div><div style="color: #666666;"><b><span style="font-size: x-small;"> lat, lng = 0, 0</span></b></div><div style="color: #666666;"><b><span style="font-size: x-small;"><br />
</span></b></div><div style="color: #666666;"><b><span style="font-size: x-small;"><span style="color: #0b5394;">return </span>lng, lat</span></b></div><div style="color: #666666;"><b><span style="font-size: x-small;">$$</span></b></div><div style="color: #666666;"><b><span style="font-size: x-small;"><span style="color: #0b5394;">LANGUAGE </span>'plpythonu' <span style="color: #0b5394;">VOLATILE</span>;</span></b></div><div style="color: #666666;"><b><span style="font-size: x-small;"><br />
</span></b></div><div style="color: #38761d;"><b><span style="font-size: x-small;">--Converts an address to PostGIS geometry</span></b></div><div style="color: #666666;"><b><span style="font-size: x-small;"><span style="color: #0b5394;">CREATE OR REPLACE FUNCTION</span> google_geocode(<span style="color: #0b5394;">IN </span>address <span style="color: #0b5394;">text</span>)<span style="color: #0b5394;"> RETURNS geometry AS</span></span></b></div><div style="color: #666666;"><b><span style="font-size: x-small;">$$</span></b></div><div style="color: #0b5394;"><b><span style="font-size: x-small;">SELECT</span></b></div><div style="color: #666666;"><b><span style="font-size: x-small;">ST_SetSRID(ST_GeomFromText('POINT(' || r.ll[1] || ' ' || r.ll[2] || ')'), 4326)</span></b></div><div style="color: #666666;"><b><span style="font-size: x-small;"><span style="color: #0b5394;">FROM </span>(</span></b></div><div style="color: #0b5394;"><b><span style="font-size: x-small;">SELECT</span></b></div><div style="color: #666666;"><b><span style="font-size: x-small;">google_latlng($1) AS ll ) r</span></b></div><div style="color: #666666;"><b><span style="font-size: x-small;">$$</span></b></div><div style="color: #666666;"><b><span style="font-size: x-small;"><span style="color: #0b5394;">LANGUAGE </span>sql <span style="color: #0b5394;">VOLATILE</span>;</span></b></div><div style="color: #666666;"><b><span style="font-size: x-small;"><br />
</span></b></div><div style="color: #38761d;"><b><span style="font-size: x-small;">--Reverse Geocoding! Converts lat, long to an address</span></b></div><div style="color: #666666;"><b><span style="font-size: x-small;"><span style="color: #0b5394;">CREATE OR REPLACE FUNCTION</span> google_addr(<span style="color: #0b5394;">IN</span> lat <span style="color: #0b5394;">double precision</span>, lng <span style="color: #0b5394;">double precision</span>)<span style="color: #0b5394;"> RETURNS text AS</span></span></b></div><div style="color: #666666;"><b><span style="font-size: x-small;">$$</span></b></div><div style="color: #666666;"><b><span style="font-size: x-small;">#If using virtualenv set path to appropriate environment</span></b></div><div style="color: #666666;"><b><span style="font-size: x-small;">import site</span></b></div><div style="color: #666666;"><b><span style="font-size: x-small;">site.addsitedir('</span></b><b><span style="font-size: x-small;">/your/virtualenv/lib/pythonX.Y/site-packages</span></b><b><span style="font-size: x-small;">')</span></b></div><div style="color: #666666;"><b><span style="font-size: x-small;"><br />
</span></b></div><div style="color: #666666;"><b><span style="font-size: x-small;">from geopy import geocoders</span></b></div><div style="color: #666666;"><b><span style="font-size: x-small;">gc = geocoders.Google('Your key here')</span></b></div><div style="color: #666666;"><b><span style="font-size: x-small;">try:</span></b></div><div style="color: #666666;"><b><span style="font-size: x-small;"> (address, point) = gc.reverse((lat,lng))</span></b></div><div style="color: #666666;"><b><span style="font-size: x-small;">except:</span></b></div><div style="color: #666666;"><b><span style="font-size: x-small;"> address = '????'</span></b></div><div style="color: #666666;"><b><span style="font-size: x-small;"><br />
</span></b></div><div style="color: #666666;"><b><span style="font-size: x-small;"><span style="color: #0b5394;">return</span> address</span></b></div><div style="color: #666666;"><b><span style="font-size: x-small;">$$</span></b></div><div style="color: #666666;"><b><span style="font-size: x-small;"><span style="color: #0b5394;">LANGUAGE</span> 'plpythonu' <span style="color: #0b5394;">VOLATILE</span>;</span></b></div><br />
------------------------------------------------------------------------------------<br />
<br />
<br />
That's it! Create those functions and your on your way. Just remember to update your Geography Columns Table if you create geometry.<br />
<br />
And here are some examples using the functions:<br />
<div style="color: #38761d;">--Reverse Geocode a lat, long</div><div style="color: #666666;"><span style="font-size: x-small;"><span style="color: #0b5394;">SELECT</span> google_addr(39.7439848, -105.0213316</span><span style="font-size: x-small;">)</span></div><br />
<div style="color: #38761d;">--Return the lat, long of an address as a text array</div><div style="color: #666666;"><span style="font-size: x-small;"><span style="color: #0b5394;">SELECT </span>google_latlng('1701 Mile High Stadium # 300 Denver, CO 80204')</span></div><br />
<div style="color: #38761d;">--Return a PostGIS geometry given an address, note this is dependent on google_latlng()</div><div style="color: #666666;"><span style="font-size: x-small;"><span style="color: #0b5394;">SELECT </span>google_geocode('1701 Mile High Stadium # 300 Denver, CO 80204')</span></div><br />
Enjoy!Anonymoushttp://www.blogger.com/profile/16259376360116803538noreply@blogger.com3tag:blogger.com,1999:blog-9041577666726226947.post-58801835979597490322010-11-24T21:42:00.003-07:002011-05-19T07:17:01.079-06:00Reclassify Shapefile using PythonReclassifying data is pretty standard in the GIS world. In the Raster realm using ESRI its pretty easy, simply create a reclass table. But occasionally (in my case often) you'll find yourself needing to reclassify vector data. Enter Python to help save your sanity. Below you'll find a simple function that accepts a value and Python dictionary. To use, first create a function to define a dictionary (your reclass table) then simply iterate through your shapfile and pass the value, and dictionary. Easy and simple but what a time saver...<br />
<br />
<pre style="background-image: URL(https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEitTzQYpR3KnwgrimJHF2i-bLvCl-N0qx7_09T6Pzu3lZiCVXnlY1g4aZOMK-YRyKXOLGh8DBYDfUm_ed1IEqT3Sth-Xs-9L79Rl2FS7gKrXIGPNEXkhqFbJZ6n5g0PnLiZFHC3JMmuQqS7/s320/codebg.gif); background: #f0f0f0; border: 1px dashed #CCCCCC; color: black; font-family: arial; font-size: 12px; height: auto; line-height: 20px; overflow: auto; padding: 0px; text-align: left; width: 99%;"><code style="color: black; word-wrap: normal;"> import arcgisscripting
gp = arcgiscripting.Create(9.3)
# DEFINE THE RECLASS TABLE
def reclassTable():
myDict:{'a': range(0, 11),
'b': range(11, 21),
'c': range(21, 31),
'd': range(31, 100)}
return myDict
# THE RECLASS FUNCTION
def reclass_Function(value, classDict):
for cls, rng in classDict.iteritems():
if value in rng:
return cls
return -9999
# THE READ AND WRITE FUNCTION
def classify_Data(inputSHP, oldFieldVals, newField):
gp.addField(inputSHP, newField, 'TEXT', 5)
rows = gp.UpdateCursor(inputSHP)
row = rows.Next()
while row:
row.SetValue(newField, reclass_Function(row.GetValue(oldFieldVals), reclassTable()))
rows.UpdateRow(row)
row = rows.Next()
# PASS THE VALUES
if__name__=='__main__':
inputSHP = 'C:/temp/myShapefile.shp'
classify_Data(inputSHP, 'KWh', 'Class')
</code></pre>Now simply edit the dictionary to fit your needs and pass the values accordingly.Anonymoushttp://www.blogger.com/profile/16259376360116803538noreply@blogger.com0tag:blogger.com,1999:blog-9041577666726226947.post-6633130161012974382010-11-20T09:21:00.005-07:002011-05-19T07:22:05.400-06:00The Hydrogen ViewerThe Hydrogen Viewer was a project I worked on several months ago. I was tasked with providing a visualization of hydrogen station installations over time. Although spatial-temporal in nature it still was not the most exciting data...so I decided to spice it up using Flex, Google Maps and a couple custom components and features. I'll cover a couple of them here.<br />
<br />
Having a reference or background is key to any map but sometimes the reference layer can take away from the data you are really trying to highlight. Enter the gray-scale mat. Applying a grayscale mat to your Google Map really makes your data pop out and it is easy. Simply apply a color matrix filter in your "onMapReady" function.<br />
<div style="text-align: center;"> <span style="font-size: large;"><b>APPLY A GRAYSCALE TO A GOOGLE MAP </b></span></div><pre style="background: url("https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEitTzQYpR3KnwgrimJHF2i-bLvCl-N0qx7_09T6Pzu3lZiCVXnlY1g4aZOMK-YRyKXOLGh8DBYDfUm_ed1IEqT3Sth-Xs-9L79Rl2FS7gKrXIGPNEXkhqFbJZ6n5g0PnLiZFHC3JMmuQqS7/s320/codebg.gif") repeat scroll 0% 0% rgb(240, 240, 240); border: 1px dashed rgb(204, 204, 204); color: black; font-family: arial; font-size: 12px; height: auto; line-height: 20px; overflow: auto; padding: 0px; text-align: left; width: 99%;"><code style="color: black; word-wrap: normal;">1: private function onMapReady(event:Event):void{
2: map.setCenter(new LatLng(38, -110), 3, MapType.NORMAL_MAP_TYPE);
3: var mat:Array = [0.3086,0.6094,0.082,0,-100,0.3086,0.6094,0.082,0,-100,0.3086,0.6094,0.082,0,-100,0,0,0,1,0];
4: var colorMat: ColorMatrixFilter = new ColorMatrixFilter(mat);
5: var s1: Sprite = map.getChildAt(1) as Sprite;
6: var s2: Sprite = s1.getChildAt(0) as Sprite;
7: s2.filters = [colorMat];}
</code></pre><br />
Creating a custom marker info window allows you maintain a consistent theme throughout your app. You can accomplish this with the following (I used one function for all of this):<br />
<br />
<div style="text-align: center;"><span style="font-size: large;"><b>CREATE CUSTOM INFO WINDOW</b></span><br />
<span style="font-size: large;"><b> </b></span> </div>First, create your marker. I generally use custom markers:<br />
<pre style="background: url("https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEitTzQYpR3KnwgrimJHF2i-bLvCl-N0qx7_09T6Pzu3lZiCVXnlY1g4aZOMK-YRyKXOLGh8DBYDfUm_ed1IEqT3Sth-Xs-9L79Rl2FS7gKrXIGPNEXkhqFbJZ6n5g0PnLiZFHC3JMmuQqS7/s320/codebg.gif") repeat scroll 0% 0% rgb(240, 240, 240); border: 1px dashed rgb(204, 204, 204); color: black; font-family: arial; font-size: 12px; height: auto; line-height: 20px; overflow: auto; padding: 0px; text-align: left; width: 99%;"><code style="color: black; word-wrap: normal;">1: var point:Sprite = new Sprite();
2: point.graphics.beginFill(color, .9);
3: point.graphics.lineStyle(1, 0xFF2F4F4F, 1);
4: point.graphics.drawCircle(0, 0, size);
5: point.graphics.endFill();
6:
7: markerOptions:MarkerOptions = new MarkerOptions();
8: markerOptions.icon = point;
9: markerOptions.clickable = true;
10: markerOptions.hasShadow = false;
11: marker = new Marker(new LatLng(yPosition,xPosition), markerOptions);
</code></pre><br />
Then create the content of the window:<br />
<pre style="background: url("https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEitTzQYpR3KnwgrimJHF2i-bLvCl-N0qx7_09T6Pzu3lZiCVXnlY1g4aZOMK-YRyKXOLGh8DBYDfUm_ed1IEqT3Sth-Xs-9L79Rl2FS7gKrXIGPNEXkhqFbJZ6n5g0PnLiZFHC3JMmuQqS7/s320/codebg.gif") repeat scroll 0% 0% rgb(240, 240, 240); border: 1px dashed rgb(204, 204, 204); color: black; font-family: arial; font-size: 12px; height: auto; line-height: 20px; overflow: auto; padding: 0px; text-align: left; width: 99%;"><code style="color: black; word-wrap: normal;">1: var content:Sprite = new Sprite()
2: var contentText:TextField = new TextField();
3: contentText.autoSize=TextFieldAutoSize.LEFT;
4: contentText.multiline = true;
5: contentText.htmlText = "INSERT YOUR INFO HERE!"
6: contentText.textColor = 0xFFFFFFFF;
7: content.graphics.clear();
8: content.graphics.beginFill(0xFF000000, 0);
9: content.graphics.drawRect(0, 0, 50, 45);
10: content.graphics.endFill();
11: content.addChild(contentText);
</code></pre><br />
Next, create a smart info window to place the content in:<br />
<pre style="background: url("https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEitTzQYpR3KnwgrimJHF2i-bLvCl-N0qx7_09T6Pzu3lZiCVXnlY1g4aZOMK-YRyKXOLGh8DBYDfUm_ed1IEqT3Sth-Xs-9L79Rl2FS7gKrXIGPNEXkhqFbJZ6n5g0PnLiZFHC3JMmuQqS7/s320/codebg.gif") repeat scroll 0% 0% rgb(240, 240, 240); border: 1px dashed rgb(204, 204, 204); color: black; font-family: arial; font-size: 12px; height: auto; line-height: 20px; overflow: auto; padding: 0px; text-align: left; width: 99%;"><code style="color: black; word-wrap: normal;">1: var infowindow:SmartInfoWindow = new SmartInfoWindow();
2: infowindow.position = infowindow.getBestAlignment(map, marker.getLatLng(), new Point(content.width, content.height));
3: infowindow.padding = 10;
4: infowindow.cornerRadius = 3;
5: infowindow.fillStyle = new FillStyle({color:0xFF000000, alpha:.8});
6: infowindow.strokeStyle = new StrokeStyle({color:0xFFC0C0C0,alpha:.8, thickness:1.5, pixelHinting:true});
7: infowindow.setContent(content);
8: infowindow.render();
</code></pre><br />
Almost there! Create info window options and assign it the custom info window created above:<br />
<pre style="background: url("https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEitTzQYpR3KnwgrimJHF2i-bLvCl-N0qx7_09T6Pzu3lZiCVXnlY1g4aZOMK-YRyKXOLGh8DBYDfUm_ed1IEqT3Sth-Xs-9L79Rl2FS7gKrXIGPNEXkhqFbJZ6n5g0PnLiZFHC3JMmuQqS7/s320/codebg.gif") repeat scroll 0% 0% rgb(240, 240, 240); border: 1px dashed rgb(204, 204, 204); color: black; font-family: arial; font-size: 12px; height: auto; line-height: 20px; overflow: auto; padding: 0px; text-align: left; width: 99%;"><code style="color: black; word-wrap: normal;">1: var options:InfoWindowOptions = new InfoWindowOptions();
2: options.customContent = infowindow;
3: options.customOffset = infowindow.anchorPoint;
4: options.hasShadow = false;
5: options.drawDefaultFrame = false;
6: options.hasCloseButton = true;
7: options.customCloseRect = infowindow.closeButtonRect;
8: options.pointOffset = new Point(0,-5);
</code></pre><br />
Lastly, create an event listener and add the marker to the map:<br />
<pre style="background: url("https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEitTzQYpR3KnwgrimJHF2i-bLvCl-N0qx7_09T6Pzu3lZiCVXnlY1g4aZOMK-YRyKXOLGh8DBYDfUm_ed1IEqT3Sth-Xs-9L79Rl2FS7gKrXIGPNEXkhqFbJZ6n5g0PnLiZFHC3JMmuQqS7/s320/codebg.gif") repeat scroll 0% 0% rgb(240, 240, 240); border: 1px dashed rgb(204, 204, 204); color: black; font-family: arial; font-size: 12px; height: auto; line-height: 20px; overflow: auto; padding: 0px; text-align: left; width: 99%;"><code style="color: black; word-wrap: normal;">1: marker.addEventListener(MapMouseEvent.CLICK, function (e:MapMouseEvent):void{map.openInfoWindow(e.target.getLatLng(),options);});
2: map.addOverlay(marker);
</code></pre><br />
There are a couple other features I think that bring this little app full circle. The roll over highlighting of the markers, the double Y axis in the chart, and overlay of the line chart onto the column chart. But I'll cover those in future posts. <br />
<div style="text-align: center;"><span style="font-size: large;"><b>HYDROGEN VIEWER</b></span></div><embed align="left" allowscriptaccess="always" height="640" pluginspage="http://www.macromedia.com/go/getflashplayer" quality="high" src="http://rpm.nrel.gov/GridView/Hydrogen/HydrogenViewer.swf" type="application/x-shockwave-flash" width="640"></embed><br />
<blockquote></blockquote>Anonymoushttp://www.blogger.com/profile/16259376360116803538noreply@blogger.com0tag:blogger.com,1999:blog-9041577666726226947.post-7913871813298914212010-10-16T13:13:00.001-06:002011-05-19T07:21:03.159-06:00Tween Google Marker Icon in FlexThere are some powerful and fast tweening engines for use in Adobe Flex. One in particular that I have come accustomed to is <a href="http://www.greensock.com/tweenlite/">TweenLite</a>. Recently I had been toying around with the idea of animating a Google marker icon by dynamically resizing its radius.<br />
<br />
But just how do you access the icon and tween its size? Well, tweening the icon is easy using TweenLite, just make sure you do a couple things first:<br />
<br />
<br />
<br />
1. Create a Sprite to use as a marker, give the sprite some color and a default size<br />
<span style="font-weight: bold;"><br />
</span><span style="font-weight: bold;"></span><br />
<pre style="background: url("https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEitTzQYpR3KnwgrimJHF2i-bLvCl-N0qx7_09T6Pzu3lZiCVXnlY1g4aZOMK-YRyKXOLGh8DBYDfUm_ed1IEqT3Sth-Xs-9L79Rl2FS7gKrXIGPNEXkhqFbJZ6n5g0PnLiZFHC3JMmuQqS7/s320/codebg.gif") repeat scroll 0% 0% rgb(240, 240, 240); border: 1px dashed rgb(204, 204, 204); color: black; font-family: arial; font-size: 12px; height: auto; line-height: 20px; overflow: auto; padding: 0px; text-align: left; width: 99%;"><span style="font-weight: bold;"><code style="color: black; word-wrap: normal;"><span style="font-weight: normal;">1: myIcon:Sprite = new Sprite();</span>
<span style="font-weight: normal;">2: myIcon.graphics.beginFill(0xC4D79B, 0.8);</span>
<span style="font-weight: normal;">3: myIcon.drawCircle(0, 0, 15);</span>
<span style="font-weight: normal;">4: myIcon.endFill();</span>
</code></span></pre><span style="font-weight: bold;"><br />
<br />
</span>2. Create marker options and set the icon as the sprite created above<br />
<pre style="background: url("https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEitTzQYpR3KnwgrimJHF2i-bLvCl-N0qx7_09T6Pzu3lZiCVXnlY1g4aZOMK-YRyKXOLGh8DBYDfUm_ed1IEqT3Sth-Xs-9L79Rl2FS7gKrXIGPNEXkhqFbJZ6n5g0PnLiZFHC3JMmuQqS7/s320/codebg.gif") repeat scroll 0% 0% rgb(240, 240, 240); border: 1px dashed rgb(204, 204, 204); color: black; font-family: arial; font-size: 12px; height: auto; line-height: 20px; overflow: auto; padding: 0px; text-align: left; width: 99%;"><code style="color: black; word-wrap: normal;">1: markerOptions:MarkerOptions = new MarkerOptions();
2: markerOptions.icon = myIcon;
3: marker = new Marker(new LatLng(39.75, -104.87), markerOptions);
4: map.addOverlay(marker);
</code></pre><br />
<br />
3. Create some function that houses and executes your tweening effect. The tween itself is:<br />
<pre style="background: url("https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEitTzQYpR3KnwgrimJHF2i-bLvCl-N0qx7_09T6Pzu3lZiCVXnlY1g4aZOMK-YRyKXOLGh8DBYDfUm_ed1IEqT3Sth-Xs-9L79Rl2FS7gKrXIGPNEXkhqFbJZ6n5g0PnLiZFHC3JMmuQqS7/s320/codebg.gif") repeat scroll 0% 0% rgb(240, 240, 240); border: 1px dashed rgb(204, 204, 204); color: black; font-family: arial; font-size: 12px; height: auto; line-height: 20px; overflow: auto; padding: 0px; text-align: left; width: 99%;"><code style="color: black; word-wrap: normal;">1: TweenLite.to(marker.getOptions().icon, 2, {scaleX: 3, scaleY: 3});
</code></pre><br />
That's it! Just make sure if you are going to be tweening multiple markers you have some way of storing and recalling them as needed.<br />
<br />
<embed align="middle" allowscriptaccess="always" height="450" pluginspage="http://www.macromedia.com/go/getflashplayer" quality="high" src="http://rpm.nrel.gov/GridView/GoogleMarkerExample/GM.swf" type="application/x-shockwave-flash" width="500"></embed>Anonymoushttp://www.blogger.com/profile/16259376360116803538noreply@blogger.com0