tag:blogger.com,1999:blog-86336234746653144662024-03-05T18:35:35.018+01:00webrianTechnical GIS snippets and related noteswebrianhttp://www.blogger.com/profile/15325993494848772844noreply@blogger.comBlogger35125tag:blogger.com,1999:blog-8633623474665314466.post-87494154683835912872022-12-18T21:51:00.004+01:002022-12-18T21:52:33.086+01:00Unzip Processing Algorithm for QGIS<p> While there is the handy QGIS processing algorithm <a href="https://docs.qgis.org/3.22/en/docs/user_manual/processing_algs/qgis/filetools.html#download-file" target="_blank">Download file</a>, often publicly availabe (geo-)data are compressed in a ZIP archive. That's why I've started to develop a "Unzip Archive" processing algorithm for use in the QGIS processing toolbox. Thanks to the already available class <a href="https://qgis.org/pyqgis/3.28/core/QgsZipUtils.html#module-QgsZipUtils" target="_blank">QgsZipUtils</a> it was not that difficult.</p><div class="separator" style="clear: both; text-align: center;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEiEY_EEg-SzwGbcbJgclS7BsRHckIT9NfIIylSDFbQ1Ic50IU1ajNly1xUvOi0I8vUtDrxhyYS-rH8755mXZqG2Yw2veJeqGCR_kvS1DI0biIyJYQoy7VZR7PumvnfA-2Q6Ke5wU6YmypFTy3XEq-XLfI9ivQ4xi7n-OVDD8eJm3VTCwIbqBHWWlv4Z/s597/model.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" data-original-height="435" data-original-width="597" height="233" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEiEY_EEg-SzwGbcbJgclS7BsRHckIT9NfIIylSDFbQ1Ic50IU1ajNly1xUvOi0I8vUtDrxhyYS-rH8755mXZqG2Yw2veJeqGCR_kvS1DI0biIyJYQoy7VZR7PumvnfA-2Q6Ke5wU6YmypFTy3XEq-XLfI9ivQ4xi7n-OVDD8eJm3VTCwIbqBHWWlv4Z/s320/model.png" width="320" /></a></div><p></p><span><a name='more'></a></span><p>In combination with the Download File algorithm it is quite easy to create a model that downloads and unzip data.<br /></p><p>The algorithm is published on Github.</p>webrianhttp://www.blogger.com/profile/15325993494848772844noreply@blogger.com0tag:blogger.com,1999:blog-8633623474665314466.post-12526466600944657862021-12-13T20:44:00.000+01:002021-12-13T20:44:13.129+01:00Updates in 2021<p>Apart from updating and rerun the scripts for the cloudless atlas, unfortunately no updates ... :(</p>webrianhttp://www.blogger.com/profile/15325993494848772844noreply@blogger.com0tag:blogger.com,1999:blog-8633623474665314466.post-73776999280818096172020-12-22T22:58:00.003+01:002020-12-22T23:19:44.032+01:00Cloudless Atlas Laos 2020<p style="text-align: left;">To track down new hydro power reservoir and schemes I've downloaded most 2020 Landsat 8 scenes of northern Laos and applied my <a href="/2014/07/imitating-mapbox-cloudless-atlas.html">cloudless script</a>.<br /></p><table align="center" cellpadding="0" cellspacing="0" class="tr-caption-container" style="margin-left: auto; margin-right: auto;"><tbody><tr><td style="text-align: center;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEiJzF_hSmiha5kfEu640LtAPwXQPBjmgbqyvGqKP144rumL1_Wxp1GSk12Te1lrdrJ0Tzi6SUq-mdSTI8DX2YqUyT30Dv_xmba6sYpUMwBkWBu0iStmaQLJ3JFCeAsRwVL8hYwRdF_bjyQ/s1496/cloudlessatlas2020_laos.jpeg" style="margin-left: auto; margin-right: auto;"><img alt="Cloudless Atlas Laos 2020" border="0" data-original-height="892" data-original-width="1496" height="239" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEiJzF_hSmiha5kfEu640LtAPwXQPBjmgbqyvGqKP144rumL1_Wxp1GSk12Te1lrdrJ0Tzi6SUq-mdSTI8DX2YqUyT30Dv_xmba6sYpUMwBkWBu0iStmaQLJ3JFCeAsRwVL8hYwRdF_bjyQ/w400-h239/cloudlessatlas2020_laos.jpeg" title="Cloudless Atlas Laos 2020" width="400" /></a></td></tr><tr><td class="tr-caption" style="text-align: center;">Northern Laos</td><td class="tr-caption" style="text-align: center;"><br /></td></tr></tbody></table><p></p><p><br /></p><span><a name='more'></a></span><p>And indeed:</p><div class="separator" style="clear: both; text-align: center;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhMk7-l-AbOOUpV-U4P-Z777O0es2HG1mz3spwNkVwpeS_jqzk3fcicGiBzgRht1n2PFHEvM746sxLjmXooTjuqTidrhdXlDoXix9b5hiDCVXWMVp8d1HxlCjNWWbRoFR-FMBuvF3Dd960/s1496/hydropower_res.jpeg" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" data-original-height="892" data-original-width="1496" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhMk7-l-AbOOUpV-U4P-Z777O0es2HG1mz3spwNkVwpeS_jqzk3fcicGiBzgRht1n2PFHEvM746sxLjmXooTjuqTidrhdXlDoXix9b5hiDCVXWMVp8d1HxlCjNWWbRoFR-FMBuvF3Dd960/s320/hydropower_res.jpeg" width="320" /></a></div><br /><p><br /></p>webrianhttp://www.blogger.com/profile/15325993494848772844noreply@blogger.com0Vientiane Province, Laos18.5705063 102.621621116.49473774498027 100.424355475 20.646274855019733 104.818886725tag:blogger.com,1999:blog-8633623474665314466.post-14957736754498027722019-12-20T22:07:00.000+01:002020-12-22T22:59:22.670+01:00The Nam Ngiep Hydropower SchemeJust recently I continued mapping the Nam Ngiep hydropower scheme, see also my previous <a href="https://www.webrian.ch/2018/12/yet-another-hydropower-dam-nam-ngiep-1.html">post about Nam Ngiep 1</a>.<br />
<br />
So far my findings:<br />
<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/AVvXsEjMyLFwJuboUenZjqNzfnxC2z90QqVQ3xpg37iL08nM6Yr0U4TMag2ozMGHsdslHxmbkCluIemncehmMB1QORNCk03T5dFH38d1i_Q1x6jGb_21cdcWl0wH8sDX3RM-Z9OIGVYuqTq9zhk/s1600/NamNgiep1.png" imageanchor="1" style="margin-left: auto; margin-right: auto;"><img border="0" data-original-height="781" data-original-width="688" height="400" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjMyLFwJuboUenZjqNzfnxC2z90QqVQ3xpg37iL08nM6Yr0U4TMag2ozMGHsdslHxmbkCluIemncehmMB1QORNCk03T5dFH38d1i_Q1x6jGb_21cdcWl0wH8sDX3RM-Z9OIGVYuqTq9zhk/s400/NamNgiep1.png" width="351" /> </a></td><td style="text-align: center;"><br /></td><td style="text-align: center;"><br /></td><td style="text-align: center;"><br /></td><td style="text-align: center;"></td></tr>
<tr><td class="tr-caption" style="text-align: center;">Nam Ngiep1</td><td class="tr-caption" style="text-align: center;"><br /></td></tr>
</tbody></table>
<div style="text-align: left;">
<br />
<a name='more'></a><br /><br />
Nam Ngiep 2 which is actually on Nam Sen: <br />
<table cellpadding="0" cellspacing="0" class="tr-caption-container" style="margin-left: 0px; margin-right: auto; text-align: left;"><tbody></tbody></table>
</div>
<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/AVvXsEiCah4M0s0VzoC-NFhJmTTZyM0pbx_9htqOabecgSLfIyaLY0LyHPGpBapPzxg-uWafVuzEHjYIuKBOc1MQ8kBWmKEB0-ts5h9CGuojqBqcJWQIFn1xneXQ8ydOiVSLkr2pIZ2N9Z5aKys/s1600/NamNgiep2.png" imageanchor="1" style="margin-left: auto; margin-right: auto;"><img border="0" data-original-height="732" data-original-width="1289" height="226" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEiCah4M0s0VzoC-NFhJmTTZyM0pbx_9htqOabecgSLfIyaLY0LyHPGpBapPzxg-uWafVuzEHjYIuKBOc1MQ8kBWmKEB0-ts5h9CGuojqBqcJWQIFn1xneXQ8ydOiVSLkr2pIZ2N9Z5aKys/s400/NamNgiep2.png" width="400" /></a></td></tr>
<tr><td class="tr-caption" style="text-align: center;"> Nam Ngiep 2 </td></tr>
</tbody></table>
<br />
Nam Ngiep 3A further north with the powerhouse in the lower left corner:<br />
<br />
<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/AVvXsEj3wJESnXtabnv264NWLRUNJJlLmDnbRtgPdJVX3WtYvYTkb-h6ENGBwJZT19NTgCiqNdC1DY1ukG5b68UMPad7cal18AXHXpdLLLPx-NuPp8urRpsBTTNU0g6pCfMS9CcUXjNcJTD4Zdw/s1600/NamNgiep3a.png" imageanchor="1" style="margin-left: auto; margin-right: auto;"><img border="0" data-original-height="828" data-original-width="921" height="358" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEj3wJESnXtabnv264NWLRUNJJlLmDnbRtgPdJVX3WtYvYTkb-h6ENGBwJZT19NTgCiqNdC1DY1ukG5b68UMPad7cal18AXHXpdLLLPx-NuPp8urRpsBTTNU0g6pCfMS9CcUXjNcJTD4Zdw/s400/NamNgiep3a.png" width="400" /></a></td></tr>
<tr><td class="tr-caption" style="text-align: center;">Nam Ngiep 3</td></tr>
</tbody></table>
<br />
More dams to come with <a href="https://siteresources.worldbank.org/INTLAOPRD/513906-1095113468080/20252949/Section%202%20Project%20Profiles%202_16%20to%202_19.pdf">the Nam Pot dam</a>:<br />
<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/AVvXsEi-34pljUGdqA4cn5vtjwREcEw_ZSc28Rri7rBAXPUPxSOn1nR4a4rNy9C-anNKpmI2pMfKeDlzQbvLrMNthL_TcvhQpGRNPseWZpcd-UfgL_Z8LoMgaJnifGP8nUDorUIax6cZhX7cpmU/s1600/NamPot.png" imageanchor="1" style="margin-left: auto; margin-right: auto;"><img border="0" data-original-height="1053" data-original-width="930" height="400" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEi-34pljUGdqA4cn5vtjwREcEw_ZSc28Rri7rBAXPUPxSOn1nR4a4rNy9C-anNKpmI2pMfKeDlzQbvLrMNthL_TcvhQpGRNPseWZpcd-UfgL_Z8LoMgaJnifGP8nUDorUIax6cZhX7cpmU/s400/NamPot.png" width="352" /></a></td></tr>
<tr><td class="tr-caption" style="text-align: center;">Nam Pot dam under contstruction</td></tr>
</tbody></table>
<br />
<br />webrianhttp://www.blogger.com/profile/15325993494848772844noreply@blogger.com0Unnamed Road, Laos18.909800472463541 103.6421558575525518.669373472463541 103.31943235755256 19.150227472463541 103.96487935755255tag:blogger.com,1999:blog-8633623474665314466.post-66210699603749459522018-12-07T19:30:00.000+01:002018-12-07T22:35:10.614+01:00Yet Another Hydropower Dam: Nam Ngiep 1Meanwhile <a href="https://namngiep1.com/" target="_blank">Nam Ngiep 1 dam</a> has been completed and the reservoir has been filled during rainy season in 2018. I.e., it's high time to put in <a href="https://www.openstreetmap.org/#map=11/18.8394/103.4359" target="_blank">on the map</a>!<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/AVvXsEg-fbqq4xXShtU3hfOnpSEvlTFg7SxnqZBsy5HddyAOPj9KFHw8dVUBpI3ZFEFv0UzyqcVpIpAMKnGjII0I8HVWBo9hz8VIiufsDtNemdpMv_jL9kthb2gn8s_fsY0jllP95_kQnfxchH0/s1600/namngiep1_jan18.png" imageanchor="1" style="margin-left: auto; margin-right: auto;"><img border="0" data-original-height="837" data-original-width="837" height="320" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEg-fbqq4xXShtU3hfOnpSEvlTFg7SxnqZBsy5HddyAOPj9KFHw8dVUBpI3ZFEFv0UzyqcVpIpAMKnGjII0I8HVWBo9hz8VIiufsDtNemdpMv_jL9kthb2gn8s_fsY0jllP95_kQnfxchH0/s320/namngiep1_jan18.png" width="320" /></a></td></tr>
<tr><td class="tr-caption" style="text-align: center;">Nam Ngiep 1 reservoir area in January 2018</td></tr>
</tbody></table>
<br />
<a name='more'></a><br />
<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/AVvXsEh7Q9c8vx2IIORruN1Lv4xL3rDlxbWGbHtfd-GSdgCdQBR3TkWSr6R3sWpP7pFtLECEkksVTNzIsI2Thju6msrwOP1r5DqoIfwBFFeLJNO8tysG9h82zcAV7JMmNVtkQ3T8KKJS-rpVmOI/s1600/namngiep1_okt18.png" imageanchor="1" style="margin-left: auto; margin-right: auto;"><img border="0" data-original-height="838" data-original-width="838" height="320" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEh7Q9c8vx2IIORruN1Lv4xL3rDlxbWGbHtfd-GSdgCdQBR3TkWSr6R3sWpP7pFtLECEkksVTNzIsI2Thju6msrwOP1r5DqoIfwBFFeLJNO8tysG9h82zcAV7JMmNVtkQ3T8KKJS-rpVmOI/s320/namngiep1_okt18.png" width="320" /></a></td></tr>
<tr><td class="tr-caption" style="text-align: center;">Nam Ngiep 1 reservoir in October 2018</td></tr>
</tbody></table>
<span id="goog_1244604741"></span><span id="goog_1244604742"></span>
<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/AVvXsEhmaVeau75BuBsYw9hkckhakd8WieIS4G39mh7Tmh6nQ0dwqxaTMTXiqi-5mYwEVtT9IZGjlw1c_YMskfmfwsUEjc-IhaGzZnqNv4eSdrdfBVmAg0z0La5przkRtkYK9QISFa4Vo7-1WCY/s1600/namngiep1_gm_overlay.jpg" imageanchor="1" style="margin-left: auto; margin-right: auto;"><img border="0" data-original-height="838" data-original-width="838" height="320" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhmaVeau75BuBsYw9hkckhakd8WieIS4G39mh7Tmh6nQ0dwqxaTMTXiqi-5mYwEVtT9IZGjlw1c_YMskfmfwsUEjc-IhaGzZnqNv4eSdrdfBVmAg0z0La5przkRtkYK9QISFa4Vo7-1WCY/s320/namngiep1_gm_overlay.jpg" width="320" /></a></td></tr>
<tr><td class="tr-caption" style="text-align: center;">Goodbye villages Na Nhao and Vang Naxay!</td></tr>
</tbody></table>
<br />webrianhttp://www.blogger.com/profile/15325993494848772844noreply@blogger.com0Vientiane Province, Laos18.714639877380719 103.4361948839787118.65448687738072 103.35551388397872 18.774792877380719 103.51687588397871tag:blogger.com,1999:blog-8633623474665314466.post-3646743942806217452017-11-01T21:05:00.000+01:002017-11-01T21:05:40.322+01:00New hydropower dam: Nam PhayJust recently I've put the Nam Phay hydropower reservoir <a href="https://www.openstreetmap.org/#map=13/19.1182/102.7116" target="_blank">on the map</a>!<br />
<br />
During Lao New Year 2014 I had the chance to visite this area, see pictures below:<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/AVvXsEi16X573knqFrWWseNRbE9bW8fVb81bHFZpipzxxLFmcl-vuC1in1K3sTnrNptE9PJgUIVxkiRWSN2c2Ao2vQFa3I-489KBc7OTBSRdKpslg7l5uAel2aI2Z_OWcSYi6u9V902C5a45LNw/s1600/NamPhay_20150113.jpg" imageanchor="1" style="margin-left: auto; margin-right: auto;"><img border="0" data-original-height="838" data-original-width="838" height="400" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEi16X573knqFrWWseNRbE9bW8fVb81bHFZpipzxxLFmcl-vuC1in1K3sTnrNptE9PJgUIVxkiRWSN2c2Ao2vQFa3I-489KBc7OTBSRdKpslg7l5uAel2aI2Z_OWcSYi6u9V902C5a45LNw/s400/NamPhay_20150113.jpg" width="400" /></a></td></tr>
<tr><td class="tr-caption" style="text-align: center;">The Nam Phay reservoir area in Januar 2015</td><td class="tr-caption" style="text-align: center;"><br /></td><td class="tr-caption" style="text-align: center;"><br /></td></tr>
</tbody></table>
<br />
<a name='more'></a><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/AVvXsEgjQLYRe8hYpDh7Ke7bjR8yac64HTeVr6I27l5OAmTLbXvpDGtH0cmvJV65mJlAfr3-Y4rN6HnJQfrdGeGRtM2E4sYfFZEXEz2z96DTldOcyNGZQ5AykwKe5FTEemGNQ5W7A_07z_6iBn0/s1600/NamPhay_20170329.jpg" imageanchor="1" style="margin-left: auto; margin-right: auto;"><img border="0" data-original-height="838" data-original-width="838" height="400" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgjQLYRe8hYpDh7Ke7bjR8yac64HTeVr6I27l5OAmTLbXvpDGtH0cmvJV65mJlAfr3-Y4rN6HnJQfrdGeGRtM2E4sYfFZEXEz2z96DTldOcyNGZQ5AykwKe5FTEemGNQ5W7A_07z_6iBn0/s400/NamPhay_20170329.jpg" width="400" /></a></td></tr>
<tr><td class="tr-caption" style="text-align: center;">The Nam Phay reservoir in March 2017</td></tr>
</tbody></table>
<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/AVvXsEhJpEl2Js3A-lGn9ao4eodJwLKQuWnzM_AFyJNir8vrRy75ZdX9B0RPOl85TJtaQw-hNSBoB53W9WRooEYToyD1DUtCko9O6rHc8Mu8L8OqwuOLrhfPyyY-9woHkZi5CIbz3ungt22i5Po/s1600/namphay.jpg" imageanchor="1" style="margin-left: auto; margin-right: auto;"><img border="0" data-original-height="1045" data-original-width="1461" height="285" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhJpEl2Js3A-lGn9ao4eodJwLKQuWnzM_AFyJNir8vrRy75ZdX9B0RPOl85TJtaQw-hNSBoB53W9WRooEYToyD1DUtCko9O6rHc8Mu8L8OqwuOLrhfPyyY-9woHkZi5CIbz3ungt22i5Po/s400/namphay.jpg" width="400" /></a></td></tr>
<tr><td class="tr-caption" style="text-align: center;">An overview map found on the website of the Department of Energy Business</td></tr>
</tbody></table>
<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/AVvXsEhJhrqdc_NM_mub7kmBCCFogDx5oQnzFydWpBRdbE815e_pnRP-lmAFtEaLsWbxZhsNC63M0P53OkbNg9DQXOAF0DbC6eUySGNwet8BvhMekkBXe1qSAM-hQ8OXD-E0Da1opzSd8tBA7ps/s1600/IMG_3407.JPG" imageanchor="1" style="margin-left: auto; margin-right: auto;"><img border="0" data-original-height="684" data-original-width="912" height="300" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhJhrqdc_NM_mub7kmBCCFogDx5oQnzFydWpBRdbE815e_pnRP-lmAFtEaLsWbxZhsNC63M0P53OkbNg9DQXOAF0DbC6eUySGNwet8BvhMekkBXe1qSAM-hQ8OXD-E0Da1opzSd8tBA7ps/s400/IMG_3407.JPG" width="400" /></a></td></tr>
<tr><td class="tr-caption" style="text-align: center;">On the way to the reservoir construction site</td></tr>
</tbody></table>
<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/AVvXsEgO7dNtCfawpuF3PkO47VPz5uDNn_HiSJrTrqBkeuw5mbtOafL2mu02QlPNc-5Zk3lmH6kWPbfv9lAf0z-gYIs9BdChjRsWPkMutYfvOeONqXz-CG2Wno9Lkk9_UHGOmPlOvzFG1kPyVTY/s1600/village1.jpg" imageanchor="1" style="margin-left: auto; margin-right: auto;"><img border="0" data-original-height="684" data-original-width="912" height="300" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgO7dNtCfawpuF3PkO47VPz5uDNn_HiSJrTrqBkeuw5mbtOafL2mu02QlPNc-5Zk3lmH6kWPbfv9lAf0z-gYIs9BdChjRsWPkMutYfvOeONqXz-CG2Wno9Lkk9_UHGOmPlOvzFG1kPyVTY/s400/village1.jpg" width="400" /></a></td></tr>
<tr><td class="tr-caption" style="text-align: center;">One of the village which has been flooded by the reservoir</td></tr>
</tbody></table>
<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/AVvXsEj9ai5sbltY9s99eVUS3YrTPdTMR2vmGC_LL6niurEnSwI_9Hq-VPIPvrbYJqTnTqetEPXl1nAO9LT48txZZMUWWn1F9gMJJHRqaoMxtM8guROm7FsJqswPrvay1eWPfuGO-DfUiFwPEJg/s1600/waterintake1.jpg" imageanchor="1" style="margin-left: auto; margin-right: auto;"><img border="0" data-original-height="684" data-original-width="912" height="300" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEj9ai5sbltY9s99eVUS3YrTPdTMR2vmGC_LL6niurEnSwI_9Hq-VPIPvrbYJqTnTqetEPXl1nAO9LT48txZZMUWWn1F9gMJJHRqaoMxtM8guROm7FsJqswPrvay1eWPfuGO-DfUiFwPEJg/s400/waterintake1.jpg" width="400" /></a></td></tr>
<tr><td class="tr-caption" style="text-align: center;">The water intake</td></tr>
</tbody></table>
<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/AVvXsEj2ejhgoRJg0Ne7zFgzHt9vx-qYBOqQnRr0LD9ojJ4kJahqbMu4qssTlsHIF9ybqNx-0ti2-JyPIGD3U-ImM1-6Frq84H3JjP7C-2Sndh08_7iKuZNqLtMh0WNanlMamVgJ3c3XX5mHmu8/s1600/waterintake2.jpg" imageanchor="1" style="margin-left: auto; margin-right: auto;"><img border="0" data-original-height="684" data-original-width="912" height="300" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEj2ejhgoRJg0Ne7zFgzHt9vx-qYBOqQnRr0LD9ojJ4kJahqbMu4qssTlsHIF9ybqNx-0ti2-JyPIGD3U-ImM1-6Frq84H3JjP7C-2Sndh08_7iKuZNqLtMh0WNanlMamVgJ3c3XX5mHmu8/s400/waterintake2.jpg" width="400" /></a></td></tr>
<tr><td class="tr-caption" style="text-align: center;">The water intake close-up</td></tr>
</tbody></table>
<br />
<br />webrianhttp://www.blogger.com/profile/15325993494848772844noreply@blogger.comUnnamed Road, Laos18.924474330741365 102.782592773437518.683994330741363 102.4598692734375 19.164954330741367 103.1053162734375tag:blogger.com,1999:blog-8633623474665314466.post-63561565330786159692016-11-26T22:36:00.002+01:002016-11-26T22:36:39.682+01:00Another Hydropower dam: Xekaman 1<div class="separator" style="clear: both; text-align: left;">
Just recently I've found by chance the Xekaman 1 hydropower reservoir on <a href="http://landsat.usgs.gov/landsat8.php">Landsat 8 imagery</a>. Time to put it on <a href="http://www.openstreetmap.org/">the map</a>!</div>
<div class="separator" style="clear: both; text-align: left;">
<br /></div>
<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/AVvXsEhfsPzwqdOmdtrk4gfOY3kMFnpXGeB88CqFQQ3kDnKvm6vlpad4E8Fs1zrhDkdMcIDGvhMdMH6MnaRWye1wiukNILLTr9Th6nGf84ZknrKuzinTWIYYq3twdRA6Gazt5OdErbxvDPlpvvw/s1600/2015-01-24_xekaman1.png" imageanchor="1" style="margin-left: auto; margin-right: auto;"><img border="0" height="208" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhfsPzwqdOmdtrk4gfOY3kMFnpXGeB88CqFQQ3kDnKvm6vlpad4E8Fs1zrhDkdMcIDGvhMdMH6MnaRWye1wiukNILLTr9Th6nGf84ZknrKuzinTWIYYq3twdRA6Gazt5OdErbxvDPlpvvw/s320/2015-01-24_xekaman1.png" width="320" /></a></td></tr>
<tr><td class="tr-caption" style="text-align: center;">The Xekaman river area recorded on 24th January 2015, just before water storage started</td></tr>
</tbody></table>
<a name='more'></a><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/AVvXsEhCDIniaYD4yxB1AUpOV5eWTUeemp6p9pHYh8KAwaYpPDFCU6b9d9NmbVcL42f5o3_FXGcwhy-jmVgfOlTKfXkk0mxPFFCjtt4T9d1gVzdz0n_3MB4zERBtsFzdiPpu0j80VgcOLBQKN9s/s1600/2016-06-19_xekaman1.png" imageanchor="1" style="margin-left: auto; margin-right: auto;"><img border="0" height="209" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhCDIniaYD4yxB1AUpOV5eWTUeemp6p9pHYh8KAwaYpPDFCU6b9d9NmbVcL42f5o3_FXGcwhy-jmVgfOlTKfXkk0mxPFFCjtt4T9d1gVzdz0n_3MB4zERBtsFzdiPpu0j80VgcOLBQKN9s/s320/2016-06-19_xekaman1.png" width="320" /></a></td></tr>
<tr><td class="tr-caption" style="text-align: center;">The exactly same area on 19th June 2016</td></tr>
</tbody></table>
<br /><div>
I've used the well-proven workflow in <a href="https://grass.osgeo.org/">GRASS GIS</a> as described in previous posts.<div>
<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/AVvXsEiA7PWAjVFjHV04Qlr6C_8sUj-8tRjhp7zMnIOlDYF-T8vx1oI-zOqp700xr5JT72bRgw3z2rWD9hgGFAClAJjOmMQGNDNXw2uY8FjVovx_TOYGbg4ZGDcj1KrCR5QOg02REcwv5lROz9k/s1600/xekaman1_ndvi.png" imageanchor="1" style="margin-left: auto; margin-right: auto;"><img border="0" height="320" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEiA7PWAjVFjHV04Qlr6C_8sUj-8tRjhp7zMnIOlDYF-T8vx1oI-zOqp700xr5JT72bRgw3z2rWD9hgGFAClAJjOmMQGNDNXw2uY8FjVovx_TOYGbg4ZGDcj1KrCR5QOg02REcwv5lROz9k/s320/xekaman1_ndvi.png" width="281" /></a></td></tr>
<tr><td class="tr-caption" style="text-align: center;">NDVI map of Xekaman 1 reservoir</td></tr>
</tbody></table>
<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/AVvXsEjxdM9FM4N1lf_XWFqLAyDpozS8wB5zapiARzhVRod23WUXQx2xyu2nDCxLQBiDVlJWoG6-8XIinDVaZxdq75QmVDVfBAPZAgAmPFx80tpnpo4cWD_fJLjdJrKWeCcwWPxxUdJGWSGTCaI/s1600/xekaman1_final.png" imageanchor="1" style="margin-left: auto; margin-right: auto;"><img border="0" height="320" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjxdM9FM4N1lf_XWFqLAyDpozS8wB5zapiARzhVRod23WUXQx2xyu2nDCxLQBiDVlJWoG6-8XIinDVaZxdq75QmVDVfBAPZAgAmPFx80tpnpo4cWD_fJLjdJrKWeCcwWPxxUdJGWSGTCaI/s320/xekaman1_final.png" width="281" /></a></td></tr>
<tr><td class="tr-caption" style="text-align: center;">The final result</td></tr>
</tbody></table>
</div>
</div>
<div>
Let's find more missing reservoirs to put those on the map!</div>
webrianhttp://www.blogger.com/profile/15325993494848772844noreply@blogger.com0Attapeu Province, Laos15.095335677611214 107.2712269531250514.604745177611214 106.62577995312505 15.585926177611213 107.91667395312506tag:blogger.com,1999:blog-8633623474665314466.post-22079165209853987382015-11-29T15:23:00.000+01:002015-11-30T16:48:39.549+01:00Raster math in GRASS GIS using numpyHouay Lamphan Gnai, yet another hydropower plant just started commercial operations recently
with a reservoir area of seven square kilometers as reported by local media. Time to put it
on OpenStreetMap using up-to-date Landsat 8 imagery and the new
<a href="http://www.gissula.eu/qgis-grass-plugin-crowdfunding/">GRASS GIS 7 plugin</a> for
<a href="http://qgis.org/">QGIS 2.12</a>. To get familiar with GRASS GIS'
<a href="https://grass.osgeo.org/grass70/manuals/libpython/script.html#module-script.array">interface to numpy</a>,
I've decided to extract the reservoir using numpy.
<br /><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/AVvXsEj55P5t0BfShqE8uurjEWybPj35AxJI3vDBBTOnaiZCKZqizFUNilMM4mUuMoUz5XLkNtpw4r1Jdsl7szBEGZsOEezu2SPCZdoYzdW9PipWclRIp9I-PCCtP3zzOcM4UOERueA7snVNn-c/s1600/houay-lamphan-landsat8_b5.jpg" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img alt="Houay Lamphan Gnai on Landsat 8 image" border="0" height="400" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEj55P5t0BfShqE8uurjEWybPj35AxJI3vDBBTOnaiZCKZqizFUNilMM4mUuMoUz5XLkNtpw4r1Jdsl7szBEGZsOEezu2SPCZdoYzdW9PipWclRIp9I-PCCtP3zzOcM4UOERueA7snVNn-c/s400/houay-lamphan-landsat8_b5.jpg" title="" width="372" /></a></td></tr>
<tr><td class="tr-caption" style="text-align: center;">Houay Lamphan Gnai on Landsat 8 image taken on the 23rd Oct, 2015</td></tr>
</tbody></table>
<br />
<a name='more'></a><br />
As with Landsat 7 I used the <a href="http://landsat.usgs.gov/band_designations_landsat_satellites.php">near infrared (NIR) channel of Landsat 8</a> with a wavelength of 0.85 - 0.88 micrometers. Whereas the infrared channel was band 4 in Landsat 7, it is now band 5 in Landsat 8.
<br />
The following script extracts pixels representing waterbodies (low values) and creates a new raster
layer with 1.0 and null values. It takes as arguments the name of the input map, name of the output map and
a threshold value.
<br />
<pre class="brush: python">#!/usr/bin/env python
from grass.script import array as garray
import numpy as np
import sys
def main(argv=None):
if argv is None:
argv = sys.argv
inputname = argv[1]
outputname = argv[2]
threshold = np.float_(argv[3])
in_arr = garray.array()
in_arr.read(inputname)
out_arr = garray.array()
out_arr[...] = in_arr
out_arr[out_arr > threshold] = np.nan
out_arr[out_arr <= threshold] = np.float_(1.0)
out_arr.write(outputname, overwrite=True)
return 0
if __name__ == '__main__':
sys.exit(main())
</pre>
After extracting the reservoir as raster layer, I continued with the <a href="http://www.webrian.ch/2012/02/armchair-mapping-using-grass-gis.html">previously described workflow</a>.
<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/AVvXsEg_3wMYQtKswj8FmXWxaGzXnfixxa4XAC1sc3bseZcW8EaOdGqLWZghVH5f71zL9-JirN7fVVbLg1Xl1JIeEFswxz_12iqi0lkHJf_yNScjzNi45Eje572USl1S5nt3IqPHchyphenhyphenlir7LsyE/s1600/houay-lamphan-vector.jpg" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img alt="Houay Lamphan Gnai smoothed vector" border="0" height="400" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEg_3wMYQtKswj8FmXWxaGzXnfixxa4XAC1sc3bseZcW8EaOdGqLWZghVH5f71zL9-JirN7fVVbLg1Xl1JIeEFswxz_12iqi0lkHJf_yNScjzNi45Eje572USl1S5nt3IqPHchyphenhyphenlir7LsyE/s400/houay-lamphan-vector.jpg" title="" width="372" /></a></td></tr>
<tr><td class="tr-caption" style="text-align: center;">Final result after extracting pixels representing water bodies, growing by one pixel, vectorizing and smoothing</td></tr>
</tbody></table>
<br />
I was wondering about the performance of this script against
GRASS GIS' own <a href="https://grass.osgeo.org/grass70/manuals/r.mapcalc.html">map calculator</a>. So
I compared above script with this one:
<br />
<pre class="brush: python">#!/usr/bin/env python
from grass.script.raster import mapcalc
import sys
def main(argv=None):
if argv is None:
argv = sys.argv
inputname = argv[1]
outputname = argv[2]
threshold = float(argv[3])
expr = '"%s" = if("%s" <= %f, 1, null())' % (outputname,inputname,threshold)
mapcalc(expr, overwrite=True)
return 0
if __name__ == '__main__':
sys.exit(main())
</pre>
For a region with 320223 cells I got the following results:
<br />
<pre class="brush: bash;">script with numpy: 0.357898100217
script with r.mapcalc: 0.270698928833
</pre>
As actually expected r.mapcalc performs better in this simple case. But still the
interface to numpy unlocks the great potential and functionality of numpy which excels
GRASS GIS' map calculator.
<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/AVvXsEhIhye2y39jY7Z-sWipS6q5hF3jkQ_IDf9baWFNr7Id6k2lzGD8IM0gCc9kPiN4fmJjm1OEj52ZFn-0QnWGU-WecgNjX26lnhyphenhyphenhOMnsZI6qs6jQazmUn0dG_R08ZINUKSb4ZFuqCaW2GXk/s1600/houay-lamphan-josm.jpg" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img alt="Editing newly created reservoir in JOSM" border="0" height="275" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhIhye2y39jY7Z-sWipS6q5hF3jkQ_IDf9baWFNr7Id6k2lzGD8IM0gCc9kPiN4fmJjm1OEj52ZFn-0QnWGU-WecgNjX26lnhyphenhyphenhOMnsZI6qs6jQazmUn0dG_R08ZINUKSb4ZFuqCaW2GXk/s400/houay-lamphan-josm.jpg" title="" width="400" /></a></td></tr>
<tr><td class="tr-caption" style="text-align: center;">Editing the newly created reservoir in JOSM</td></tr>
</tbody></table>
<br />
Is there any other map or dataset than OpenStreetMap that contains this reservoir right now? I don't think so!webrianhttp://www.blogger.com/profile/15325993494848772844noreply@blogger.com0Salavan, Laos15.424558071072397 106.364135742187515.179804571072397 106.0414122421875 15.669311571072397 106.6868592421875tag:blogger.com,1999:blog-8633623474665314466.post-49336496415803127172014-08-14T15:14:00.000+02:002014-08-15T17:08:08.505+02:00Web Map Service with PyCairo and PyramidFrom time to time I have the chance to work with <a href="http://cairographics.org/documentation/pycairo/2/">PyCairo</a>, the Python bindings to the <a href="http://cairographics.org/">cairo</a> graphic library. I appreciate this library because it allows to draw graphics primitives so super fast and thus it is very suitable to use it in web projects.<br />
After my trials with <a href="http://www.webrian.ch/2013/12/grass-gis-web-map-service-with-pyramid.html">GRASS GIS as WMS backend</a> I thought I could create a raster web map service with improved performance by using PyCairo instead of GRASS GIS as renderer.<br />
<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/AVvXsEgANSBhR1Of7968JuGpvmrMFw1j3_UnRqg-IUcyn_qjWnCSVEym3rHfR68FRgjYxAnMaflGXcWHWlbCxRN34O6iIXqX9pVy6T0BLTXmKv-AcZL5BkPt2IsSa8qvHbj8yQ1ye2_K1cVDfD4/s1600/LC812807_namngum.png" imageanchor="1" style="margin-left: auto; margin-right: auto;"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgANSBhR1Of7968JuGpvmrMFw1j3_UnRqg-IUcyn_qjWnCSVEym3rHfR68FRgjYxAnMaflGXcWHWlbCxRN34O6iIXqX9pVy6T0BLTXmKv-AcZL5BkPt2IsSa8qvHbj8yQ1ye2_K1cVDfD4/s320/LC812807_namngum.png" /></a></td></tr>
<tr><td class="tr-caption" style="text-align: center;">Zoom in - it is a super fast WMS!</td></tr>
</tbody></table><a name='more'></a><br />
Since I wanted to limit the dependencies as much as possible I had the following constraints:<br />
<ul><li>Input raster dataset must be a <a href="https://en.wikipedia.org/wiki/Portable_Network_Graphics">PNG</a> file, since cairo natively does not read nor write another raster file format</li>
<li>The input raster datasets needs to be georeferenced with an accompanying <a href="https://en.wikipedia.org/wiki/World_file">world file</a> to avoid an additional GIS library like <a href="http://gdal.org/">GDAL</a>. A world file contains exactly the same transformation parameters that are needed to place the raster image to the cairo context thus it can directly applied as transformation matrix.</li>
</ul>As with my GRASS GIS trials I wrapped everything in a simple <a href="http://docs.pylonsproject.org/projects/pyramid/">Pyramid</a> application.<br />
<br />
To avoid to read the input raster on every request it is read at server start up together with the world file and stored in a global variable as <a href="http://cairographics.org/documentation/pycairo/2/reference/patterns.html#class-surfacepattern-pattern">cairo.SurfacePattern</a>.<br />
<pre class="brush: python">def main(global_config, **settings):
"""
This function returns a Pyramid WSGI application.
"""
config = Configurator(settings=settings)
config.add_static_view('static', 'static', cache_max_age=3600)
config.add_route("index", "/")
config.add_route("wms", "/wms")
config.scan()
# Check pycairo capabilities
if not (cairo.HAS_IMAGE_SURFACE and cairo.HAS_PNG_FUNCTIONS):
raise SystemExit('cairo was not compiled with ImageSurface and PNG support')
# Read the input files
raster_data = cairo.ImageSurface.create_from_png(os.path.join(os.path.dirname(__file__), "data", "LC8128047_w_MASK.png"))
world_file = open(os.path.join(os.path.dirname(__file__), "data", "LC8128047_w_MASK.wld"), 'r')
world_file_content = world_file.read()
# Create a cairo pattern from the input raster layer
# and declare the pattern global
global raster_pattern
raster_pattern = cairo.SurfacePattern(raster_data)
# Parse the world file
lines = world_file_content.split("\n")
georef_scale_x = float(lines[0])
georef_scale_y = float(lines[3])
georef_trans_x = float(lines[4])
georef_trans_y = float(lines[5])
# Setup a new transformation matrix for the georeferenced raster layer
scaler = cairo.Matrix()
scaler.scale(georef_scale_x, georef_scale_y)
scaler.translate(georef_trans_x / georef_scale_x, georef_trans_y / georef_scale_y)
scaler.invert()
raster_pattern.set_matrix(scaler)
# Set resampling filter
raster_pattern.set_filter(cairo.FILTER_FAST)
# Start the WSGI application
return config.make_wsgi_app()
</pre>Having this SurfacePattern a web map service view is very easy. The request parameters need to be read and a new cairo context created and transformed. Then I just paint the SurfacePattern to the context and return the newly created in-memory image. It is not even necessary to write it to the hard disk which saves time as well.<br />
<pre class="brush: python">@view_config(route_name="wms")
def wms_view(request):
start_time = time()
# Read the WMS parameters
layers = request.params.get("LAYERS", "").split(",")
bbox = [float(b) for b in request.params.get("BBOX", "").split(",")]
width = int(request.params.get("WIDTH"))
height = int(request.params.get("HEIGHT"))
# Create a new cairo surface with requested size
canvas = cairo.ImageSurface(cairo.FORMAT_ARGB32, width, height)
ctx = cairo.Context(canvas)
# Calculate the scale in x and y direction
scale_x = width / (bbox[2] - bbox[0])
scale_y = height / (bbox[1] - bbox[3])
# Create a forward transformation matrix
ctx_matrix = cairo.Matrix(y0=(-1) * bbox[3] * scale_y,
x0=(-1) * bbox[0] * scale_x,
yy=scale_y,
xx=scale_x)
ctx.set_matrix(ctx_matrix)
# Set the raster pattern as source and paint the canvas
ctx.set_source(raster_pattern)
ctx.paint()
# Create a new response object
response = Response(content_type="image/png", status="200 OK")
try:
# Write the canvas as png to the response body
canvas.write_to_png(response.body_file)
return response
except:
pass
finally:
log.debug("'cairo' took: %ss" % (time() - start_time))
</pre>Of course it is finally necessary to wrap up everything in a simple <a href="http://www.openlayers.org/">OpenLayers</a> web map:<br />
<pre class="brush: python">@view_config(route_name='index', renderer='templates/mytemplate.pt')
def index_view(request):
body = """
<html xml:lang="en" xmlns:tal="http://xml.zope.org/namespaces/tal" xmlns="http://www.w3.org/1999/xhtml">
<head>
<title>PyCairo Web Map Service</title>
<script src="http://openlayers.org/api/OpenLayers.js" type="text/javascript"></script>
<script type="text/javascript">
var map, wms_layer;
function init(){
var map = new OpenLayers.Map("map-div",{
projection: "EPSG:32648",
maxExtent: new OpenLayers.Bounds(208800.0, 1967189.0, 430230.0, 2188289.0),
numZoomLevels: 8,
});
var wms_layer = new OpenLayers.Layer.WMS("GRASS WMS", '/wms', {
layers: "LC8128047_w_MASK"
});
map.addLayer(wms_layer);
map.zoomToMaxExtent();
}
</script>
</head>
<body onload="init()">
<div id="map-div" style="bottom: 0px; left: 0px; position: absolute; right: 0px; top: 0px;">
</div>
</body>
</html>"""
return Response(body=body, content_type="text/html", status=200)
</pre>Conclusions: Once more PyCairo didn't disappoint! To return a single raster tile it only takes around 0.01 to 0.04 seconds which is really fast.<br />
As with the previous GRASS GIS service reprojection on the fly is not supported.<br />
<br />
Finally another screenshot of the resulting map with a <a href="http://www.webrian.ch/2014/07/imitating-mapbox-cloudless-atlas.html">cloudless Landsat 8 scene</a>:<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/AVvXsEhtsaGBOuNya7_ootNBmGnydwbNeXRslXdvQjfs_-0zA2p4igOJ1tmKwe4ffKgkUjbvRdxUqDD6fR3J15dyAP87fE3v__Rg4yUuIkxhRMUSzjo6P5KENpP0kq3fuGT-ybluiUV_xGTtOUk/s1600/LC812807_ovrvw.png" imageanchor="1" style="margin-left: auto; margin-right: auto;"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhtsaGBOuNya7_ootNBmGnydwbNeXRslXdvQjfs_-0zA2p4igOJ1tmKwe4ffKgkUjbvRdxUqDD6fR3J15dyAP87fE3v__Rg4yUuIkxhRMUSzjo6P5KENpP0kq3fuGT-ybluiUV_xGTtOUk/s320/LC812807_ovrvw.png" /></a></td></tr>
<tr><td class="tr-caption" style="text-align: center;">The cloudless Landsat 8 scene</td></tr>
</tbody></table>webrianhttp://www.blogger.com/profile/15325993494848772844noreply@blogger.com0Bolikhamsai, Laos18.158901213171138 103.0840301513671918.098549213171136 103.00334915136719 18.21925321317114 103.16471115136719tag:blogger.com,1999:blog-8633623474665314466.post-65901360539765496632014-07-22T09:56:00.000+02:002014-07-22T10:02:37.795+02:00Imitating MapBox' Cloudless AtlasInspired by <a href="http://sgillies.github.io/scipy-2014-rasterio/#/">several</a> <a href="http://www.wired.com/2013/05/a-cloudless-atlas/">posts</a> which explain <a href="https://www.mapbox.com/">MapBox</a>' marvellous Cloudless Atlas, I wanted to imitate a cloudless atlas by myself. Limited in time and computer power I had to do it much simpler of course. But using first class spatial Python libraries I could achieve quite nice results with a simple Python script.<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/AVvXsEjfpAuCH7X4SeBThy4LQCdMtfcbyeUTCrVfiBv_fTVXVsdD0KnSTvLunxVcukpYmP-MHZa018ptjb7NO9v2ipeHh3i9ZvtrWo4IOSkZrLEdE90hf0gjhccacrDa-rVWOA2zTqm0gaoIOMA/s1600/LC8128047_sm.jpg" imageanchor="1" style="margin-left: auto; margin-right: auto;"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjfpAuCH7X4SeBThy4LQCdMtfcbyeUTCrVfiBv_fTVXVsdD0KnSTvLunxVcukpYmP-MHZa018ptjb7NO9v2ipeHh3i9ZvtrWo4IOSkZrLEdE90hf0gjhccacrDa-rVWOA2zTqm0gaoIOMA/s320/LC8128047_sm.jpg" height="320" width="320" /></a></td></tr>
<tr><td class="tr-caption" style="text-align: center;">The result: a cloudless Landsat 8 scene</td></tr>
</tbody></table><a name='more'></a><br />
As input images I use all available <a href="http://landsat.usgs.gov/landsat8.php">Landsat 8</a> images for one scene. I store the image names in an array and define the upper left and lower right corner:<br />
<pre class="brush: python;">filenames = [
"LC81280472013103LGN01",
"LC81280472013119LGN01",
"LC81280472013151LGN00",
"LC81280472013167LGN00",
"LC81280472013279LGN00",
"LC81280472013295LGN00",
"LC81280472013327LGN00",
"LC81280472013343LGN00",
"LC81280472014010LGN00",
"LC81280472014026LGN00",
"LC81280472014058LGN00",
"LC81280472014090LGN00",
"LC81280472014106LGN00",
"LC81280472014154LGN00"
]
# upper left corner for LC8128047:
ul = (208800, 2188289)
# lower right corner for LC8128047:
lr = (430215, 1967195)
</pre><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/AVvXsEhpYyhYrUS9Lqg2SI1uU7OsIy5D_ruhSw2uBMb-_-P6XP0YxuO_SOYDDq1qZMLpZXLw9HR1r-eis6cdsx1H_rBtet3H1_dVTfUkpdL__0sv4uzC3GXULyuViIrJDgobXxokB31GjjCth9o/s1600/lc8_extent_original.jpg" imageanchor="1" style="margin-left: auto; margin-right: auto;"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhpYyhYrUS9Lqg2SI1uU7OsIy5D_ruhSw2uBMb-_-P6XP0YxuO_SOYDDq1qZMLpZXLw9HR1r-eis6cdsx1H_rBtet3H1_dVTfUkpdL__0sv4uzC3GXULyuViIrJDgobXxokB31GjjCth9o/s320/lc8_extent_original.jpg" height="90" width="320" /></a></td></tr>
<tr><td class="tr-caption" style="text-align: center;">A small extent in chronological order</td></tr>
</tbody></table>Then all input images are read and stored aa <a href="http://www.numpy.org/">numpy</a> array to a list. This is done for each color band:<br />
<pre class="brush: python;">for file in filenames:
ds = gdal.Open("%s/%s/%s.jpg" % (datadir, file, file), gdalconst.GA_ReadOnly)
geotransform = ds.GetGeoTransform()
projection = ds.GetProjection()
originX = geotransform[0]
originY = geotransform[3]
pixelWidth = geotransform[1]
pixelHeight = geotransform[5]
ulx = int((ul[0] - originX) / pixelWidth)
uly = int((ul[1] - originY) / pixelHeight)
lrx = int((lr[0] - originX) / pixelWidth)
lry = int((lr[1] - originY) / pixelHeight)
redband = ds.GetRasterBand(1)
redpixels.append(redband.ReadAsArray(ulx, uly, (lrx - ulx), (lry - uly)))
greenband = ds.GetRasterBand(2)
greenpixels.append(greenband.ReadAsArray(ulx, uly, (lrx - ulx), (lry - uly)))
blueband = ds.GetRasterBand(3)
bluepixels.append(blueband.ReadAsArray(ulx, uly, (lrx - ulx), (lry - uly)))
</pre>I also prepare the output image with the same dimensions as the input images:<br />
<pre class="brush: python;">height = len(redpixels[0])
width = len(redpixels[0][0])
out_redpixels = np.empty((height, width), dtype=int)
out_greenpixels = np.empty((height, width), dtype=int)
out_bluepixels = np.empty((height, width), dtype=int)
</pre>Now comes the main part where the pixels from each band and image are sorted. Then the second darkest pixel is selected for the output image. Since the darkest pixels often are cloud shadows, the second or even third darkest pixel can be chosen, depending on the number of available images and cloud cover of the processed scene.<br />
<pre class="brush: python;">for col in range(width):
for row in range(height):
# Choose the second darkest pixel
j = 1
out_redpixels[row][col] = np.sort(np.array([rp[row][col] for rp in redpixels]))[j]
out_greenpixels[row][col] = np.sort(np.array([gp[row][col] for gp in greenpixels]))[j]
out_bluepixels[row][col] = np.sort(np.array([bp[row][col] for bp in bluepixels]))[j]
</pre>After that the output bands are written and georeferenced to the final image using <a href="http://www.gdal.org/">GDAL</a>.<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/AVvXsEgMmB7SuD8Fz1l5hLaJDhoB50ZpAwOFdVaXIhBNcCVUEPYMEsjGKabKf43oQmCPJJOivub7m09lr4OqOTcBr9U4lBhReuCYaLJK8QlNKPMQdFTwFn_YdeUozTNc1wZZ49M5-sLQUuFDtVs/s1600/lc8_extent_reshuffled.jpg" imageanchor="1" style="margin-left: auto; margin-right: auto;"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgMmB7SuD8Fz1l5hLaJDhoB50ZpAwOFdVaXIhBNcCVUEPYMEsjGKabKf43oQmCPJJOivub7m09lr4OqOTcBr9U4lBhReuCYaLJK8QlNKPMQdFTwFn_YdeUozTNc1wZZ49M5-sLQUuFDtVs/s320/lc8_extent_reshuffled.jpg" height="90" width="320" /></a></td></tr>
<tr><td class="tr-caption" style="text-align: center;">A small extent with pixels ordered by value</td></tr>
</tbody></table>In another test I also tried to exclude cloud pixels using the cloud mask and calculating the mean or the median from the remaining pixels, but the result was not satisfying with this limited number of images.webrianhttp://www.blogger.com/profile/15325993494848772844noreply@blogger.com0Vientiane, Laos18.986817585497494 102.8237915039062518.746430585497496 102.50106800390625 19.227204585497493 103.14651500390625tag:blogger.com,1999:blog-8633623474665314466.post-91912030321210693122014-02-22T16:43:00.001+01:002014-02-22T16:44:31.558+01:00Building ZOO with Java on Ubuntu 12.04I tried to build the <a href="http://zoo-project.org/">ZOO WPS platform</a> with Java support on Ubuntu 12.04 and it turned out to be not so easy. This is how I managed it with <a href="http://openjdk.java.net/">OpenJDK 7</a> (packages <i>openjdk-7-jdk</i> and <i>openjdk-7-jre</i>):<br />
<pre class="brush: bash">./configure LDFLAGS=-L/usr/lib/jvm/java-7-openjdk-amd64/jre/lib/amd64/server --with-java=/usr/lib/jvm/java-7-openjdk-amd64</pre>For unknown reasons the configuration writes wrong JAVALDFLAGS to the Makefile, thus the Makefile needs to be altered as follows. Replace the line where the JAVALDFLAGS are defined:<br />
<pre class="brush: plain">JAVALDFLAGS=-L/usr/lib/jvm/java-7-openjdk-amd64/jre/lib/amd64/client/ -ljvm -lpthread</pre>with<br />
<pre class="brush: plain">JAVALDFLAGS=-L/usr/lib/jvm/java-7-openjdk-amd64/jre/lib/amd64/server/ -ljvm -lpthread</pre>then follow the normal <a href="http://zoo-project.org/docs/kernel/install-debian.html#installation-workflow">instructions</a> with<br />
<pre class="brush: plain">make</pre>webrianhttp://www.blogger.com/profile/15325993494848772844noreply@blogger.com0tag:blogger.com,1999:blog-8633623474665314466.post-73562140304561978172013-12-19T10:51:00.001+01:002014-08-15T16:23:52.663+02:00GRASS GIS Web Map Service with PyramidLately I explored the capabilities of the <a href="http://grasswiki.osgeo.org/wiki/Python">Python bindings</a> in <a href="http://grass.osgeo.org/">GRASS GIS</a>. Since I'm quite familiar with Python web frameworks it was an easy task to wrap the <a href="http://grass.osgeo.org/grass70/manuals/d.mon.html">d.mon module</a> and create a <a href="https://en.wikipedia.org/wiki/Web_Map_Service">Web Map Service</a> based on a GRASS location.<br />
<br />
Here I'll show how I implemented the WMS using the Pyramid framework in a single file.<br />
<a name='more'></a><br />
Below is the core function of the WMS, it is basically a wrapper for the d.mon module.<br />
First of all the GRASS Python modules are imported and a new GRASS environment is started. After setting the current region from the bounding box parameters in the WMS request, a new <a href="http://grass.osgeo.org/grass70/manuals/cairodriver.html">Cairo</a> monitor is started with d.mon. The output is written to a temporary file which will be returned.<br />
<pre class="brush: python;">def _grass_wms(layers=[], bbox=[], width=256, height=256):
# Path to GRASS installation
gisbase = "/usr/local/grass-7.0.svn"
# Set the environment variable
os.environ["GISBASE"] = gisbase
# Append the GRASS Python directory to the path
sys.path.append("%s/etc/python" % gisbase)
# Import the GRASS scripts
from grass.script import core as grass
from grass.script import setup as gsetup
# Open the mapset
gsetup.init(gisbase,
"/home/adrian/Data/GISBASE7",
"Landsat_128047",
"adrian")
vector_layers = grass.list_strings("vect")
raster_layers = grass.list_strings("rast")
grass.run_command("g.region", w=bbox[0], s=bbox[1], e=bbox[2], n=bbox[3])
# Create a temp file
tempfile = grass.tempfile()
filename = "%s.png" % tempfile
grass.run_command("d.mon",
start="cairo",
width=width,
height=height,
output=filename)
for layer in layers:
if layer in raster_layers:
grass.run_command("d.rast", map=layer, quiet=1)
elif layer in vector_layers:
grass.run_command("d.vect", map=layer, quiet=1, fcolor="0:0:255", color=None)
grass.run_command("d.mon", stop="cairo")
return filename
</pre><br />
Above method is called by a Pyramid view that gets the layers, bounding box and image width and height parameters from the request and returns the image:<br />
<pre class="brush: python;">@view_config(route_name="wms")
def wms_view(request):
layers = request.params.get("LAYERS", "").split(",")
bbox = request.params.get("BBOX", "").split(",")
width = request.params.get("WIDTH")
height = request.params.get("HEIGHT")
try:
filename = _grass_wms(layers=layers,
bbox=bbox,
width=width,
height=height)
f = open(filename, "r+")
return Response(body=f.read(), content_type="image/png")
except:
pass
</pre><br />
Furthermore I added a super simple index page which uses OpenLayers to show the WMS layer:<br />
<pre class="brush: python;">@view_config(route_name="index")
def index_view(request):
body = """
<html xml:lang="en"
xmlns:tal="http://xml.zope.org/namespaces/tal"
xmlns="http://www.w3.org/1999/xhtml">
<head>
<title>GRASS GIS Web Map Service</title>
<script src="http://openlayers.org/api/OpenLayers.js" type="text/javascript"></script>
<script type="text/javascript">
var map, wms_layer;
function init(){
var map = new OpenLayers.Map("map-div",{
projection: "EPSG:32648",
maxExtent: new OpenLayers.Bounds(235080, 2132280, 257100, 2147400),
numZoomLevels: 6,
});
var wms_layer = new OpenLayers.Layer.WMS("GRASS WMS", '/wms', {
layers: "namngum5_patch,namngum5_bw,namngum5_smooth"
},{
singleTile: true,
});
map.addLayer(wms_layer);
map.zoomToMaxExtent();
}
</script>
</head>
<body onload="init()">
<div id="map-div" style="height: 600px; width: 800px;">
</div>
</body>
</html>"""
return Response(body=body, content_type="text/html", status=200)
</pre><br />
Last but not least it needs some Pyramid / WSGI code to get the server running:<br />
<pre class="brush: python;">def main(global_config, ** settings):
""" This function returns a Pyramid WSGI application.
"""
config = Configurator(settings=settings)
config.add_static_view('static', 'static', cache_max_age=3600)
config.add_route("index", '/')
config.add_route("wms", '/wms')
config.scan()
return config.make_wsgi_app()
</pre><br />
Since this is a very quick implementation, there are some unsolved issues:<br />
<ul><li>Reprojection on the fly is not supported, i.e. a WMS request must be in the same projection as the GRASS location</li>
<li>Pyramid requests are not thread-safe i.e. if there are several requests setting a different region in GRASS the monitor output gets mixed up</li>
<li>probably others as well</li>
</ul><div>Last but not least two screenshots to illustrate the Web Map Service:</div><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/AVvXsEjjxedvfE5yJW_J3jHTFuDLaaQaO3TNJ7Em1C9h1iXjltwuEZ-GlpNgIU6OkPXKNB54uI5FC5UHuC02JpdAbwZx3iHMQks4wGtW2je2qnG27gTgDI08dVUcpz6SFO4G7zwAfkjRhn5N0W4/s1600/grass-wms-openlayers-preview.png" imageanchor="1" style="margin-left: auto; margin-right: auto;"><img border="0" height="281" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjjxedvfE5yJW_J3jHTFuDLaaQaO3TNJ7Em1C9h1iXjltwuEZ-GlpNgIU6OkPXKNB54uI5FC5UHuC02JpdAbwZx3iHMQks4wGtW2je2qnG27gTgDI08dVUcpz6SFO4G7zwAfkjRhn5N0W4/s400/grass-wms-openlayers-preview.png" width="400" /></a></td></tr>
<tr><td class="tr-caption" style="text-align: center;">The index page shows the WMS layer</td></tr>
</tbody></table><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/AVvXsEjS1FMFkjf-qizrZepTXtvQ01lZQl5n3vH8ZQG-ECusrUFGCf_PcnsNdw0TeHd8Fy6y_ggnsmI3UYICdw2crI10QxhvVTe1lb0zOeA5FrcL02Of3Ac9NEmkAKI8PRl44Os_syc8aLnghTE/s1600/grass-preview.png" imageanchor="1" style="margin-left: auto; margin-right: auto;"><img border="0" height="292" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjS1FMFkjf-qizrZepTXtvQ01lZQl5n3vH8ZQG-ECusrUFGCf_PcnsNdw0TeHd8Fy6y_ggnsmI3UYICdw2crI10QxhvVTe1lb0zOeA5FrcL02Of3Ac9NEmkAKI8PRl44Os_syc8aLnghTE/s400/grass-preview.png" width="400" /></a></td></tr>
<tr><td class="tr-caption" style="text-align: center;">The same location and layers in GRASS GIS</td></tr>
</tbody></table>webrianhttp://www.blogger.com/profile/15325993494848772844noreply@blogger.com0Vientiane, Laos19.220396465588731 102.5903320312518.98049596558873 102.26760853125 19.460296965588732 102.91305553125tag:blogger.com,1999:blog-8633623474665314466.post-37239840342946249142012-11-28T16:09:00.000+01:002013-09-24T16:17:50.355+02:00New downloads availableI'm happy to announce additional Shapefile downloads on <a href="http://www.openstreetmap.la/">openstreetmap.la</a>:<br />
those include buildings, national parks, provincial and national boundaries for Laos and Cambodia.<br />
Please find all files at <a href="http://www.openstreetmap.la/downloads">openstreetmap.la/downloads</a>.<br />
<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/AVvXsEhWS8WKWUE53Pr8qhlC5szhCMYmwFmbv031h-PQpRp9bXQJn2SbHhi3x1Etzs3sZyyr5qVNPEeIlMSeXyaAC5sk3xSmB7xJ_y_hXCen_9GD5aIMe-0pOY2GQUbqazMUMltnPgaQKGPt47Y/s1600/buildings.png" imageanchor="1" style="margin-left: auto; margin-right: auto;"><img border="0" height="313" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhWS8WKWUE53Pr8qhlC5szhCMYmwFmbv031h-PQpRp9bXQJn2SbHhi3x1Etzs3sZyyr5qVNPEeIlMSeXyaAC5sk3xSmB7xJ_y_hXCen_9GD5aIMe-0pOY2GQUbqazMUMltnPgaQKGPt47Y/s400/buildings.png" width="400" /></a></td></tr>
<tr><td class="tr-caption" style="text-align: center;">Two villages close to Pakxan with buildings displayed in QGIS.<br />
Data © OpenStreetMap contributors</td></tr>
</tbody></table>
webrianhttp://www.blogger.com/profile/15325993494848772844noreply@blogger.com0tag:blogger.com,1999:blog-8633623474665314466.post-36139724362645216872012-11-19T09:10:00.000+01:002012-11-20T09:16:01.735+01:00Nam Ngum 2 dam three year agoYesterday I found by chance a picture that is exactly three years old. I took it when I flew from Phonsavan back to Vientiane and it shows the <a href="http://en.wikipedia.org/wiki/Mekong_River_Basin_Hydropower#Existing_hydropower_infrastructure">Nam Ngum 2 dam</a> under construction.<br />
<br />
<div class="separator" style="clear: both; text-align: center;">
</div>
<div class="separator" style="clear: both; text-align: center;">
</div>
<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/AVvXsEjWql2HdqxB2YQDd4LEHQ4VAXibXgDZX_xBiyXrX5kUhjkEC9CCc7XF4-WWMaEfk2aXR7cmrUyBfvUYMbohNor3jqQXHOuk7Pw4M0Ni5FYDQgVKsu6g_QXqPEntQmFjGCU88SAYYO6MHEc/s1600/P1010708_enhanced.jpg" imageanchor="1" style="margin-left: auto; margin-right: auto;"><img border="0" height="300" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjWql2HdqxB2YQDd4LEHQ4VAXibXgDZX_xBiyXrX5kUhjkEC9CCc7XF4-WWMaEfk2aXR7cmrUyBfvUYMbohNor3jqQXHOuk7Pw4M0Ni5FYDQgVKsu6g_QXqPEntQmFjGCU88SAYYO6MHEc/s400/P1010708_enhanced.jpg" width="400" /></a></td></tr>
<tr><td class="tr-caption" style="text-align: center;">Nam Ngum 2 dam under construction on 21. November 2009</td></tr>
</tbody></table>
<br />
<a name='more'></a><br />
This is how it looks on the latest <a href="http://www.bing.com/maps/?v=2&cp=18.751736~102.769915&lvl=15&dir=0&sty=h&form=LMLTCC">Bing maps</a> with overlaid <a href="http://openstreetmap.la/?lat=18.747302592269225&lon=102.78087615966797&zoom=13&mlat=18.753073158054768&mlon=102.77503967285155&lang=en">OpenStreetMap data</a>:<br />
<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/AVvXsEgkM6upIVHN0F8iaRAG-0Cqc-uATb7pCOMDNLMPGt_UGpal1AvK2A6KBGOUCJfHLVTlpoCV5IZLNS_wGhKrn7wppvY9OGLZvGEMerwg7zXIiGRB1_F95OLA-noHif0Av11fTZnftN2-nds/s1600/josm_namngum2.jpg" imageanchor="1" style="margin-left: auto; margin-right: auto;"><img border="0" height="300" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgkM6upIVHN0F8iaRAG-0Cqc-uATb7pCOMDNLMPGt_UGpal1AvK2A6KBGOUCJfHLVTlpoCV5IZLNS_wGhKrn7wppvY9OGLZvGEMerwg7zXIiGRB1_F95OLA-noHif0Av11fTZnftN2-nds/s400/josm_namngum2.jpg" width="400" /></a></td></tr>
<tr><td class="tr-caption" style="text-align: center;">Nam Ngum 2 dam on Bing satellite maps with OpenStreetMap data</td></tr>
</tbody></table>
<br />webrianhttp://www.blogger.com/profile/15325993494848772844noreply@blogger.com0Nam Ngum 2 dam, Laos18.746896206961409 102.7758979797363318.731860206961411 102.75615697973633 18.761932206961408 102.79563897973632tag:blogger.com,1999:blog-8633623474665314466.post-53045465365659786222012-11-07T22:48:00.000+01:002012-11-08T08:07:17.170+01:00Heat maps with PythonHeat maps are a <a href="http://gis.stackexchange.com/search?q=heatmap">quite often discussed</a> topic on <a href="http://gis.stackexchange.com/">gis.stackexchange.com</a>, probably because they are easy to understand and very appealing when showing spatial distributions of features. Working with a <a href="http://pyramid.readthedocs.org/">Python web framework</a>, I looked for a Python way to create heat maps. Finally I came up with a script that uses <a href="http://www.gdal.org/ogr/index.html">OGR</a> to read an input Shapefile and <a href="http://matplotlib.org/">matplotlib</a> to do the actual work.<br />
<a name='more'></a><br />
Snippets from the script:<br />
<br />
Read the input layer using OGR. Of course it can be in any format, that is supported by OGR.<br />
<pre class="brush: python;"># get the shape file driver
shpDriver = osgeo.ogr.GetDriverByName('ESRI Shapefile')
# Open shape file with points
shpFile = 'activities.shp'
dataSource = shpDriver.Open(shpFile, 0)
if dataSource is None:
print 'Could not open file ' + shpFile
sys.exit(1)
# Get the layer
layer = dataSource.GetLayer()</pre>
<br />
A bounding box and the number of cells need to be defined to calculate the resolution.<br />
<pre class="brush: python;"># The global bounding box
xmin = -180.0
ymin = -90.0
xmax = 180.0
ymax = 90.0
# Number of columns and rows
nbrColumns = 72.0
nbrRows = 36.0
# Caculate the cell size in x and y direction
csx = (xmax - xmin) / nbrColumns
csy = (ymax - ymin) / nbrRows</pre>
<br />
Knowing the cell size the script loops the whole extent and counts the number of feature in each cell. The number of features is stored in a two-dimensional array.<br />
<pre class="brush: python;">rows = []
i = ymax
while i > ymin:
j = xmin
cols = []
while j < xmax:
# Set a spatial filter
layer.SetSpatialFilterRect(j, (i-csy), (j+csx), i)
# And count the features
cols.append(layer.GetFeatureCount())
j += csx
rows.append(cols)
i -= csy</pre>
After preparing the data it's time for plotting this data with matplotlib. First a new figure is instantiated and the size is set. Then a new <a href="http://matplotlib.org/api/axes_api.html?highlight=imshow#matplotlib.axes.Axes">Axes</a> object is created and added to the figure.<br />
<pre class="brush: python;"># Create a new figure
import matplotlib.pyplot as plt
fig = plt.figure()
# Set the size
fig.set_size_inches((outputWidth / 100.0), ((outputWidth / 2.0) / 100.0))
ax = plt.Axes(fig, [0., 0., 1., 1.])
ax.set_axis_off()
fig.add_axes(ax)</pre>
<br />
Finally the array containing the number of features per cell is plotted using the <a href="http://matplotlib.org/api/axes_api.html?highlight=imshow#matplotlib.axes.Axes.imshow">imshow</a> method. The result is saved to an image.<br />
<pre class="brush: python;"># Create the heatmap
ax.imshow(rows, extent=[xmin, xmax, ymin, ymax], interpolation='nearest')
# Write the image to a file
plt.savefig('%s.png' % outputName)</pre>
<br />
Optionally a <a href="https://en.wikipedia.org/wiki/World_file">world file</a> can be written to georeference this result and use it in a GIS.<br />
<pre class="brush: python;"># Write the world file to georeference the image
worldFile = open('%s.wld' % outputName, 'w')
# Pixel size in x direction
psx = (xmax - xmin) / outputWidth
worldFile.write('%s\n' % psx)
worldFile.write('0\n')
worldFile.write('0\n')
# Pixel size in y direction
pixelSizeY = (ymin - ymax) / (outputWidth / 2)
worldFile.write('%s\n' % pixelSizeY)
# Origin i.e. upper left corner
worldFile.write('%s\n' % (xmin + psx/2))
worldFile.write('%s\n' % (ymax + pixelSizeY/2))</pre>
<br />
Using data from the <a href="http://landportal.info/landmatrix">Land Matrix database</a>, an online public database that collects and visualize data on international land deals, I got the following results:<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/AVvXsEj9Z3x6lPn2Oym9YDOczYHmsDmM8zhCXa4owbNBEhfGVSuMEgYu5byVW-jkrHNIUWeoVd9NYDEgv4_sfmXL2MqWA4fPllRSu019Epygmdae9RniVErcB1Bfum4gjDRJrfOUiKAvkM4x8B0/s1600/activities.png" imageanchor="1" style="margin-left: auto; margin-right: auto;"><img border="0" height="200" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEj9Z3x6lPn2Oym9YDOczYHmsDmM8zhCXa4owbNBEhfGVSuMEgYu5byVW-jkrHNIUWeoVd9NYDEgv4_sfmXL2MqWA4fPllRSu019Epygmdae9RniVErcB1Bfum4gjDRJrfOUiKAvkM4x8B0/s400/activities.png" width="400" /></a></td></tr>
<tr><td class="tr-caption" style="text-align: center;">A heat map showing world wide land deals</td></tr>
</tbody></table>
Instead of the nearest, a bilinear interpolation can be used to get a fuzzier heat map. With some more matplotlib functions, features can also be plotted on top.<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/AVvXsEgU15Lc2xidA4Lbu2ZtgexMNw7RPyGj-VOEyV5sBE0hyv8Gbpa1maDatOKlPtj0bJrMSop1m_FHhKpEEYaVlJnyyc4mqawy7erfEYAHklNgx4xMKZgiWKaRwyqWAqbAb2Y_Jk9Ek71TQH4/s1600/worldwide_deals.png" imageanchor="1" style="margin-left: auto; margin-right: auto;"><img border="0" height="200" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgU15Lc2xidA4Lbu2ZtgexMNw7RPyGj-VOEyV5sBE0hyv8Gbpa1maDatOKlPtj0bJrMSop1m_FHhKpEEYaVlJnyyc4mqawy7erfEYAHklNgx4xMKZgiWKaRwyqWAqbAb2Y_Jk9Ek71TQH4/s400/worldwide_deals.png" width="400" /></a></td></tr>
<tr><td class="tr-caption" style="text-align: center;">Plotted land deals on top of a bilinear interpolated heat map</td></tr>
</tbody></table>
Last but not least the georeferenced map can be used for further mapping or analysis in a GIS.<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/AVvXsEgHbARPos8Q7NPHoO5Gqi6LCv4yB5-5PynibyUO5KvCelEq8Dgg1qEW2W1vkSozL5cQrEdiVCPa371g7YmXWY-MNGCPnS1bLnZ4GsTs5Rl0EJ-T9vHPnlOeQlfXZ9SSbHJuq8HbVMFRkBg/s1600/deals_with_countries.png" imageanchor="1" style="margin-left: auto; margin-right: auto;"><img border="0" height="198" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgHbARPos8Q7NPHoO5Gqi6LCv4yB5-5PynibyUO5KvCelEq8Dgg1qEW2W1vkSozL5cQrEdiVCPa371g7YmXWY-MNGCPnS1bLnZ4GsTs5Rl0EJ-T9vHPnlOeQlfXZ9SSbHJuq8HbVMFRkBg/s400/deals_with_countries.png" width="400" /></a></td></tr>
<tr><td class="tr-caption" style="text-align: center;">A heat map showing land deals with country borders from <a href="http://www.naturalearthdata.com/">Natural Earth</a></td></tr>
</tbody></table>
To me matplotlib seems a very powerful and comprehensive plotting library, and of course it offers a lot more options than I could show here.<br />
<br />webrianhttp://www.blogger.com/profile/15325993494848772844noreply@blogger.com0tag:blogger.com,1999:blog-8633623474665314466.post-4626054114773198492012-05-23T11:30:00.000+02:002012-05-24T08:43:45.885+02:00Editing GPX files with Garmin extensionsDuring my outdoor activities I often (not to say always) carry a <a href="http://www.garmin.com/">Garmin</a> GPS device with me. Nowadays current Garmin devices record also heart rate, pedal cadence or even temperature and write these data into a <a href="http://www.topografix.com/gpx.asp">GPX</a> file as <a href="http://www.topografix.com/GPX/1/1/#type_extensionsType">extensions</a>.<br />
<br />
After my last cycling trip I edited (like I normally do) my tracks in <a href="http://sourceforge.net/apps/mediawiki/viking/index.php?title=Main_Page">viking</a>, when I realized that all the additional Garmin data are lost. I tried to find another workflow without losing the data in the GPX extensions. Finally <a href="http://www.gdal.org/ogr/">OGR</a> was once more my big saver. For the following workflow OGR is required at least in version 1.8.<br />
<a name='more'></a><br />
The idea is to convert the GPX track to a Shapefile, edit it in a GIS and convert it back to GPX.<br />
<pre class="brush: bash;">ogr2ogr -sql "SELECT track_fid,track_seg_id,track_seg_point_id,ele,CAST(time AS character(32)),gpxtpx_TrackPointExtension FROM track_points" track_points.shp 20120517_Cycling.gpx
</pre>It is important to name the new layer <i>track_points</i> to make sure that OGR converts the points back to a track again, see also the <a href="http://www.gdal.org/ogr/drv_gpx.html">GPX driver documentation</a>.<br />
<br />
Then the track can be easily edited in any GIS:<br />
<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/AVvXsEijJl_FTjLPLzacNH9rg3FldkK0NrYFvDGvyxZh8DSQwHA_YmMHC1nDZEoHAbmD7UJuq-lO4po3cI3qVZejjB8aIPcsTfV4OTinPHXqdLhE_md-WaBkcLLw89jW7f9QwUJA816BFBzlTWk/s1600/track_points.png" imageanchor="1" style="margin-left: auto; margin-right: auto;"><img border="0" height="217" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEijJl_FTjLPLzacNH9rg3FldkK0NrYFvDGvyxZh8DSQwHA_YmMHC1nDZEoHAbmD7UJuq-lO4po3cI3qVZejjB8aIPcsTfV4OTinPHXqdLhE_md-WaBkcLLw89jW7f9QwUJA816BFBzlTWk/s320/track_points.png" width="320" /></a></td></tr>
<tr><td class="tr-caption" style="text-align: center;">Edit the Shapefile in <a href="http://www.qgis.org/">QGIS</a></td><td class="tr-caption" style="text-align: center;"><br />
</td></tr>
</tbody></table><br />
After editing the layer is converted back to a GPX file:<br />
<pre class="brush: bash;">ogr2ogr -f GPX -dsco GPX_USE_EXTENSIONS=YES -dsco GPX_EXTENSIONS_NS="gpxtpx" -dsco GPX_EXTENSIONS_NS_URL="http://www.garmin.com/xmlschemas/TrackPointExtensionv1.xsd" -sql "SELECT track_fid,track_seg_ AS track_seg_id,track_se_1 AS track_seg_point_id,ele,time,gpxtpx_Tra AS 'gpxtpx:TrackPointExtension' FROM track_points" 20120517_Cyling_edited.gpx track_points.shp
</pre><br />
Last but not least the GPX file can be formatted with <a href="http://xmlstar.sourceforge.net/">xmlstarlet</a>:<br />
<pre class="brush: bash;">xmlstarlet fo 20120517_Cycling_edited.gpx</pre><br />
<br />webrianhttp://www.blogger.com/profile/15325993494848772844noreply@blogger.com2tag:blogger.com,1999:blog-8633623474665314466.post-24130458493587840692012-04-18T13:09:00.001+02:002012-04-19T08:15:24.493+02:00Native SLD support in QGISThe current Quantum GIS master (i.e. the latest development version) <a href="http://hub.qgis.org/projects/quantum-gis/repository/revisions/84be523ba7722b540c9af0a18930ae905c65cd98">supports now</a> loading and saving of <a href="https://en.wikipedia.org/wiki/Styled_Layer_Descriptor">Styled Layer Descriptor (SLD)</a> styles in the layer properties dialog.<br />
<br />
The QGIS API has been extended by the new methods saveSldStyle and loadSldStyle in <a href="http://doc.qgis.org/api/classQgsMapLayer.html">QgsMapLayer</a> and writeSld and loadSld in <a href="http://doc.qgis.org/api/classQgsVectorLayer.html">QgsVectorLayer</a>. Currently these methods are not yet available in the Python API, but it is simple to add them in the corresponding <a href="http://riverbankcomputing.co.uk/software/sip/intro">SIP</a> files.
<br />
<a name='more'></a>
<br />
To write the current layer style to a file, type in the Python console:
<br />
<pre class="brush: python;">l = qgis.utils.iface.activeLayer()
l.saveSldStyle("/path/to/style.sld")
</pre>
It's is also possible to get the current layer style as string:
<br />
<pre class="brush: python;"># Get the active layer
l = qgis.utils.iface.activeLayer()
# Create a new XML document
document = QDomDocument()
header = document.createProcessingInstruction( "xml", "version=\"1.0\" encoding=\"UTF-8\"" )
document.appendChild( header )
# Create the SLD root element
root = document.createElementNS( "http://www.opengis.net/sld", "StyledLayerDescriptor" )
root.setAttribute( "version", "1.1.0" )
root.setAttribute( "xsi:schemaLocation", "http://www.opengis.net/sld http://schemas.opengis.net/sld/1.1.0/StyledLayerDescriptor.xsd" )
root.setAttribute( "xmlns:ogc", "http://www.opengis.net/ogc" )
root.setAttribute( "xmlns:se", "http://www.opengis.net/se" )
root.setAttribute( "xmlns:xlink", "http://www.w3.org/1999/xlink" )
root.setAttribute( "xmlns:xsi", "http://www.w3.org/2001/XMLSchema-instance" )
document.appendChild( root )
# Create the NamedLayer element
namedLayerNode = document.createElement( "NamedLayer" )
root.appendChild( namedLayerNode )
errorMsg = QString("")
l.writeSld(namedLayerNode, document, errorMsg)
print document.toString( 4 )
</pre>
It's still necessary to load SLD styles manually. SLD styles, unlike QGIS styles (.qml), saved with the same name in the same directory like a Shapefile are not (yet?) immediately loaded.webrianhttp://www.blogger.com/profile/15325993494848772844noreply@blogger.com3tag:blogger.com,1999:blog-8633623474665314466.post-27721920797598723382012-04-04T18:37:00.000+02:002012-04-06T00:11:27.288+02:00Line extraction with GRASS GISThis post presents another way how to use <a href="http://en.wikipedia.org/wiki/Landsat_7">Landsat imagery</a> to map features for <a href="http://www.openstreetmap.org/">OpenStreetMap</a> in an easy way with free and open source software only. After remapping the <a href="http://www.webrian.ch/2012/02/armchair-mapping-using-grass-gis.html">Nam Ngum 1</a>, <a href="http://www.webrian.ch/2012/01/armchair-mapping-with-landsat-imagery.html">Nam Leuk</a> reservoir and other lakes using a similar approach, I wanted to improve the <a href="http://www.openstreetmap.la/?lat=18.2657&lon=102.5518&zoom=13&mlat=18.2657&mlon=102.5518&lang=lo">Nam Ngum river</a>. The second largest river in Laos and an important tributary to the Mekong was in OpenStreetMap only mapped as a single line instead of an area, although it's more than hundred metre wide downstream of the Nam Ngum 1 reservoir.
<br />
<a name='more'></a><br />
As approach I chose to apply <a href="http://en.wikipedia.org/wiki/Edge_detection">edge detection</a> algorithms to the satellite images, namely the <a href="http://en.wikipedia.org/wiki/Sobel_operator">Sobel filter</a>. To get used to the rather new Python scripting in <a href="http://grass.osgeo.org/">GRASS GIS</a>, I wrote a Python script to preprocess the images. The Python scripting ability is a great addition to GRASS GIS, since it facilitates considerably the entry point to scripting.<br />
<br />
First all necessary modules needed to be imported.
<br />
<pre class="brush: python;">#!/usr/bin/env python
from grass.script import core as gcore
from grass.script import raster as graster
from string import join
import sys
</pre>
In the main function I defined a list of regions to reduce the computational area and thus processing time.
<br />
<pre class="brush: python;">def main():
# Create a list of regions along the river to reduce the computation time
regions = []
regions.append({'w': 279850, 'e': 300648, 's': 2003246, 'n': 2017059})
regions.append({'w': 260146, 'e': 280944, 's': 2001343, 'n': 2015156})
regions.append({'w': 242189, 'e': 262987, 's': 2005578, 'n': 2019392})
regions.append({'w': 234104, 'e': 243982, 's': 2017410, 'n': 2036064})
regions.append({'w': 231550, 'e': 245964, 's': 2033931, 'n': 2052559})
</pre>
As in my <a href="http://www.webrian.ch/2012/02/armchair-mapping-using-grass-gis.html">previous approach</a> I used again the near infrared band from Landsat images to distinct water bodies from land.
<br />
<pre class="brush: python;">
# The input map
inputMap = "20011130_B.4"
# Make sure the region especially the resolution is correctly set
gcore.run_command("g.region", rast=inputMap)
</pre>
A loop through the five regions applied the Sobel filter to the subimages. Alternatively it would also possible to use the GRASS GIS <a href="http://grass.osgeo.org/grass64/manuals/html64_user/i.zc.html">i.zc module</a>, that also serves the purpose to detect edges. The results of the Sobel filter and the i.zc module are very similar.
<br />
<pre class="brush: python;">
# A list to store the b/w maps
mapList = []
# Loop all regions
for region in regions:
# Temporary name for the current map
currentMap = "tmp.%s" % len(mapList)
gcore.debug("Current map: %s" % currentMap)
# Set the (sub-)region
gcore.run_command("g.region",
w = region['w'],
e = region['e'],
s = region['s'],
n = region['n'])
# The Sobel edge detection filter
sobelCmd = "%s = sqrt(pow((-1)*20011130_B.4[-1,-1] + 20011130_B.4[-1,1] - 2*20011130_B.4[0,-1] + 2*20011130_B.4[0,1] - 20011130_B.4[1,-1] + 20011130_B.4[1,1],2) + pow((-1)*20011130_B.4[-1,-1] - 2*20011130_B.4[-1,0] - 20011130_B.4[-1,1] + 20011130_B.4[1,-1] + 2*20011130_B.4[1,0] + 20011130_B.4[1,1],2))" % currentMap
graster.mapcalc(sobelCmd, overwrite = True)
#gcore.run_command("i.zc", input = "20011130_B.4", output= currentMap )
</pre>
In the next step I separated the image to 1 and null values using the map calculator.
<br />
<pre class="brush: python;">
bwCmd = "{m}.bw = if({m} >= 120, 1, null())".format(m=currentMap)
graster.mapcalc(bwCmd)
# Append the current map to patchInput
mapList.append("%s.bw" % currentMap)
</pre>
Then I set the region again to the full extent, patched the images together and applied the thinning to the raster file to prepare the vectorization.
<br />
<pre class="brush: python;">
# Set the region to the full extent
gcore.run_command("g.region", rast=join(mapList, ','))
# Merge the raster maps
patchedMap = "namngum_river"
gcore.run_command("r.patch", input = join(mapList, ','), output = patchedMap, overwrite=True)
# Prepare the map to vectorize with r.thin
thinMap = "namngum_river.thin"
gcore.run_command("r.thin", input = patchedMap, output = thinMap, overwrite=True)
# Vectorize the map
vectorMap = "namngum_river_raw"
gcore.run_command("r.to.vect", input = thinMap, output = vectorMap, type = 'line', overwrite=True)
</pre>
The final step in the main function was to clean up all temporary raster maps.
<br />
<pre class="brush: python;"> # Clean up all temporary raster maps
for i in range(len(mapList)):
maps = "tmp.%s,tmp.%s.bw" % (i, i)
gcore.run_command("g.remove", rast = maps)
</pre>
The main entry point for the whole script.
<br />
<pre class="brush: python;">if __name__ == "__main__":
options, flags = gcore.parser()
sys.exit(main())
</pre>
As next step I had to manually clean the resulting vector file. This can be done in different ways, I edited it <a href="http://grass.osgeo.org/grass64/manuals/html64_user/wxGUI.Vector_Digitizer.html">interactively in GRASS GIS</a> by deleting the false line scraps.<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/AVvXsEhKBD-wcXhvvfWVBnrG356H3l_OSGdQjlw-MZhdDl1kqxQZUQq0PkhqNaeUOpTSTXGIyhr3-NlIuYPkcbKz4RqZAd80cKtdKaC7aud7hFgRFcoRvjSlTUU7us1H6H38DGIDwx9p5Bpsxiw/s1600/GRASS_InteractiveEditing.png" imageanchor="1" style="margin-left: auto; margin-right: auto;"><img border="0" height="194" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhKBD-wcXhvvfWVBnrG356H3l_OSGdQjlw-MZhdDl1kqxQZUQq0PkhqNaeUOpTSTXGIyhr3-NlIuYPkcbKz4RqZAd80cKtdKaC7aud7hFgRFcoRvjSlTUU7us1H6H38DGIDwx9p5Bpsxiw/s320/GRASS_InteractiveEditing.png" width="320" /></a></td></tr>
<tr><td class="tr-caption" style="text-align: center;">The interactive editing in GRASS GIS</td></tr>
</tbody></table>
Last but not least I smoothed the lines with <a href="http://grass.osgeo.org/grass64/manuals/html64_user/v.generalize.html">v.generalize</a>, the module that advances to one of my favourite in GRASS GIS.
<br />
<pre class="brush: shell;">v.generalize input=namngum_river output=namngum_river_smooth method=snakes alpha=0.5 beta=0.5 threshold=0
</pre>
<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/AVvXsEheU6Fxb1Xt4o8ReN0iuVMBJVccWrUHgf80Vy7fcvU7_bHMolJNWM1QCGm78I9fHebwwkLJTWE3FN6mto-vCXDBM6ivnIJS_i_R6In6wQTDTLgQn86yPEZHCgh_MlZ5tx4MtWpP-an82d4/s1600/NamNgum_smooth.png" imageanchor="1" style="margin-left: auto; margin-right: auto;"><img border="0" height="240" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEheU6Fxb1Xt4o8ReN0iuVMBJVccWrUHgf80Vy7fcvU7_bHMolJNWM1QCGm78I9fHebwwkLJTWE3FN6mto-vCXDBM6ivnIJS_i_R6In6wQTDTLgQn86yPEZHCgh_MlZ5tx4MtWpP-an82d4/s320/NamNgum_smooth.png" width="320" /></a></td></tr>
<tr><td class="tr-caption" style="text-align: center;">Nam Ngum river before and after the smoothing</td></tr>
</tbody></table>
<br />
Finally I exported the vector map, converted it to OpenStreetMap's XML format and tagged it correctly in <a href="http://josm.openstreetmap.de/">JOSM</a>.
<br />
<pre class="brush: shell;">v.out.ogr input=namngum_river_smooth type=line dsn=~/namngum_river.shp
python ogr2osm.py -e 32648 -o ~/namngum_river.osm -a ~/namngum_river.shp
</pre>
Overlaying the extracted river on Bing imagery in JOSM shows that the result is very accurate.
<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/AVvXsEhs1nxUu_KTgS5jRJ5ZfDuf7dXgNIn0qeYPNubs_w822uJuJIgSn8CbZIO2midjXCFjuPOezlfKOyeT7_G59whPRJLEGWoutbuH3y5v3ry-gKC6XcFlHZff8C6PjTf8VJoI0NE932aJShM/s1600/NamNgumRiver_Bing.png" imageanchor="1" style="margin-left: auto; margin-right: auto;"><img border="0" height="187" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhs1nxUu_KTgS5jRJ5ZfDuf7dXgNIn0qeYPNubs_w822uJuJIgSn8CbZIO2midjXCFjuPOezlfKOyeT7_G59whPRJLEGWoutbuH3y5v3ry-gKC6XcFlHZff8C6PjTf8VJoI0NE932aJShM/s320/NamNgumRiver_Bing.png" width="320" /></a></td></tr>
<tr><td class="tr-caption" style="text-align: center;">The resulting river verified with Bing imagery</td></tr>
</tbody></table>webrianhttp://www.blogger.com/profile/15325993494848772844noreply@blogger.com2Thangone, Laos18.147972372237067 102.6299858093261718.132883872237066 102.61024480932618 18.163060872237068 102.64972680932617tag:blogger.com,1999:blog-8633623474665314466.post-23657550668151829992012-02-19T22:04:00.000+01:002012-02-19T22:04:40.808+01:00Armchair mapping using GRASS GISAfter reading my <a href="http://www.webrian.ch/2012/01/armchair-mapping-with-landsat-imagery.html">previous post</a> a friend asked me why I didn't use just <a href="http://grass.osgeo.org/">GRASS GIS</a>. That's really a valid question. Since the <a href="http://www.gdal.org/">GDAL/OGR</a> utilities are some of my daily tools, I'm more familiar with them than with GRASS. Out of interest I developed a workflow how to extract and vectorize waterbodies from Landsat images using (almost solely) GRASS GIS. The resulting features I wanted to upload to <a href="http://www.openstreetmap.org/">OpenStreetMap</a>.
<br />
This time the <a href="http://en.wikipedia.org/wiki/Nam_Ngum_Dam">Nam Ngum 1 hydro power reservoir</a> in Laos was my test area. Like Nam Leuk last time it was only roughly digitized and some of its islands needed to be remapped.
<br />
<a name='more'></a><br />
Since GRASS uses own raster and vector data formats, it was necessary to import the images:
<br />
<pre class="brush: bash;">GRASS 6.4.1 (Landsat_128047):~ > r.in.gdal input=Laos/Landsat7/LE71280472011042PFS00/L71128047_04720110211_B40.TIF output=20110211_B40
GRASS 6.4.1 (Landsat_128047):~ > r.in.gdal input=Laos/Landsat7/LE71280472011330EDC00/L71128047_04720111126_B40.TIF output=20111126_B40
GRASS 6.4.1 (Landsat_128047):~ > r.in.gdal input=Laos/Landsat7/LE71280472011026PFS00/L71128047_04720110126_B40.TIF output=20110126_B40
</pre>
Set the computational region to the area of interest:
<br />
<pre class="brush: bash;">GRASS 6.4.1 (Landsat_128047):~ > g.region w=231000 n=2082000 e=280000 s=2038000
</pre>
I merged the three images to remove the stripes. It was necessary to indicate the -z option to set zero values transparent.
<br />
<pre class="brush: bash;">GRASS 6.4.1 (Landsat_128047):~ > r.patch -z input=20111126_B40,20110211_B40,20110126_B40 output=NamNgum_B40
</pre>
<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/AVvXsEgcif4cbyvxKVkZCGpaScbuY4wQzp2tdSHpYoPhCu2maKClFA5Z8SJ8LS8wZF6BY0BwXcI2_56ggYzmuLP-SB4XS-T9Cs9bCh8czkWOt82RJc68d1HjT-P8mX26uipz07AvdXDryHQ5ad0/s1600/namngum_grass.png" imageanchor="1" style="margin-left: auto; margin-right: auto;"><img border="0" height="317" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgcif4cbyvxKVkZCGpaScbuY4wQzp2tdSHpYoPhCu2maKClFA5Z8SJ8LS8wZF6BY0BwXcI2_56ggYzmuLP-SB4XS-T9Cs9bCh8czkWOt82RJc68d1HjT-P8mX26uipz07AvdXDryHQ5ad0/s320/namngum_grass.png" width="320" /></a></td></tr>
<tr><td class="tr-caption" style="text-align: center;">Merged Landsat images in GRASS</td></tr>
</tbody></table>
<br />
Then I separated the image values with the map calculator. Alternatively the <i>r.classify</i> module is also an option.
<br />
<pre class="brush: bash;">GRASS 6.4.1 (Landsat_128047):~ > r.mapcalc 'NamNgum.bw = if(NamNgum_B40 <= 30, 1, null())'
</pre>
I took value 30 as threshold, but this depends on the image and needs to be evaluated. I set the threshold rather low to make sure that no land was classified as water. To assign mixed pixel along the lakeshore to water I expanded the waterbodies for one pixel with <i>r.grow</i>.
<br />
<pre class="brush: bash">GRASS 6.4.1 (Landsat_128047):~ > r.grow input=NamNgum.bw output=NamNgum.bw.grow
</pre>
Like in my previous workflow there were noisy pixels I wanted to get rid of. I had to export the image to edit it in GIMP. Because I set the <i>nodata</i> option to zero, no data pixels are shown transparent in GIMP.
<br />
<pre class="brush: bash;">GRASS 6.4.1 (Landsat_128047):~ > r.out.gdal input=NamNgum.bw.grow format=PNG output=~/NamNgum.bw.png nodata=0
</pre>
In GIMP I cleaned the image with the <i>eraser tool</i>, then I imported it again:<br />
<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/AVvXsEgINFHcMh4gWU01Aymyf25rn17JPKmQxCDMr1UtjEf7LeQtA5CRbkVuwRBLEn9aL1zHXOvtc1f_ybvXoKG_VGu8GcgRv0tam9ZuzNktib6_6GCiZmWmxty17QrOGAlLtYbTpxjGQOLbfos/s1600/namngum_gimp.png" imageanchor="1" style="margin-left: auto; margin-right: auto;"><img border="0" height="314" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgINFHcMh4gWU01Aymyf25rn17JPKmQxCDMr1UtjEf7LeQtA5CRbkVuwRBLEn9aL1zHXOvtc1f_ybvXoKG_VGu8GcgRv0tam9ZuzNktib6_6GCiZmWmxty17QrOGAlLtYbTpxjGQOLbfos/s320/namngum_gimp.png" width="320" /></a></td></tr>
<tr><td class="tr-caption" style="text-align: center;">Clean image in GIMP</td></tr>
</tbody></table>
<br />
<pre class="brush: bash;">GRASS 6.4.1 (Landsat_128047):~ > r.in.gdal input=~/NamNgum.bw.png output=NamNgum.clean
</pre>
The next step was to vectorize the raster data:
<br />
<pre class="brush: bash;">GRASS 6.4.1 (Landsat_128047):~ > r.to.vect input=NamNgum.clean output=NamNgum_raw feature=area
</pre>
As a last processing step I wanted to smooth the newly created vector geometry. The <i>v.generalize</i> module that generalizes and smooths geometries was the only processing step I did in GRASS in my <a href="http://www.webrian.ch/2012/01/armchair-mapping-with-landsat-imagery.html">previous workflow</a>.
<br />
<pre class="brush: bash;">GRASS 6.4.1 (Landsat_128047):~ > v.generalize input=NamNgum_raw output=NamNgum_smooth method=snakes alpha=0.5 beta=0.5
</pre>
To convert the data to the osm format I needed a Shapefile. I exported the layer with:
<br />
<pre class="brush: bash;">GRASS 6.4.1 (Landsat_128047):~ > v.out.ogr -c input=NamNgum_smooth dsn=~/namngum.shp type=area
</pre>
The conversion I did like last time with the <a href="http://wiki.openstreetmap.org/wiki/Ogr2osm">ogr2osm</a> python script.
<br />
<pre class="brush: bash;">python ogr2osm.py -e 32648 -a -o ~/namngum.osm ~/namgnum.shp
</pre>
Again it was necessary to finally tag the geometries in JOSM and correct the multipolygon relation. A task that needed to be done manually.
<br />
<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/AVvXsEhHx7yRRivLvKwBV4z5pS1DAjgheJtdPOTmEzOMzjmRIzVkW90-o2PrBCikJIk7yx4cXiNfH0k7l5ob3B4Mf6pRhaU9YQX2ukS4UnosTOgQPuoQdj0rYb48Q8Kg1eam9FLVKdfFqVHGVrA/s1600/namngum_comparison.png" imageanchor="1" style="margin-left: auto; margin-right: auto;"><img border="0" height="188" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhHx7yRRivLvKwBV4z5pS1DAjgheJtdPOTmEzOMzjmRIzVkW90-o2PrBCikJIk7yx4cXiNfH0k7l5ob3B4Mf6pRhaU9YQX2ukS4UnosTOgQPuoQdj0rYb48Q8Kg1eam9FLVKdfFqVHGVrA/s320/namngum_comparison.png" width="320" /></a></td></tr>
<tr><td class="tr-caption" style="text-align: center;">Newly created Nam Ngum in the background</td></tr>
</tbody></table>
<br />
<b>Conclusion</b><br />
The final result is basically the same like in the <a href="http://www.webrian.ch/2012/01/armchair-mapping-with-landsat-imagery.html">previous presented workflow</a>, but this one is probably even easier:<br />
With GRASS it is not necessary to clip the images, it's enough to set the computational region. Another needless step is to delete polygons that represents land area because I separated the raster pixels to value 1 or null (i.e. no data) and GRASS ignores the null data pixels during the vectorization.<br />
A disadvantage in GRASS are the additional import and export steps.webrianhttp://www.blogger.com/profile/15325993494848772844noreply@blogger.com0Nam Ngum hydropower reservoir18.526491895773912 102.628784179687518.044836395773913 101.9970701796875 19.008147395773911 103.2604981796875tag:blogger.com,1999:blog-8633623474665314466.post-52813254170549562342012-01-11T20:46:00.000+01:002012-01-11T22:49:44.230+01:00Armchair mapping with Landsat imageryIn this post I want to show step by step how I extracted hydro power reservoirs and lakes in Laos from up-to-date <a href="http://en.wikipedia.org/wiki/Landsat_7">Landsat imagery</a> for further use in <a href="http://www.openstreetmap.org/">OpenStreetMap</a>. The extraction is done straightforward with common GIS tools but without sophisticated
and/or proprietary remote-sensing software.<br />
As an example I take the <a href="http://www.openstreetmap.la/?mlat=18.47602&mlon=102.93159&zoom=12">Nam Leuk hydro power reservoir</a> in Laos. This reservoir needs to be <a href="http://wiki.openstreetmap.org/wiki/Remapping">remapped</a>, since the origin user declined the upcoming <a href="http://wiki.openstreetmap.org/wiki/Open_Database_License">license change</a>. Currently the reservoir is rather roughly generalized.
<br />
<a name='more'></a><br />
<b>Obtain the images</b><br />
There are different ways how to get Landsat images, my preferred way is the <a href="http://glovis.usgs.gov/">USGS Global Visualization Viewer</a> (short: GLOVIS), a web application that allows comfortable browsing through the Landsat images. Use the GLOVIS preview to make sure that the desired region is not covered by clouds.<br />
The Landsat Product Level 1 contains six visible bands, a thermal infrared and a panchromatic band. For my purpose, to extract water bodies, band 4 seems to be appropriate and sufficient.
<br />
<br />
<b>Remove the stripes</b><br />
Due to technical problems with the Landsat sensor there are undesirable stripes in all images. It is possible to remove these stripes using overlapping images taken on two different dates.<br />
First of all I clipped the GeoTIFF images to the required extent with GDAL to accelerate the whole image processing.
<br />
<pre class="brush: bash;">gdal_translate -projwin 272500 2050000 287000 2038000 L71128047_04720111126_B40.TIF L71128047_04720111126_B40_NamLeuk.tif
</pre>
<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/AVvXsEjNjxU8P-TuZXY_5YrNaJc78kCKEAcSRXg9DHqEDEDPMVXzhGAHilsBXdceZCoZ3iR4GysUk9H7iBEaJaa-HFCt6EqApbB92pJ3pXIW1mTsc4ekFuNZAFOJ3dQ7oRM40zsQ_9krKU-9gmc/s1600/namleuk1.png" imageanchor="1" style="margin-left: auto; margin-right: auto;"><img border="0" height="240" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjNjxU8P-TuZXY_5YrNaJc78kCKEAcSRXg9DHqEDEDPMVXzhGAHilsBXdceZCoZ3iR4GysUk9H7iBEaJaa-HFCt6EqApbB92pJ3pXIW1mTsc4ekFuNZAFOJ3dQ7oRM40zsQ_9krKU-9gmc/s320/namleuk1.png" width="320" /></a></td></tr>
<tr><td class="tr-caption" style="text-align: center;">Nam Leuk reservoir with stripes</td></tr>
</tbody></table>
<br />
I merged both images using <i><a href="http://www.gdal.org/gdal_merge.html">gdal_merge.py</a></i>, here it is important to specify the no-data value in the images, in this case 0.
<br />
<pre class="brush: bash;">gdal_merge.py -n 0 -o B40_NamLeuk.tif L71128047_04720110211_B40_NamLeuk.tif L71128047_04720111126_B40_NamLeuk.tif
</pre>
<b>Convert to black-white</b><br />
Next step is to convert the gray value image to a black-white 2-bit image to separate water from land. A threshold has to be set, all pixel values below the threshold are set to 0 (black), all pixels above to 1 (white). I used <a href="http://www.imagemagick.org/">ImageMagick</a>, but of course there are different ways how to perform that.
The threshold value had to be evaluated by <a href="http://en.wikipedia.org/wiki/Trial_and_error">trial and error</a>, the <i>display</i> command is very helpful to preview the image.
<br />
<pre class="brush: bash;">convert B40_NamLeuk.tif -threshold 11% png:- | display png:-
</pre>
After I had found a reasonable threshold I wrote it to a new PNG file
<br />
<pre class="brush: bash;">convert B40_NamLeuk.tif -threshold 11% png:namleuk_bw_dirty.png
</pre>
<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/AVvXsEgymrXq3gOZ5ufVDfU6NvHxDJwpO_2Ijp6CheAUGmWOxCcrfBsInNyEUcmr5ePoeNjVdXLlvvqbL_vfjmf0cPyi8GjLgf5vK6wOYi3wKu8vUq70alHKt-ekOXm1X33xJ674YHnjxGd0PHE/s1600/out3.png" imageanchor="1" style="margin-left: auto; margin-right: auto;"><img border="0" height="265" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgymrXq3gOZ5ufVDfU6NvHxDJwpO_2Ijp6CheAUGmWOxCcrfBsInNyEUcmr5ePoeNjVdXLlvvqbL_vfjmf0cPyi8GjLgf5vK6wOYi3wKu8vUq70alHKt-ekOXm1X33xJ674YHnjxGd0PHE/s320/out3.png" width="320" /></a></td></tr>
<tr><td class="tr-caption" style="text-align: center;">Still noisy black-white image</td></tr>
</tbody></table>
<br />
Now there were still unwanted, noisy black pixels that needed to be set to white. It was necessary to edit the image manually in <a href="http://www.gimp.org/">GIMP</a> using the pencil tool. The next image shows the clean black and white image.
<br />
<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/AVvXsEjAnvOz1Bdv89879BMjhtKgzZh7JwRp5iWBU8VldIJIR_Jx4URhJIPCCNNKS9dR1NKTtHAnDsrfuUYwea2gbt1v1pY621H1EQ0HL8U8fCXf17tdUMmJ2MFkAHBpmGXuDJeDPpYw1aS1vZA/s1600/namleuk_bw_clean.png" imageanchor="1" style="margin-left: auto; margin-right: auto;"><img border="0" height="265" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjAnvOz1Bdv89879BMjhtKgzZh7JwRp5iWBU8VldIJIR_Jx4URhJIPCCNNKS9dR1NKTtHAnDsrfuUYwea2gbt1v1pY621H1EQ0HL8U8fCXf17tdUMmJ2MFkAHBpmGXuDJeDPpYw1aS1vZA/s320/namleuk_bw_clean.png" width="320" /></a></td></tr>
<tr><td class="tr-caption" style="text-align: center;">Clean black-white image</td></tr>
</tbody></table>
<br />
ImageMagick dropped the georeference, that's why it was necessary to reassign it with GDAL and converting the image back to GeoTIFF.
<br />
<pre class="brush: bash">gdal_translate -a_srs EPSG:32648 -a_ullr 272475 2050005 286965 2038005 -of PNG namleuk_bw_clean.png namleuk_bw_clean.tif
</pre>
<b>Vectorization and smoothing</b><br />
Having a clean black and white referenced image I used <i><a href="http://www.gdal.org/gdal_polygonize.html">gdal_polygonize</a></i> to vectorize the image and create a new Shapefile.
<br />
<pre class="brush: bash;">gdal_polygonize.py namleuk_bw_clean.tif -f "ESRI Shapefile" namleuk.shp
</pre>
A new polygon was created for every contiguous region with the same value, i.e. there were polygons with attribute value 0 representing water and value 1 representing land.
<br />
<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/AVvXsEg-icbT2pZ1IgMusjX-Ej6dmRymlk58o09KedJYfAH29aDWK63Biddv0qgEAMXO-IblNJG8QgI4cGcbKDYuNhR3G4URCTNmqtnK_Nsq_S1Md-nF6D46c1pPMzNvR00IcuhQ6kj9zhNt1xg/s1600/namleuk_vect.png" imageanchor="1" style="margin-left: auto; margin-right: auto;"><img border="0" height="259" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEg-icbT2pZ1IgMusjX-Ej6dmRymlk58o09KedJYfAH29aDWK63Biddv0qgEAMXO-IblNJG8QgI4cGcbKDYuNhR3G4URCTNmqtnK_Nsq_S1Md-nF6D46c1pPMzNvR00IcuhQ6kj9zhNt1xg/s320/namleuk_vect.png" width="320" /></a></td></tr>
<tr><td class="tr-caption" style="text-align: center;">Vectorized Nam Leuk shape</td></tr>
</tbody></table>
<br />
Since I was only interested in water polygons I dropped all land polyons
<br />
<pre class="brush: bash;">ogr2ogr -where "DN=0" namleuk_clean.shp namleuk.shp</pre>
The polygons still showed the original pixels from the Landsat images, that's why I wanted to smooth the geometries.
<br />
<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/AVvXsEi-vox0sPCY1N9zf3pv4iNcQqOfoKy2DHKEinCn0XZMCZsVtUwXtYAcH5bDL1gepwAn2OQgG1yPCdIH2FxWRHWqGaOPpR4HZLAwlOXQPEFsgLX7f_C2YGADITBukXpE54uRVHgjjpxeTdw/s1600/aftersmooth.png" imageanchor="1" style="margin-left: auto; margin-right: auto;"><img border="0" height="203" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEi-vox0sPCY1N9zf3pv4iNcQqOfoKy2DHKEinCn0XZMCZsVtUwXtYAcH5bDL1gepwAn2OQgG1yPCdIH2FxWRHWqGaOPpR4HZLAwlOXQPEFsgLX7f_C2YGADITBukXpE54uRVHgjjpxeTdw/s320/aftersmooth.png" width="320" /></a></td></tr>
<tr><td class="tr-caption" style="text-align: center;">Before and after smoothing</td></tr>
</tbody></table>
<br />
The following steps I did in GRASS GIS. I created a new GRASS location named <i>namleuk</i> using the georeferenced image and imported the Shapefile.
<br />
<pre class="brush: bash;">GRASS 6.4.1 (namleuk):~ > v.in.ogr dsn=~/Desktop/namleuk_clean.shp output=namleuk_clean</pre>
To perform the smoothing I used module v.generalize:
<br />
<pre class="brush: bash;">GRASS 6.4.1 (namleuk):~ > v.generalize input=namleuk_clean output=namleuk_smooth method=snakes alpha=0.5 beta=0.5</pre>
After completing the smoothing I exported the layer again. It was necessary to specify the type as area.
<br />
<pre class="brush: bash;">GRASS 6.4.1 (namleuk):~/Data/GISBASE > v.out.ogr -s input=namleuk_smooth type=area dsn=~/Desktop/namleuk_smooth.shp</pre>
<b>Uploading to OSM</b><br />
For the final processing steps it was necessary to convert the Shapefile to an OSM file, clean up the OSM file and upload it.<br />
To convert the layer I used <a href="http://wiki.openstreetmap.org/wiki/Ogr2osm">ogr2osm</a>, a very powerful and helpful Python script, that gives full control of the attribute mapping. The script also reproject the input layer if it is not yet in geographic coordinates.
<br />
<pre class="brush: bash;">python ogr2osm.py -e 32648 -o ~/Desktop/namleuk.osm ~/Desktop/namleuk_smooth.shp
</pre>
After opening in <a href="http://josm.openstreetmap.de/">JOSM</a> I had to clean up manually the tags and the multipolygon relation, replace the existing reservoir with the newly extracted one and finally, finally upload it to the OpenStreetMap database.
<br />
<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/AVvXsEjc3olJrIcM6_5D7gXX9iFpeZcv_dMt2uVdH4yyy_cGZpkRYXE7V6C4_oLqsYZaoRcz26mOEz4-ojURxIbWsxboqr6rusMn0A6XGK2C_gUpmGZZ-o4uHHPzmMO2xgaqFx7vr8cRS6zf88Y/s1600/josm_compare.png" imageanchor="1" style="margin-left: auto; margin-right: auto;"><img border="0" height="188" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjc3olJrIcM6_5D7gXX9iFpeZcv_dMt2uVdH4yyy_cGZpkRYXE7V6C4_oLqsYZaoRcz26mOEz4-ojURxIbWsxboqr6rusMn0A6XGK2C_gUpmGZZ-o4uHHPzmMO2xgaqFx7vr8cRS6zf88Y/s320/josm_compare.png" width="320" /></a></td></tr>
<tr><td class="tr-caption" style="text-align: center;">Comparison hand-digitized (blue) and newly extracted reservoir (gray)</td></tr>
</tbody></table>webrianhttp://www.blogger.com/profile/15325993494848772844noreply@blogger.com0Nam Leuk, Laos18.457465477473214 102.9460144042968818.336969977473213 102.78808590429688 18.577960977473214 103.10394290429687tag:blogger.com,1999:blog-8633623474665314466.post-14105927299606384062011-12-10T20:19:00.001+01:002011-12-14T11:39:51.782+01:00Two years of mapping in LaosAfter more than two years of OpenStreetMap<i>ping</i> in Laos I want to look back and give a short review. Two years ago the OpenStreetMap in Vientiane but also in whole Laos was quite empty, as I illustrated <a href="http://www.webrian.ch/2011/09/vientiane-downtown-and-lonely-planet.html">in a previous post</a>.<br />
I started mapping Vientiane downtown rather unsystematic till I met another mapper and we decided to initiate the probably first mapping party in Vientiane end of October 2009!<br />
<br />
<a name='more'></a><br />
<br />
<div class="separator" style="clear: both; text-align: center;">
</div>
<div class="separator" style="clear: both; text-align: center;">
</div>
<div class="separator" style="clear: both; text-align: center;">
</div>
<div class="separator" style="clear: both; text-align: center;">
</div>
<div class="separator" style="clear: both; text-align: center;">
<iframe allowfullscreen='allowfullscreen' webkitallowfullscreen='webkitallowfullscreen' mozallowfullscreen='mozallowfullscreen' width='320' height='266' src='https://www.youtube.com/embed/trg9Qwbvc28?feature=player_embedded' frameborder='0'></iframe></div>
<br />
Thanks to this party I got in touch with the <a href="http://www.laozaa.com/">Laozaa</a> IT community, which proved to be a big support of OSM. At the Ubuntu 10.04 release and Laozaa birthday party I could introduce OSM to a wider audience of young and IT-interested people. The interest in OSM was big and a lot of questions raised concerning the <a href="http://www.laozaa.com/viewthread.php?tid=10084">use of OSM on mobile devices</a>. We agreed to do another mapping party soon. Shortly after this event I registered the domain <a href="http://www.openstreetmap.la/">openstreetmap.la</a> with the intention to create a Lao OSM portal.<br />
<br />
The next <a href="http://www.laozaa.com/viewthread.php?tid=10365">mapping party in June 2010</a> was rather a technical introduction to OSM and the <a href="http://josm.openstreetmap.de/">JOSM editor</a> than an outdoor data collection.<br />
<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/AVvXsEgkLXvsRDZIEeydgm7aA3u7jNFcw31V_U0XWPpAGT4mYSq5dHnP4jpyXvlqwxSIiaC9V7qgZiXV_gH6ZX70z-EVMGHWOqy8s3sM0yEyECg2sjPNd_KrziilYu47FxLzVUWNky5Y6NATyUY/s1600/mappingparty.jpg" imageanchor="1" style="margin-left: auto; margin-right: auto;"><img border="0" height="213" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgkLXvsRDZIEeydgm7aA3u7jNFcw31V_U0XWPpAGT4mYSq5dHnP4jpyXvlqwxSIiaC9V7qgZiXV_gH6ZX70z-EVMGHWOqy8s3sM0yEyECg2sjPNd_KrziilYu47FxLzVUWNky5Y6NATyUY/s320/mappingparty.jpg" width="320" /></a></td></tr>
<tr><td class="tr-caption" style="text-align: center;">Introduction to the JOSM editor</td></tr>
</tbody></table>
<br />
During Summer 2010 I developed the first version of the openstreetmap.la website, that I presented later to a wider audience at the <a href="http://barcampvientiane.org/">Barcamp in Vientiane</a>. Since this first version I implemented a lot of improvements like minimization of data traffic (due to slow internet connections in Laos), <a href="http://www.webrian.ch/2011/08/search-bar-on-openstreetmapla.html">places and points of interest search</a>, <a href="http://www.webrian.ch/2011/09/find-your-way-with-openstreetmapla.html">routing</a> and <a href="http://validator.w3.org/check?uri=http%3a%2f%2fwww.openstreetmap.la%2f">W3C compliance</a>.<br />
<br />
For reasons of limited server capacity I hardly accomplished a regularly updated map rendering, that's why I got in touch (first online, later also offline) with <a href="http://www.openstreetmap.org/user/stephankn">user Stephankn</a>, maintainer and developer of <a href="http://osm-tools.org/">osm-tools.org</a>, who rendered at this time a <a href="http://thaimap.osm-tools.org/">bilingual map of Thailand</a> in Thai and English. He kindly agreed to expand the bilingual map rendering to Laos (and even to Vietnam and Cambodia), therefore openstreetmap.la currently uses these map tiles. The tiles are now updated within a few minutes after data edits, thanks to Stephankn for this big support!<br />
<br />
In 2011 I increasingly started armchair mapping in Laos when I realized that a lot of bigger important features, mainly waterbodies, are still missing on the map. For the first time I downloaded <a href="http://glovis.usgs.gov/">Landsat imagery</a> to map newly built hydro power reservoirs like <a href="http://www.openstreetmap.la/?lat=18.88354892151393&lon=102.77572631835938&zoom=12&mlat=18.908560083449792&mlon=102.7657699584961">Nam Ngum 2</a> or <a href="http://www.openstreetmap.la/?lat=17.79184300887134&lon=105.11512756347656&zoom=11&mlat=17.800996047666967&mlon=105.1226806640625">Nam Theun 2</a> with its countless islands. OSM is as far as I know still the only dataset that maps the currently largest reservoir in Laos.<br />
<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/AVvXsEg6IPlM1z2buVsDOL2VjTj_7W0xIzWjWW_Eh-G6Qlpl402pG8FHdWS6aKH9_TcwwhFhGKEnDCpA8uEv7Qfc2bNq5IIFDExWge45TdFe77QNxfxZsuza_2meA5nWppGYhR_BCMSi-IK3vZ8/s1600/namtheun2.png" imageanchor="1" style="margin-left: auto; margin-right: auto;"><img border="0" height="190" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEg6IPlM1z2buVsDOL2VjTj_7W0xIzWjWW_Eh-G6Qlpl402pG8FHdWS6aKH9_TcwwhFhGKEnDCpA8uEv7Qfc2bNq5IIFDExWge45TdFe77QNxfxZsuza_2meA5nWppGYhR_BCMSi-IK3vZ8/s320/namtheun2.png" width="320" /></a></td></tr>
<tr><td class="tr-caption" style="text-align: center;">Nam Theun 2</td></tr>
</tbody></table>
<br />
Last but not least the OpenStreetMap wiki states on <a href="http://wiki.openstreetmap.org/wiki/Main_Page">its main page</a>, that the project was started to allow people using maps (and map data) <i>in creative, productive, or unexpected ways.</i> The <a href="http://www.fluid-forms.com/">wall clock</a> that pictures Vientiane downtown uses OpenStreetMap in one of those unexpected ways.<br />
<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/AVvXsEgQLq153ucWutTOty0rhoSnfeOe2yzVGsQ9TXXSm1VxCGsr0acAsIGQ97yC2OButf9mYBPh8ceAQGXn2XMTRM3pNZReJ4JkVjK5D4cnhIdxHMfO_SAFQUm9Q-DhvXlUK0_Quge2GN9gU8M/s1600/osm_vte_clock2.jpg" imageanchor="1" style="margin-left: auto; margin-right: auto;"><img border="0" height="320" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgQLq153ucWutTOty0rhoSnfeOe2yzVGsQ9TXXSm1VxCGsr0acAsIGQ97yC2OButf9mYBPh8ceAQGXn2XMTRM3pNZReJ4JkVjK5D4cnhIdxHMfO_SAFQUm9Q-DhvXlUK0_Quge2GN9gU8M/s320/osm_vte_clock2.jpg" width="320" /></a></td></tr>
<tr><td class="tr-caption" style="text-align: center;">Vientiane wall clock</td></tr>
</tbody></table>
<div class="separator" style="clear: both; text-align: center;">
</div>
<br />
<br />
<br />
<br />
<br />
<br />webrianhttp://www.blogger.com/profile/15325993494848772844noreply@blogger.com2Vientiane, Laos17.962769 102.61442917.9023505 102.535465 18.023187500000002 102.693393tag:blogger.com,1999:blog-8633623474665314466.post-33372059207895612392011-10-13T05:42:00.001+02:002011-10-14T04:21:37.785+02:00Save As SLD 0.3.0 releasedWhen I developed the Save As <a href="http://www.opengeospatial.org/standards/sld">SLD</a> <a href="http://www.qgis.org/">QGIS</a> plugin some months ago, I did it for my own purpose to publish simple <a href="http://en.wikipedia.org/wiki/Choropleth_map">choropleth maps</a> with <a href="http://www.geoserver.org/">GeoServer</a> and I implemented solely renderers and styling options I needed (and a little bit more). That's the reason why this plugin still lacks a lot of features.<br />
<a name='more'></a><br />
Later on I thought about a tighter integration with GeoServer than to save first a SLD file and then upload or copy-and-paste to GeoServer. Using <a href="http://docs.geoserver.org/stable/en/user/restconfig/overview.html">GeoServer's REST API</a> the SLD preview dialog features now the possibility to upload the style directly to GeoServer.<br />
<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/AVvXsEh9Flc13HA22L98EY9q6qmZkAQ5SGRUnISMhjXkIUZrWncIg1_Q01xUKZm_tJ6b6pns9FR8M20i7a7VCMkmBFOB2C8gYlnlQ1IAI5JDUU4aGHfnLRP57IFror916NS0SUwlgcHnPGR103U/s1600/Screenshot-SLD+Preview.png" imageanchor="1" style="margin-left: auto; margin-right: auto;"><img border="0" height="305" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEh9Flc13HA22L98EY9q6qmZkAQ5SGRUnISMhjXkIUZrWncIg1_Q01xUKZm_tJ6b6pns9FR8M20i7a7VCMkmBFOB2C8gYlnlQ1IAI5JDUU4aGHfnLRP57IFror916NS0SUwlgcHnPGR103U/s320/Screenshot-SLD+Preview.png" width="320" /></a></td></tr>
<tr><td class="tr-caption" style="text-align: center;">SLD Preview dialog</td></tr>
</tbody></table><br />
As parameters the GeoServer URL, user name and password and a layer name are required. The URL has to point to the REST API usually ending with <i>/geoserver/rest/</i>. The REST API was introduced with GeoServer version 2.1, for older version it can be installed as an extension.<br />
From the technical side: The network communication is done using <a href="http://doc.qt.nokia.com/stable/qnetworkaccessmanager.html">QNetworkAccessManager</a> and <a href="http://doc.qt.nokia.com/stable/qnetworkrequest.html">QNetworkRequest</a>s. A new style is created by <a href="http://en.wikipedia.org/wiki/POST_%28HTTP%29">POST</a>ing the SLD to <i>/geoserver/rest/styles.sld?name=</i>. Authentication is done using <a href="http://en.wikipedia.org/wiki/Basic_access_authentication">basic access authentication</a> as required by GeoServer. The server returns a 403 Forbidden status code if there is already a style with the same name.<br />
<br />
Important: Consider it as a proof-of-concept, it is far from complete! And as always, feedback is highly appreciated.<br />
<br />webrianhttp://www.blogger.com/profile/15325993494848772844noreply@blogger.com16tag:blogger.com,1999:blog-8633623474665314466.post-78006914873253630912011-09-30T13:58:00.000+02:002011-10-03T12:37:34.699+02:00Vientiane downtown and Lonely Planet updateWhen I arrived in <a href="http://en.wikipedia.org/wiki/Vientiane">Vientiane</a> in August 2009 the <a href="http://www.openstreetmap.org/">OpenStreetMap</a> was quite empty, there were only some unconnected roads traced from <a href="http://wiki.openstreetmap.org/wiki/Yahoo">Yahoo imagery</a> and a few points of interest. It was my second or third weekend when I took my GPS, a pencil and a sketchbook and started mapping the downtown. <br />
<a name='more'></a>Things change fast in Laos and restaurants come and go. From time to time I added a new hotel here and a new shop there but never systematically. That's why I decided some days ago to map again the downtown and update all the points of interest.<br />
See the following screenshots that shows the map development in Vientiane downtown:<br />
<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 rel="lightbox" href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjD_NIZVpqCYk1KVYcYMz7r8XkEQQzno6LCudrV4kmRQrp2nWTaWR0grZxNcAjXUXWa5hVJKbU28njNs0ebc83M0MeBr5NYXJ2AiMXUHzH6XCk3GHD92onOapOCNmUUJJdevBEVJUU6Cgw/s1600/vte_downtown_20090829.png" imageanchor="1" style="margin-left: auto; margin-right: auto;"><img alt="Vientiane downtown in August 2009" border="0" height="184" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjD_NIZVpqCYk1KVYcYMz7r8XkEQQzno6LCudrV4kmRQrp2nWTaWR0grZxNcAjXUXWa5hVJKbU28njNs0ebc83M0MeBr5NYXJ2AiMXUHzH6XCk3GHD92onOapOCNmUUJJdevBEVJUU6Cgw/s320/vte_downtown_20090829.png" width="320" /></a></td></tr>
<tr><td class="tr-caption" style="text-align: center;">Vientiane downtown in August 2009</td></tr>
</tbody></table><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 rel="lightbox" href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEi9Q_Sp658LcOte5aGe4vtxG3-gfVHw6koCUQnpi3CSdyBEN9eaGQODZIREiQVrdGKK8_rrD0sf4UymE3ZFUoX6cKu4yk4SdmbZvzAx23gDmqAv3hbMS-EhqRspvI_CC-OfvTrhvSf3qT4/s1600/vte_downtown_20110929.png" imageanchor="1" style="margin-left: auto; margin-right: auto;"><img alt="Vientiane downtown in September 2011" border="0" height="182" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEi9Q_Sp658LcOte5aGe4vtxG3-gfVHw6koCUQnpi3CSdyBEN9eaGQODZIREiQVrdGKK8_rrD0sf4UymE3ZFUoX6cKu4yk4SdmbZvzAx23gDmqAv3hbMS-EhqRspvI_CC-OfvTrhvSf3qT4/s320/vte_downtown_20110929.png" width="320" /></a></td></tr>
<tr><td class="tr-caption" style="text-align: center;">Vientiane downtown in September 2011</td></tr>
</tbody></table><br />
<br />
On an other day I took the latest<a href="http://www.lonelyplanet.com/"> Lonely Planet</a> edition to check if all the recommended eating, drinking and sleeping places are on OpenStreetMap. Unsurprisingly, there were still a few missing (and some are already gone again). So I put the GPS on my bicycle and rode across the city to search all these missing places. As always while mapping, I've discovered also this time new places and parts of the town, and that's definitively the fun part of mapping.<br />
Except one restaurant I found all missing points of interest. Now I'm happy to say that all Lonely Planet recommendations are on OpenStreetMap! Isn't that one reason more to use OpenStreetMap?webrianhttp://www.blogger.com/profile/15325993494848772844noreply@blogger.com0Rue Samsenthai, Vientiane, Laos17.962769 102.61442917.9023505 102.535465 18.023187500000002 102.693393tag:blogger.com,1999:blog-8633623474665314466.post-12984570934168625332011-09-21T18:15:00.001+02:002011-10-13T05:45:17.251+02:00Find your way with openstreetmap.laI'm really happy to announce the latest update on <a href="http://www.openstreetmap.la/">openstreetmap.la</a>: Routing has been integrated to the main page.<br />
<br />
As illustrated in a <a href="http://www.webrian.ch/2011/08/open-source-routing-machine-yet-another.html">previous post</a> I evaluated the capabilities of the Open Source Routing Machine <a href="http://www.webrian.ch/2011/07/set-up-pgrouting-with-openstreetmap.html">amongst</a> <a href="http://www.webrian.ch/2011/06/routing-with-openstreetmap-data-in.html">other</a> routing engines. The speed and the web-oriented architecture of the Open Source Routing Machine convinced me to use it as a backend to implement routing on openstreetmap.la.<br />
<a name='more'></a><br />
In the left panel there are two text fields to enter start and destination. Similar to the <a href="http://www.webrian.ch/2011/08/search-bar-on-openstreetmapla.html">search bar</a> (in the toolbar) these are auto-complete fields: When you start typing, any points of interest like villages, restaurants, shop etc. are suggested.<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/AVvXsEgy5oT2b1IXhFFNdiS6YhrAMPMctr2ZZOPE1v49fop6BA5lJSsRsnoDVR7H2WvwaQQZje3kd3ggtgEQYoqNjRdE87isuV6P_qMDAHsALGy5OaM7dqrNPOhSaMIUZE9nzcz-DXOT6PsWCvc/s1600/getdirections.png" imageanchor="1" rel="lightbox" style="margin-left: auto; margin-right: auto;"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgy5oT2b1IXhFFNdiS6YhrAMPMctr2ZZOPE1v49fop6BA5lJSsRsnoDVR7H2WvwaQQZje3kd3ggtgEQYoqNjRdE87isuV6P_qMDAHsALGy5OaM7dqrNPOhSaMIUZE9nzcz-DXOT6PsWCvc/s1600/getdirections.png" /></a></td></tr>
<tr><td class="tr-caption" style="text-align: center;">Enter start and destination</td></tr>
</tbody></table>The shortest route is calculated immediately and displayed on the map.<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/AVvXsEhuePsif84yrbIqRhsfZa9AqGkMFg8_IC5HDZ_lB_GoWCGW5GENynWvpq2EF9pZvgWb2sI4PzM1Bpf0u4aiIVx_I1dKcQ8A2h2uhV0RaWWXn_iBNkCEREDp4QTlYdFJ1kjfQeFnM97qX-M/s1600/calcRoute.png" imageanchor="1" rel="lightbox" style="margin-left: auto; margin-right: auto;"><img border="0" height="240" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhuePsif84yrbIqRhsfZa9AqGkMFg8_IC5HDZ_lB_GoWCGW5GENynWvpq2EF9pZvgWb2sI4PzM1Bpf0u4aiIVx_I1dKcQ8A2h2uhV0RaWWXn_iBNkCEREDp4QTlYdFJ1kjfQeFnM97qX-M/s320/calcRoute.png" width="320" /></a></td></tr>
<tr><td class="tr-caption" style="text-align: center;">Shortest route</td></tr>
</tbody></table><br />
To fine-adjust or change the start or destination, drag and drop one or both markers and the route is instantly updated.<br />
In the left panel detailed routing instructions are shown with heading directions, road names and trip distances. <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/AVvXsEh7u-nli1Ro_QyQUXEo2JaIfsDBpKhGkStk1GnSFF6tRF6CxgR_Y5dsb2J0Bxq5dWbE69EFBReNskhRKsef1YQ1NSCx3Jfs71svA2dWoWykb9X8hbNUTreN_cmKss_B4y2PEYrb4vCIhXQ/s1600/routeInstructions.png" imageanchor="1" rel="lightbox" style="margin-left: auto; margin-right: auto;"><img alt="Distance: 23.04km - Duration: 0 h 23 mins" border="0" height="260" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEh7u-nli1Ro_QyQUXEo2JaIfsDBpKhGkStk1GnSFF6tRF6CxgR_Y5dsb2J0Bxq5dWbE69EFBReNskhRKsef1YQ1NSCx3Jfs71svA2dWoWykb9X8hbNUTreN_cmKss_B4y2PEYrb4vCIhXQ/s320/routeInstructions.png" width="250" /></a></td></tr>
<tr><td class="tr-caption" style="text-align: center;">Detailed route instructions</td></tr>
</tbody></table>webrianhttp://www.blogger.com/profile/15325993494848772844noreply@blogger.com0tag:blogger.com,1999:blog-8633623474665314466.post-7186006154721586512011-09-12T16:34:00.000+02:002013-04-10T18:03:08.420+02:00Rivers and lakes downloads available on openstreetmap.laDuring my two-week trip through Laos I collected some tracks and <a href="http://en.wikipedia.org/wiki/Point_of_interest">points of interest</a>, of course ;). Last week I started to add them to OpenStreetMap when I realized that there are still a lot of rivers missing. I traced (rather randomly) some parts of rivers from <a href="http://en.wikipedia.org/wiki/Landsat_7">Landsat</a> or where available <a href="http://www.bing.com/community/site_blogs/b/maps/archive/2010/11/23/bing-engages-open-maps-community.aspx">Bing</a> imagery. <br />
<a name='more'></a><br />
To get a better overview of already traced rivers and lakes, I wanted to have these features as Shapefiles. So I added two more Shapefiles to the available downloads on <a href="http://www.openstreetmap.la/downloads">openstreetmap.la</a>, one contains rivers (line geometries) and the other contains lakes and larger rivers (polygon geometries). <br />
<br />
The downloads are available as usual on <a href="http://www.openstreetmap.la/">openstreetmap.la</a> or directly <a href="http://www.openstreetmap.la/downloads/laos/waterway_lines.shp.zip">from here (rivers)</a> or <a href="http://www.openstreetmap.la/downloads/laos/waterway_polygons.shp.zip">from here (lakes and larger rivers)</a>. <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/AVvXsEgizNK1eSQ4mGHAFGnVLLCQlLFLCSqBerYV2xJK-6n0qESfqTDOoyAkCd-1sVkfMiHmrj_W1AxdWkC3KPOg1q1wsKRLbYfOTP_j1l9adyKhVGdMw3Z2qtpMrCaIa_ILPtrmhSJ6NoV1gcg/s1600/waterways_osm.png" imageanchor="1" style="margin-left: auto; margin-right: auto;"><img border="0" height="259" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgizNK1eSQ4mGHAFGnVLLCQlLFLCSqBerYV2xJK-6n0qESfqTDOoyAkCd-1sVkfMiHmrj_W1AxdWkC3KPOg1q1wsKRLbYfOTP_j1l9adyKhVGdMw3Z2qtpMrCaIa_ILPtrmhSJ6NoV1gcg/s320/waterways_osm.png" width="320" /></a></td></tr>
<tr><td class="tr-caption" style="text-align: center;">Rivers and lakes from OpenStreetMap data as Shapefiles in QGIS </td></tr>
</tbody></table>
webrianhttp://www.blogger.com/profile/15325993494848772844noreply@blogger.com0