[OSGeoJapan-discuss] GRASSについて

Junta TAGUSARI j.tagusari @ eng.hokudai.ac.jp
2022年 1月 19日 (水) 20:27:14 PST


面白そうな話題なのでコメントさせてください。

処理のトレーサビリティ,たいへん重要な点だと思います。
「Rおじさん」である私は,主なGIS処理はRで行い(gdal/ogrを参照しているようです),Rmarkdownファイルからhtml形式のレポートを出力しています。拙いレポートですがサンプルを添付します。
このようにしておくと,処理と結果が一目瞭然でよいのではないかとおもいます。
なお,GIS以外のRで行う処理もすべて記録しておけます。

pythonにも似たような機能があったように思いますが,よく知りません。。。

QGISは,画像の出力時だとか,地図を見ながら地物を設定したいときに重宝しています。

横から失礼しました。

〒060-8628 札幌市北区北13条西8丁目
北海道大学大学院工学研究院 環境工学部門
地域環境研究室
田鎖順太 Junta Tagusari
Tel&Fax: 011-706-6872
Email: j.tagusari @ eng.hokudai.ac.jp

On 2022/01/20 13:18, y-iwasa @ limecon.co.jp wrote:
> しまださま
>
> 脱線しているかもしれませんがある意味自分も最近直面したトラブルだったのでトレーサビリティについては
> 痛感しています。正直シェルスクリプトのちゃんとした意味は知りませんがなんとなくコマンドプロンプトと同じような感じかと。
> 自分の場合、昨年からQGISを使いだしたレベルでGRASSについては今週から使いだした現状です。
>
> 特にQGISでラスタファイルの何らかの処理を行う時にパラメータ部分をデフォルトから変更していることを忘れており
> 同じ処理を行ったのに結果があまりにも違い一瞬あせりました。ただし変更したパラメータはネットで操作方法が
> 紹介されていた説明通りにやっていたので特に覚えていませんでした...
>
> 逆にコマンドプロンプトだとパラメータ等もすべて1行で実行できるのでその場合、どこかにメモを保存しているので
> その行の意味は理解していなくても毎回同じ内容でファイル名のみ変更して再利用するかと思います。
> なのでトレーサビリティ(再現性)という観点だと安心かもしれませんね。
>
> いわさ@松山
>
> -----Original Message-----
> From: Nobusuke Iwasaki <wata909 @ gmail.com>
> Sent: Thursday, January 20, 2022 6:24 AM
> To: shimada yasu <yas.shimada35 @ gmail.com>
> Cc: 岩佐 <y-iwasa @ limecon.co.jp>; osgeojapan-discuss <osgeojapan-discuss @ lists.osgeo.org>
> Subject: Re: [OSGeoJapan-discuss] GRASSについて
>
> いわささま
>
>> 個人的にはGRASSでの作業はGRASSで行いファイルが完成したらQGISで開くというのが初心者にとっては一番いいかと感じました。
> そうですね。QGISとGRASSの連携が完全にうまくいくとは言えないので、やはり、GRASSでできることは、GRASSでしたほうが手間がかからないのでは、と思います。
>
>> しまださま
> 昔はCYGWINでしたね!
> 今は、Windowsのコマンドプロンプト(cmd.exe)から立ち上がります。
> shell scriptであれ、batファイルであれ、定型処理を簡単にできるの(かつ、再現できる)のがいいなと思います。
>
> QGIS等でもPythonで書けばいいと思うのですが、簡単な定型処理であれば、Shell Scriptのほうが、学習コストが低いかなと、思います。
> あと、QGISのprocessingも、直接コマンドからいじれるようになったようです。Windowsの場合、OSGeo4Wで
>
>   qgis_process list
>
> と入力すると、使用可能なコマンドのリストが出てきます。
> たくさん出てきすぎて、どれを使ったらいいのかわからないのですが(汗
>
> 以下、参考URLです。
> https://bnhr.xyz/2020/07/17/run-qgis-algorithms-command-line.html
> https://www.bing.com/videos/search?q=Deploying+QGIS+using+command+line+options.&docid=608025278916339108&mid=2F86019AFBA7BDD9CF1C2F86019AFBA7BDD9CF1C&view=detail&FORM=VIRE
>
> 2022年1月19日(水) 19:48 shimada yasu <yas.shimada35 @ gmail.com>:
>> 些末な投稿ですので、無視してもらって構いません
>>
>> 以前grassは、単独で使用してました(cygwin上ですね)
>>
>> 筑波大学の奈佐原先生のHPにgrassの使い方、講義資料が公開されていて、それで勉強しました
>>
>> 先生はlinuxユーザーのようで、シェルスクリプトの重要性を強調してました
>>
>> grassもシェルスクリプトで走るところが「推し」ポイント(もちろん無料も)のようで、スクリプトで作成するから、時間が経過しても再現可能とし
>> 、それをトレーサビリティーと呼んでいたと記憶しています
>>
>> それ(スクリプトで記述する)は、初心者の私でよとっても賛同できます、はい
>>
>> 最近のgrassは、Windowsに対応し、スクリプトもPythonでも可になったと聞いていますが、詳細は分かりません
>>
>>
>> しまだ
>>
>> 2022年1月19日(水) 11:00 <y-iwasa @ limecon.co.jp>:
>>> いわさきさま
>>>
>>> 回答ありがとうございます。
>>>
>>> 1(GRASSでないとできない機能)
>>> 機能的にはQGISでもできそうですがQGIS単独だとうまくできなさそうですね。
>>> プロセッシングツールの場合だと多分ロケーションとマップセットがないというところで止まってしまうかと。
>>> またそこはGRASSで作成してあっても複数ある場合どこまで対応できるのが疑問ですね。
>>> プラグインだと新規にロケーションやマップセット等の設定はできるのですがこちらは入力ラスタの選択ができませんでした。
>>> ここもGRASSで一度処理した履歴があると候補に出てくるようです。
>>> 個人的にはGRASSでの作業はGRASSで行いファイルが完成したらQGISで開くというのが初心者にとっては一番いいかと感じました。
>>>
>>> あとは実際の業務で使う用途が増えてくるとQGIS側でできた方が便利なのかもしれませんね。
>>>
>>> 2(入力ラスタ選択)
>>> 処理対象が一つならばどれでもいいという感じですね。
>>>
>>> いわさ@松山
>>>
>>> -----Original Message-----
>>> From: Nobusuke Iwasaki <wata909 @ gmail.com>
>>> Sent: Wednesday, January 19, 2022 9:24 AM
>>> To: 岩佐 <y-iwasa @ limecon.co.jp>
>>> Cc: osgeojapan-discuss <osgeojapan-discuss @ lists.osgeo.org>
>>> Subject: Re: [OSGeoJapan-discuss] GRASSについて
>>>
>>> いわささま
>>>
>>>> 1.GRASSでないとできない機能(QGISではできない)
>>>> GRASSのメニューバー画面を見るといろいろな解析機能があるように感じましたがここら辺の機能はQGISでもあるものなのでしょうか?
>>>> もしくはGRASSしかない便利な機能があるのでしょうか。コマンドプロンプトと異なりGUI操作なのでGRASSだととっつきやすいかと感じました。
>>> ちゃんと精査したわけではありませんが、GRASSでしか使えない機能というのはおそらくなくて、QGISのprocessing
>>> から使うことはできると思います
>>> 私の場合、たくさんのファイル、もしくは定型処理をするときにGRASSを使うことが多いです。(コマンドラインからの使用になってしまいますが…)
>>>
>>>> 2.計算時の対象となる入力ラスタはどれを選べばいいのでしょうか?
>>> これはうろ覚えなのですがGRASSの場合、データの頭文字に数字を使えなかったような気がします。
>>> なので、503452の頭が自動的にxに変更されたのではないかと思います。
>>> x03452のうちどれを使うかですが、どちらでも同じだと思います。
>>> mapsetがたくさんあって、表示しているデータがわからない時のために「マップディスプレイ」という項目があるのだと思いますが、参照してあるデータは同じだと思います。
>>>
>>> こんなかんじで、よろしいでしょうか?
>>>
>>> 2022年1月18日(火) 13:22 <y-iwasa @ limecon.co.jp>:
>>>> いわさき様
>>>>
>>>> お世話になります。
>>>> 何となくGRASS自体がわかってきました。
>>>>
>>>> ”地理参照ファイルを選択”してGRASS単独でもできるようになりました。
>>>> ロケーションとマップセットについてはいわさき様が書かれたようにパーマネントでも十分ですね。
>>>>
>>>> 2点追加質問があります。
>>>> 1.GRASSでないとできない機能(QGISではできない)
>>>> GRASSのメニューバー画面を見るといろいろな解析機能があるように感じましたがここら辺の機能はQGISでもあるものなのでしょうか?
>>>> もしくはGRASSしかない便利な機能があるのでしょうか。コマンドプロンプトと異なりGUI操作なのでGRASSだととっつきやすいかと感じました。
>>>>
>>>> 2.計算時の対象となる入力ラスタはどれを選べばいいのでしょうか?
>>>> 添付ファイルのようにマップセットで2つマップディスプレイで1つが候補として出てくるのですが
>>>> どれを選べばいいのでしょうか?マップセットの2つはオリジナルとインポート時にファイ名の最初がxになっているものなので
>>>> 同じかと思っています。ただし3つのうちどれを選べばいいのかわからない状態です。
>>>>
>>>> まずはGRASS単独で使えるようになり、事前準備段階でのロケーションとマップセットの関係については理解できました。
>>>>
>>>> いわさ@松山
>>>>
>>>> -----Original Message-----
>>>> From: Nobusuke Iwasaki <wata909 @ gmail.com>
>>>> Sent: Tuesday, January 18, 2022 10:12 AM
>>>> To: 岩佐 <y-iwasa @ limecon.co.jp>
>>>> Cc: osgeojapan-discuss <osgeojapan-discuss @ lists.osgeo.org>
>>>> Subject: Re: [OSGeoJapan-discuss] GRASSについて
>>>>
>>>> いわささん、みなさん
>>>>
>>>> いわさきです。
>>>>
>>>>> ロケーションとマップセットについての考え方について教えて頂ければ助かります。
>>>>>
>>>>>  ロケーション:shikoku(CRSは6672)
>>>>>  マップセット:komi(現場名)
>>>>>
>>>>> 自分の場合、作業対象は四国内=CRS6672なのでロケーションは1つとしマップセットを現場ごとに増やしていくような
>>>> 感じでやっていくのがいいのでしょうか?
>>>>
>>>> 実は、その辺のやり方は人それぞれだったりします。
>>>> かなり古い資料ですが、私の作った資料が以下になるのですが、私の場合は以下のように使い分けています。
>>>> https://www.slideshare.net/wata909/20gisgrass/11
>>>>
>>>> ・LOCATION:座標系、分析対象範囲、ラスタのデフォルト解像度を定義。異なる座標系は混在できない
>>>> ・MAPSET:実際の分析や処理の対象や、データの型によってmapsetを作り分ける。たとえば、DEMの処理結果、ベクタデータ等々
>>>>   ただし、共通して使うものは、PARMANENTに格納する
>>>>
>>>> ということで、もし自分がやればですが、以下のようになるかと思います。
>>>> ・LOCATION:komi(現場名)
>>>> ・MAPSET:PERMANENT(もとになるオルソ画像、DSM)
>>>>        mapset1(例えば、DMSを処理したデータ)
>>>>        mapset2(例えば、ベクタデータ)
>>>> r.fillnullsをするだけであれば、個別のマップセットは作らず、PERMANENTだけでもいいと思います。
>>>>
>>>> また作業に詰まっている部分の「3.GRASSだと計算対象範囲指定方法」ですが、locationを作るときに、処理対象とするDSMを読みこんでつくると、自動的にそのDSMの範囲、解像度が処理対象として設定されます。
>>>> 以下のスライドや、GRASSビギナーズマニュアルの「ロケーションとマップセットの作成」でいうと、「新しいロケーションの作成方法を選択」のところで、「地理参照ファイルを選択」にチェックを入れて、処理対象となるDSMを選択すると、DSMにあわせてロケーションが作成されると思います。
>>>> https://www.slideshare.net/wata909/20gisgrass/25
>>>> https://gis-oer.github.io/gitbook/book/materials/GRASS/GRASS.html
>>>> https://gis-oer.github.io/gitbook/book/materials/GRASS/pic/grasspic_5.
>>>> png
>>>>
>>>> ここで、GRASSのメニューの「設定」→「Computatinal region」→「Show current region [g.region -p]」を選んでいただくと、計算範囲や解像度の数値が表示されます。
>>>> 以上の手順、ざっくりキャプチャしたものを共有します。
>>>> https://www.dropbox.com/s/ohbdtiehy4fpedk/GRASS_Mapset%E4%BD%9C%E6%
>>>> 88%
>>>> 90.docx?dl=0
>>>>
>>>> QGISから扱うのは、ちょっと時間がないので、また機会があれば…
>>>>
>>>> 2022年1月17日(月) 15:18 <y-iwasa @ limecon.co.jp>:
>>>>> 岩崎様
>>>>>
>>>>> やっと時間ができてGRASSを今日試したのですが午前中一杯かけてどうにかできるように
>>>>> なったのですがロケーションとマップセットについての考え方について教えて頂ければ助かります。
>>>>>
>>>>> ロケーション:shikoku(CRSは6672)
>>>>> マップセット:komi(現場名)
>>>>>
>>>>> 自分の場合、作業対象は四国内=CRS6672なのでロケーションは1つとしマップセットを現場ごとに増やしていくような
>>>>> 感じでやっていくのがいいのでしょうか?
>>>>>
>>>>> なおGRASSでの処理計算ですが以下のようか通りでした(QGISは3.22です)
>>>>> ①プロセッシングツールで"r.fillnulls"を実行:エラー
>>>>> ②プラグインのGRASSで実行:入力ファイル(処理対象のファイル)が候補に出ないも一度GRASSの単独プログラム側でラスタファイル
>>>>> を対象
>>>>> マッ
>>>>> プセットで
>>>>> 表示させるとその後はQGIS側で変換対象で表示されます、また対象範囲はQGISのプラグイン側で設定しないとできませんでした。
>>>>> ③GRASS単独:範囲指定方法がわからず、そのまま実行するとエラーで終了、ただし事前にQGISのプラグインで範囲が指定されいると処
>>>>> 理計算
>>>>> は問
>>>>> 題なく完了
>>>>>
>>>>> なので現在の状況では②のQGISプラグインだとどうにかGRASSでの計算ができた状況です。ただし正しい手順での作業かは自信がありません。
>>>>> 何となく何も設定をしていない状態で①(プロセッシングツール)だとうまく動かないかと、結局②か③になるのですが②プラグインだと入力対
>>>>> 象ファ
>>>>> イル
>>>>> がQGIS上では登録できない
>>>>> ③GRASSだと計算対象範囲指定方法がわからず結局QGISとGRASSを行ったり来たりという状況です。
>>>>>
>>>>> 何かやり方が間違っていただアドバイスを頂ければ助かります。
>>>>>
>>>>> いわさ@松山
>>>>>
>>>>> -----Original Message-----
>>>>> From: IWASAKI Nobusuke <niwasaki @ naro.affrc.go.jp>
>>>>> Sent: Thursday, January 6, 2022 8:12 PM
>>>>> To: y-iwasa @ limecon.co.jp; osgeojapan-discuss @ lists.osgeo.org
>>>>> Subject: Re: [OSGeoJapan-discuss] DEMデータ加工(水面部の処理修正)
>>>>>
>>>>> いわささま
>>>>>
>>>>> GRASSは単独でも使う事が出来ますし、QGISのprocessingから使うこともできます。
>>>>> 個人的にはGRASS単独で使う方が安定しているからなと思うのですが、数枚のデータを処理するだけであれば、QGISのprocesingから使った方が簡単でいいと思います。
>>>>>
>>>>> 以下、参考のサイトです
>>>>> GRASSビギナーズマニュアル
>>>>> https://gis-oer.github.io/gitbook/book/materials/GRASS/GRASS.html
>>>>>
>>>>> チュートリアル: GRASSツールを用いて河川と集水域の境界線を求める
>>>>> (使っている処理は違いますが、QGISからやる場合、同じ手順でできると思います)
>>>>> https://courses.gisopencourseware.org/mod/book/tool/print/index.php?
>>>>> id
>>>>> =213
>>>>>
>>>>>
>>>>> On 2022/01/06 19:57, y-iwasa @ limecon.co.jp wrote:
>>>>>> 補足です。
>>>>>>
>>>>>> GRASSは別のプログラムなんですね。
>>>>>> まずはGRASSの設定をしないといけないようなのでまずはGRASSの設定&使い方を
>>>>>> 調べてから週末にこの機能について確認させて頂きます。
>>>>>>
>>>>>> いわさ@松山
>>>>>>
>>>>>> -----Original Message-----
>>>>>> From: y-iwasa @ limecon.co.jp <y-iwasa @ limecon.co.jp>
>>>>>> Sent: Thursday, January 6, 2022 7:52 PM
>>>>>> To: 'IWASAKI Nobusuke' <niwasaki @ naro.affrc.go.jp>;
>>>>>> 'osgeojapan-discuss @ lists.osgeo.org'
>>>>>> <osgeojapan-discuss @ lists.osgeo.org>
>>>>>> Subject: RE: [OSGeoJapan-discuss] DEMデータ加工(水面部の処理修正)
>>>>>>
>>>>>> 岩崎様
>>>>>>
>>>>>> 追加コメントありがとうございます。
>>>>>> 早速確認させて頂きます。
>>>>>>
>>>>>> ただしGRASSのコマンド自体使うのが初めてなのでちょっと不安です。
>>>>>>
>>>>>> いわさ@松山
>>>>>>
>>>>>> -----Original Message-----
>>>>>> From: OSGeoJapan-discuss
>>>>>> <osgeojapan-discuss-bounces @ lists.osgeo.org>
>>>>>> On Behalf Of IWASAKI Nobusuke
>>>>>> Sent: Thursday, January 6, 2022 11:11 AM
>>>>>> To: osgeojapan-discuss @ lists.osgeo.org
>>>>>> Subject: Re: [OSGeoJapan-discuss] DEMデータ加工(水面部の処理修正)
>>>>>>
>>>>>> いわささま、田鎖さま
>>>>>>
>>>>>> 岩崎です。
>>>>>> GRASSのコマンドで、r.fillnullsというのがあります。
>>>>>> https://grass.osgeo.org/grass78/manuals/r.fillnulls.html
>>>>>>
>>>>>> みんなの自動翻訳で翻訳させると、
>>>>>> r.fillnullsは、入力ラスターマップのNULLピクセル(データ領域なし)を塗りつぶし、塗りつぶされたデータを新しい出力ラスターマップに格納します。塗りつぶし領域は、v.surf.rst正則化スプライン補間(方法=rst)またはr.resamp.bsplineキュービックまたは線形スプライン補間(Tykhonov正則化)を使用して、非データ領域境界バッファから補間されます。
>>>>>>
>>>>>> という感じです。
>>>>>> GRASSの機能なので走らせるまでがちょっと手間かもしれませんが、参考までに。
>>>>>>
>>>>>>
>>>>>> On 2022/01/05 6:31, y-iwasa @ limecon.co.jp wrote:
>>>>>>> 田鎖様
>>>>>>>
>>>>>>> 実はコメントのリンクは自分も今一度調べた時に出てきて確認させて頂きました。
>>>>>>>
>>>>>>> 前半部分の補完処理はうまくいったりいかなかったりで結局挫折してしまいました。
>>>>>>>
>>>>>>> もしかすると対象が比較的川幅が細かったかもしれません。
>>>>>>>
>>>>>>> 読んでみて意外と10mDEMならば粗いですが欠測がないのでドローンでの地形把握ならば
>>>>>>>
>>>>>>> これが一番簡単なのかもしれません。
>>>>>>>
>>>>>>> 5mと10mの合成も試したのですが境目の部分がまたー9999もしくは0のままでした。
>>>>>>>
>>>>>>> いわさ@松山
>>>>>>>
>>>>>>> *From:*OSGeoJapan-discuss
>>>>>>> <osgeojapan-discuss-bounces @ lists.osgeo.org>
>>>>>>> *On Behalf Of *Junta TAGUSARI
>>>>>>> *Sent:* Tuesday, January 4, 2022 3:02 PM
>>>>>>> *To:* osgeojapan-discuss @ lists.osgeo.org
>>>>>>> *Subject:* Re: [OSGeoJapan-discuss] DEMデータ加工(水面部の処理修正)
>>>>>>>
>>>>>>> いわさ様
>>>>>>>
>>>>>>> はじめまして,田鎖と申します。
>>>>>>> 地理情報が専門というわけではないのですが,業務で必要に迫られてGISを触っている者です。
>>>>>>>
>>>>>>> 水面データはー9999で入っていますが,これは「欠測」として扱い,周囲の値を使って補間するのが良いのではないかと思います。
>>>>>>> 過去によく似たやりとりが行われているのを見つけましたのでご参照ください。
>>>>>>>
>>>>>>> https://groups.google.com/g/qgisshitumon01/c/mTnnGJGk0-U?pli=1
>>>>>>> <https://groups.google.com/g/qgisshitumon01/c/mTnnGJGk0-U?pli=
>>>>>>> 1>
>>>>>>>
>>>>>>> ちなみに,このやり取りの中でも触れられていますが,海面の高さまで勝手に補間されるので注意が必要です。
>>>>>>> たとえば,海が陸地で挟まれている湾のような地形で,かつ海に接している陸地の標高がそこそこ大きいと,湾内の標高が数mとかで補間されてしまいます。
>>>>>>> このような場合,目的によっては,水面を0mで設定した方が良いかもしれません。
>>>>>>>
>>>>>>> ご参考になれば幸いです。
>>>>>>>
>>>>>>> 〒060-8628 札幌市北区北13条西8丁目
>>>>>>>
>>>>>>> 北海道大学大学院工学研究院 環境工学部門
>>>>>>>
>>>>>>> 地域環境研究室
>>>>>>>
>>>>>>> 田鎖順太 Junta Tagusari
>>>>>>>
>>>>>>> Tel&Fax: 011-706-6872
>>>>>>>
>>>>>>> Email: j.tagusari @ eng.hokudai.ac.jp
>>>>>>> <mailto:j.tagusari @ eng.hokudai.ac.jp>
>>>>>>>
>>>>>>> On 2022/01/03 9:00, y-iwasa @ limecon.co.jp <mailto:y-iwasa @ limecon.co.jp> wrote:
>>>>>>>
>>>>>>>       いわさ@松山です。
>>>>>>>
>>>>>>>       本来の本年最初の投稿をさせて頂きます。
>>>>>>>
>>>>>>>       仕事はじめは1/5なのでいきなり現場観測なので...
>>>>>>>
>>>>>>>       実は最近ドローンでの撮影する頻度が増えてきて事前に可能ならばDEMデータを地理院からダウンロードして
>>>>>>>
>>>>>>>       確認しているのですが川などの水面部分をできれば川辺の高さに表示できるとイメージがわきやすいのですが。
>>>>>>>
>>>>>>>       いまだと単純に水面部が“-9999”にしているので渓谷に見えてしまいます...
>>>>>>>
>>>>>>>       基盤地図情報の変換ソフトでは水面は0かー9999に指定できるのですが例えば後者(-9999)に指定すれば
>>>>>>>
>>>>>>>       ラスタ計算機でできるのかなと思っていますがその方法がわからず。
>>>>>>>
>>>>>>>       自分で思いついたのはその場所の標高あたりに適当に設定してif文で“-9999”の場合“標高試算値”に置き換えですが
>>>>>>>
>>>>>>>       
>>>>>>> これだと対象画像が同じ高さに修正されています。できれば水面と陸地の境の標高値にスムーズにかつ自動的に合わせてくれると
>>>>>>>
>>>>>>>       いいのですがこのような機能はありませんでしょうか?
>>>>>>>
>>>>>>>       もしかしたらFUSIONなどであるかもしれないのですがご存知の方がいらっしゃいましたら教えて頂けると助かります。
>>>>>>>
>>>>>>>       いわさ@松山
>>>>>>>
>>>>>>>
>>>>>>>
>>>>>>>       _______________________________________________
>>>>>>>
>>>>>>>       OSGeoJapan-discuss mailing list
>>>>>>>
>>>>>>>       OSGeoJapan-discuss @ lists.osgeo.org
>>>>>>> <mailto:OSGeoJapan-discuss @ lists.osgeo.org>
>>>>>>>
>>>>>>>       
>>>>>>> https://lists.osgeo.org/mailman/listinfo/osgeojapan-discuss
>>>>>>> <https://lists.osgeo.org/mailman/listinfo/osgeojapan-discuss>
>>>>>>>
>>>>>>>
>>>>>>> _______________________________________________
>>>>>>> OSGeoJapan-discuss mailing list
>>>>>>> OSGeoJapan-discuss @ lists.osgeo.org
>>>>>>> https://lists.osgeo.org/mailman/listinfo/osgeojapan-discuss
>>>>>> --
>>>>>> ☆通常勤務中です☆
>>>>>> #staypositive
>>>>>>
>>>>>> ---
>>>>>> Nobusuke IWASAKI
>>>>>> Representative, OSGeo Japan Chapter
>>>>>>
>>>>>> Institute for Agro-Environmental Sciences, NARO 3-1-1,
>>>>>> Kannondai, Tsukuba-shi, Ibaraki-ken 305-8604, JAPAN Tel / Fax
>>>>>> +81298388227 _______________________________________________
>>>>>> OSGeoJapan-discuss mailing list
>>>>>> OSGeoJapan-discuss @ lists.osgeo.org
>>>>>> https://lists.osgeo.org/mailman/listinfo/osgeojapan-discuss
>>>>> --
>>>>> ☆通常勤務中です☆
>>>>> #staypositive
>>>>>
>>>>> ---
>>>>> Nobusuke IWASAKI
>>>>> Representative, OSGeo Japan Chapter
>>>>>
>>>>> Institute for Agro-Environmental Sciences, NARO 3-1-1, Kannondai,
>>>>> Tsukuba-shi, Ibaraki-ken 305-8604, JAPAN Tel / Fax +81298388227
>>>>> _______________________________________________
>>>>> OSGeoJapan-discuss mailing list
>>>>> OSGeoJapan-discuss @ lists.osgeo.org
>>>>> https://lists.osgeo.org/mailman/listinfo/osgeojapan-discuss
>>>>
>>>>
>>>> --
>>>> 岩崎 亘典
>>>
>>>
>>> --
>>> 岩崎 亘典
>>> _______________________________________________
>>> OSGeoJapan-discuss mailing list
>>> OSGeoJapan-discuss @ lists.osgeo.org
>>> https://lists.osgeo.org/mailman/listinfo/osgeojapan-discuss
>
>
> --
> 岩崎 亘典
> _______________________________________________
> OSGeoJapan-discuss mailing list
> OSGeoJapan-discuss @ lists.osgeo.org
> https://lists.osgeo.org/mailman/listinfo/osgeojapan-discuss
-------------- next part --------------
HTML$B$NE:IU%U%!%$%k$rJ]4I$7$^$7$?(B...
URL: <http://lists.osgeo.org/pipermail/osgeojapan-discuss/attachments/20220120/d322aff9/attachment-0001.html>
-------------- next part --------------
---
title: "Define sokushin area"
author: "Junta Tagusari"
output: 
  html_document:
    toc: false
    toc_depth: 2
    toc_float:
      collapsed: false
    number_sections: true
editor_options: 
  chunk_output_type: console
---


# なにがしたいか
 - 洋上風車促進区域の出力
 - 「由利本荘」と,「八峰町及び能代市沖」.
 
# 設定

## ライブラリ
```{r loadlib, message = F}
library(tidyverse)
library(sf)
library(leaflet)
```

## セッション情報
```{r sessioninfo}
sessionInfo()
Sys.time()
```

## プロット設定

```{r}
setThemeCustomStyle <- function(fontsize = 10, fontfamily = "serif"){
  theme_set(
    theme_bw() + 
      theme(
        panel.grid.minor  = element_blank(),
        text              = element_text(family = fontfamily),
        legend.title      = element_text(size = fontsize),
        legend.text       = element_text(size = fontsize),
        axis.ticks.length = unit(-fontsize / 5,"points"),
        axis.text.x       = element_text(size = fontsize, margin = margin(t = fontsize * 0.6)),
        axis.text.y       = element_text(size = fontsize, margin = margin(r = fontsize * 0.4)),
        axis.title        = element_text(size = fontsize),
        strip.placement   = "outside",
        strip.background  = element_blank()
      )
  )
  
}

setThemeCustomStyle(fontfamily = "Yu Mincho")
```
# 由利本荘市沖の促進区域

## 線で定義されている

```{r}
sokushin_coords_north <- tribble(
  ~ctg,    ~lat,                 ~lng,
  "north", 39 + 35/60 + 52/3600, 140 +  3/60 + 42/3600,
  "north", 39 + 35/60 + 52/3600, 140 +  0/60 + 50/3600,
  "north", 39 + 31/60 +  7/3600, 140 +  0/60 +  5/3600,
  "north", 39 + 27/60 + 25/3600, 139 + 58/60 + 52/3600,
  "north", 39 + 27/60 +  5/3600, 140 +  2/60 +  4/3600
)

sokushin_coords_south <- tribble(
  ~ctg,    ~lat,                 ~lng,
  "south", 39 + 27/60 +  5/3600, 140 +  2/60 +  4/3600,
  "south", 39 + 27/60 + 25/3600, 139 + 58/60 + 52/3600,
  "south", 39 + 19/60 + 34/3600, 139 + 56/60 + 17/3600,
  "south", 39 + 18/60 + 43/3600, 139 + 59/60 +  4/3600
)

sokushin_lines <- list(
  st_linestring(matrix(c(sokushin_coords_north$lng, sokushin_coords_north$lat), ncol = 2), dim = "XY"),
  st_linestring(matrix(c(sokushin_coords_south$lng, sokushin_coords_south$lat), ncol = 2), dim = "XY")
) %>% 
  st_as_sfc(crs = 6668) %>% 
  st_as_sf() %>% 
  dplyr::mutate(label = c("north", "south")) %>% 
  st_transform(4326)

leaflet(sokushin_lines) %>% 
  addTiles() %>% 
  addPolylines(label = ~label)
```


南北まとめてしまう

```{r}
sokushin_coords_join <- bind_rows(
  sokushin_coords_north[1:4,],
  sokushin_coords_south[3:4,]
) 
  
sokushin_line_join <- list(
  st_linestring(matrix(c(sokushin_coords_join$lng, sokushin_coords_join$lat), ncol = 2), dim = "XY")
) %>% 
  st_as_sfc(crs = 6668) %>%
  st_transform(4326)

leaflet(sokushin_line_join) %>% 
  addTiles() %>% 
  addPolylines()
```

## ポリゴンにする

いろいろ試行錯誤した結果を含む.

```{r}
sf_coast <- read_sf("geo/akita-wt-hrisk.gpkg", layer = "natural-coastline") %>% 
  st_crop(st_bbox(sokushin_lines)) %>% 
  .[c(1,3,2),] %>% 
  st_coordinates() %>% 
  .[,1:2] %>% 
  st_linestring(dim = "XY") %>% 
  list() %>% 
  st_as_sfc(crs = 4326) 


leaflet(sf_coast) %>% 
  addTiles() %>% 
  addPolylines(color = "red") %>% 
  addPolylines(data = sokushin_line_join)

sf_coast_to_plg <- lwgeom::st_split(sf_coast, sokushin_line_join) %>% 
  st_collection_extract("LINESTRING") %>% 
  st_as_sf() %>% 
  dplyr::mutate(id = row_number(), len = as.double(st_length(.))) %>% 
  dplyr::filter(len > 1e3)

sokushin_line_to_plg <- lwgeom::st_split(sokushin_line_join, sf_coast) %>% 
  st_collection_extract("LINESTRING") %>% 
  st_as_sf() %>% 
  dplyr::mutate(id = row_number(), len = as.double(st_length(.))) %>% 
  dplyr::filter(len > 1e3)

leaflet(sf_coast_to_plg) %>% 
  addTiles() %>% 
  addPolylines(color = "red") %>%
  addPolylines(data = sokushin_line_to_plg)

sokushin_plg <- rbind(
  st_coordinates(sokushin_line_to_plg),
  st_coordinates(sf_coast_to_plg) %>% 
    .[nrow(.):1,]
) %>% 
  .[,1:2] %>%
  list() %>% 
  st_polygon() %>% 
  list() %>% 
  st_as_sfc(crs=4326)

leaflet(sokushin_plg) %>% 
  addTiles() %>% 
  addPolygons()

# write_sf(sokushin_plg, "geo/akita-wt-hrisk.gpkg", layer = "yurihonjo-sokushin")
```


# 「八峰町及び能代市沖」の促進区域

## 線で定義されている

```{r}
sokushin_coords <- tribble(
  ~lat,                 ~lng,
  40 + 17/60 + 50/3600, 140 +  1/60 + 22/3600,
  40 + 17/60 + 52/3600, 140 +  0/60 + 32/3600,
  40 + 18/60 + 28/3600, 139 + 59/60 +  7/3600,
  40 + 18/60 + 59/3600, 139 + 58/60 + 39/3600,
  40 + 15/60 + 51/3600, 139 + 58/60 +  8/3600,
  40 + 12/60 + 12/3600, 139 + 56/60 + 53/3600,
  40 + 11/60 + 44/3600, 139 + 59/60 + 55/3600
)

sokushin_lines <- list(
  st_linestring(matrix(c(sokushin_coords$lng, sokushin_coords$lat), ncol = 2), dim = "XY")
) %>% 
  st_as_sfc(crs = 6668) %>% 
  st_as_sf() %>% 
  st_transform(4326)

leaflet(sokushin_lines) %>% 
  addTiles() %>% 
  addPolylines()
```

## ポリゴンにする

いろいろ試行錯誤した結果を含む.

```{r}
sf_coast <- read_sf("geo/akita-wt-hrisk.gpkg", layer = "natural-coastline") %>% 
  st_crop(st_bbox(sokushin_lines)) %>% 
  st_cast("LINESTRING")


leaflet(sf_coast) %>% 
  addTiles() %>% 
  addPolylines(color = "red") %>% 
  addPolylines(data = sokushin_lines)

sf_coast_to_plg <- lwgeom::st_split(sf_coast, sokushin_lines) %>% 
  st_collection_extract("LINESTRING") %>% 
  st_as_sf() %>% 
  dplyr::mutate(id = row_number(), len = as.double(st_length(.))) %>% 
  dplyr::filter(len > 1e3)

sokushin_line_to_plg <- lwgeom::st_split(sokushin_lines, sf_coast) %>% 
  st_collection_extract("LINESTRING") %>% 
  st_as_sf() %>% 
  dplyr::mutate(id = row_number(), len = as.double(st_length(.))) %>% 
  dplyr::filter(len > 1e3)

leaflet(sf_coast_to_plg) %>% 
  addTiles() %>% 
  addPolylines(color = "red") %>%
  addPolylines(data = sokushin_line_to_plg)

sokushin_plg <- rbind(
  st_coordinates(sokushin_line_to_plg),
  st_coordinates(sf_coast_to_plg) %>% 
    .[nrow(.):1,]
) %>% 
  .[,1:2] %>%
  list() %>% 
  st_polygon() %>% 
  list() %>% 
  st_as_sfc(crs=4326)

leaflet(sokushin_plg) %>% 
  addTiles() %>% 
  addPolygons()

# write_sf(sokushin_line_to_plg, "geo/akita-wt.gpkg", layer = "noshiro-sokushin-lin")
```


OSGeoJapan-discuss メーリングリストの案内