GIS奮闘記

元GISエンジニアの技術紹介ブログ。主にPythonを使用。

スポンサーリンク

ArcGIS for Personal Useを使ってみよう!

さて、本日はArcGIS for Personal Useについて書いてみようと思います。

おそらくこのブログを読んでいただいている方の大半はArcGISのことはご存知かと思いますが、ご存知ではないという方は以下のエントリーに概要を記載していますので、興味がありましたら読んでみてください。
www.gis-py.com

ArcGISは商用のGISエンジンです。学校などには無償で提供されていますが、基本的には使用するにはお金がかかります。金額も個人で購入するには少しお高いかもしれません。しかし、なんとArcGIS for Personal Useを購入すると以下のソフトウェアが1 年間 18,000 円(税別)で利用することができます(非商用の個人利用に限定)。これはお買い得だと思います。

ArcGIS for Personal Use の購入方法

以下サイトで購入することができます。
ArcGIS for Personal Use | ESRIジャパン

ArcGIS for Personal Use に含まれるソフトウェア

ArcGIS Desktop Advanced (最新バージョン)
ArcGIS Pro
ArcGIS Online 1 指定ユーザー 100 サービス クレジット
ArcGIS 3D Analyst
ArcGIS Data Reviewer
ArcGIS Geostatistical Analyst
ArcGIS Image Analyst
ArcGIS Network Analyst
ArcGIS Publisher
ArcGIS Spatial Analyst
ArcGIS Tracking Analyst
ArcGIS Data Interoperability
ArcGIS Schematics
ArcGIS Workflow Manager

うーむ、豪華なラインナップですね。

休日もArcGISを使用しよう

ArcGISは多くの機能を搭載しており、ArcGISを毎日使用している私もその機能のほんの一部しか知りません。業務で使うソフトウェアは決まっていますし、使用する機能の範囲も決まっています。ですので、それを続けているだけでは自分の知識の幅を広げることができません。ArcGISは様々な業種で使われており、その使われ方も様々です。自分のスキルレベルを上げるには自分の使ったことがない機能を積極的に使う必要があるかと思います。自宅でArcGISを使うことができればそれも可能になるかと思います。

ArcGISの進化

ArcGISは日々進化しています。特に最近はArcGIS OnlineやArcGIS Proへのシフトがどんどん進んでいますし、その動きがますます加速するかと思います。しかし、業務上の都合などにより、こういった最新のソフトウェアへのキャッチアップができていない方もまだまだたくさんいらっしゃるかと思います。私も最近、ようやくArcGIS Proを使い始めました。Proと言っても基本的な考え方は今までと変わっていませんが、ArcMapとは劇的にUIが異なり、四苦八苦している最中です。その一方で「こういった考え方もあるのか」と思ったりすることで、自分の考え方の幅も広がっているような気がします。ITの世界は技術革新が早く、最新技術も数年で陳腐化してしまうことも珍しくありません。ArcGISの進化にあわせて自分も進化していかないといけないと感じています。


今日はプログラムの紹介とかは一切なく、ArcGISの宣伝のようになってしまいましたね。。もしかしたらしばらくはArcGIS関連のことを書くことが多くなるかもしれませんが、できるだけ偏らないように色々なことに挑戦してみようと思います。「こういったことを書いてほしい」、「こういったことはできないのか?」などありましたら、ぜひコメントを残していただければと思います。そちらにもチャレンジしてみようと思います。



ArcPyレシピ集⑥ ~重複した属性値を抽出~

さて、本日はArcPyについて書いてみようと思います。

以前、以下のエントリーで空間的に重複したフィーチャの抽出方法を紹介しました。

www.gis-py.com

今回は空間的ではなく、重複した属性、つまり、同カラムに同じ属性値を持つレコードの抽出方法を紹介します。

ArcPyとは

言わずもがなですが、ArcGISのデータを扱うためのPythonライブラリです。ArcPyに関しては当ブログで色々紹介していますので、興味のある方は以下エントリーを読んでみてください。

www.gis-py.com

www.gis-py.com

環境

Windows10 64bit
ArcGIS10.4.1 for Desktop
Python2.7.3

サンプルデータ

以下のようなデータを用意しました。testというカラムに重複した属性値をいくつか設定しています。
f:id:sanvarie:20190312105939p:plain

f:id:sanvarie:20190312110003p:plain

f:id:sanvarie:20190312112111p:plain

1,2,7,9が重複してますね!

サンプルコード

上記データの中にある重複した属性値を抽出するサンプルコードです。

# -*- coding: utf-8 -*-
import arcpy
from collections import Counter

arcpy.env.workspace = ur"C:\ArcPySample\Duplicate.gdb"

feature_class = "test"
column = "test"
c = []

for f in arcpy.da.SearchCursor(feature_class, column):
    c.append(f[0])

print [key for key,val in Counter(c).items() if val > 1]


結果を見ると、重複している属性値を抽出できたことが確認できました。ものすごい簡単ですね。
f:id:sanvarie:20190312112803p:plain

こんなのどういう時に使うの?と思われる方もいらっしゃるかもしれませんが、実際私は結構使います。例えば、キー項目として使っているカラムのデータに重複が発生してしまっている場合、こういうことをして対象データを抽出しています(そもそもキーカラムに重複を発生させるなよ、という話ですが)。ArcPyに興味がある方で「こんなことできないか?」という疑問がある方はぜひコメントを残してください。私にできることであればお答えできればと思います。本日は以上です。

Pythonで国土数値情報のWeb APIを使って全データを一括でダウンロードする方法

さて、今回も国土数値情報のWeb APIについて書いてみようと思います。先日、本件に関するエントリーを書きましたが、この時は一部データのダウンロードしかしませんでした。その時に全データを一括でダウンロードできたらいいなぁと思ったので、今回はそれについて書いてみようと思います。

国土数値情報のWeb APIとは?

以下エントリーで紹介しています。興味のある方はぜひ読んでみてください。

www.gis-py.com

ダウンロード対象

国土数値情報は以下の画像のように各カテゴリごとにデータをダウンロードできます。

f:id:sanvarie:20190215154354p:plain

そして、各カテゴリごとに年度、都道府県ごとにデータをダウンロードすることができます(カテゴリによっては一つの年度のみ、または、全国単位でのみダウンロード可能)。
f:id:sanvarie:20190215154622p:plain

今回は・・・

①大分類(例:1. 国土(水・土地))
②中分類(例:<水域>)
③小分類(例:海岸線)
④年度(最新の年度)

ごとにフォルダを作成してその中にダウンロードしたデータを格納していこうと思います(カテゴリによっては中分類がなかったりします)。

また、最初はすべての年度のデータをダウンロードしようと思ったのですが、どんでもないデータ量になりそうだったのであきらめました(笑)。サンプルコードには全ての年度のデータをとってくるコードを記載しますので(コメントアウトしています)、興味がある方はぜひ試してみてください。

環境

Python2.7.3
Windows10

サンプルコード

DATA_DESTINATIONに任意のディレクトリを指定してください。そこにデータをダウンロードします。

# -*- coding: utf-8 -*-
import zipfile
import os
import os.path
import requests
from lxml import etree

#任意の保存先
DATA_DESTINATION = ur"D:\gis-py\WebAPI"

#国土数値情報の概要情報取得
URL_OVERVIEW = "http://nlftp.mlit.go.jp/ksj/api/1.0b/index.php/app/getKSJSummary.xml" \
               "?appId=ksjapibeta1&lang=J&dataformat=1"

#国土数値情報取得のURL情報取得
url_detail = "http://nlftp.mlit.go.jp/ksj/api/1.0b/index.php/app/getKSJURL.xml" \
              "?appId=ksjapibeta1&lang=J&dataformat=1"

dic = {}

def get_url(url, identifier = 0):
    """
    国土数値情報のWEBAPIを使用して各情報をXMLで取得します。
    """
    if identifier == 0:
        resp = requests.get(url, timeout=10)
    else:
        url = url + "&identifier=%s" % (identifier)
        resp = requests.get(url, timeout=10)

    return etree.fromstring(resp.content)

def download_zip():
    """
    指定したURLからZIPをダウンロードします。
    """
    year = ""
    for key, value in dic.items():
        xml = get_url((url_detail), key)

        for x in xml.iter():

            if x.tag == "year":
                year = x.text

            if x.tag == "zipFileUrl":

                if year == value[0]:
                    year = ""

                    filename = x.text.split('/')[-1]
                    r = requests.get(x.text, stream=True)
                    with open(filename, 'wb') as f:
                        for chunk in r.iter_content(chunk_size=1024):
                            if chunk:
                                f.write(chunk)
                                f.flush()

                    uncompress_zip(filename, value[1])

def uncompress_zip(filename, destination):
    """
    ZIP ファイルを指定したディレクトリに展開します。
    """
    zfile = zipfile.ZipFile(filename)
    zfile.extractall(destination)

def make_directory(xml):
    """
    DATA_DESTINATIONの下に各フォルダを作成します。
    """
    element_list = []
    new_directry_list = []

    for x in xml.iter():

        if x.tag == "identifier":
            element_list.append(x.text)
        elif x.tag == "title":
            element_list.append(x.text)
        elif x.tag == "field1":
            element_list.append(x.text)
        elif x.tag == "field2":
            element_list.append(x.text)
        elif x.tag == "areaType":
            element_list.append(x.text)
            new_directry_list.append(element_list)
            element_list = []
            continue

    for n in new_directry_list:
        new_directry = ""
        year_list = []

        new_directry = os.path.join(DATA_DESTINATION, n[2])

        if n[3] != " " and n[3] != "-": #field2が不要なデータがあるため
            new_directry = os.path.join(new_directry, n[3])

        new_directry = os.path.join(new_directry, n[1])

        for x in get_url((url_detail), n[0]).iter("year"):
            year_list.append(x.text)

        year_list = list(set(year_list))
        max_year = max(year_list)

        new_dir = os.path.join(new_directry, max_year)

        if os.path.exists(new_dir) == 0: #なかったら作る
            os.makedirs(new_dir)
        dic[n[0]] = [max_year,new_dir]

        # ↓全ての年度のデータを取得したい場合は以下を使用する
        #year_list = list(set(year_list))
        #year_list.sort()
        #
        #for l in max(year_list):

        #    new_dir = os.path.join(new_directry, l)

        #    if os.path.exists(new_dir) == 0: #なかったら作る
        #        os.makedirs(new_dir)
        #    dic[n[0]] = [l,new_dir]

if __name__ == '__main__':
    make_directory(get_url(URL_OVERVIEW, 0))
    download_zip()

結果はこのようになりました。

f:id:sanvarie:20190215165917p:plain

f:id:sanvarie:20190215170009p:plain

f:id:sanvarie:20190215170120p:plain

ちゃんとダウンロードできている感じですね。データ量なのですが、なんと83GBになりました。全年度をダウンロードしたら一体どのくらいまでいくのでしょうか(汗)

本日は以上です。国土数値情報のデータを一括でダウンロードしたいという方はぜひ試してみてください。

Pythonで国土数値情報のWeb APIを使ってデータをダウンロードしてみよう

さて、今日は国土数値情報のWeb APIを使ってデータをダウンロードしてみようと思います。

国土数値情報とは?

以下エントリーでも紹介していますが、全国の河川、避難施設、鳥獣保護区、鉄道 など様々なデータをShape形式などで提供しているWEBサイトです。

www.gis-py.com

国土数値情報のAPIとは?

GISを日頃使っている方は国土数値情報に関してはすでにご存じかと思いますが、実は国土数値情報のWEB APIというものがあります。色々なデータをダウンロードする場合、何度もダウンロード作業をする必要があったり、アンケートに答える必要があったりで、正直面倒くさいなと思った方もいるのではないかと思います。そんな時の救世主がこのAPIですまた、この!APIを利用することで、国土数値情報を組み込んだWebサイトや アプリケーションの開発を行うことが可能になります。

API仕様

http://nlftp.mlit.go.jp/ksj/api/specification_api_ksj.pdfに仕様の詳細が記載されていますので、興味のある方はぜひ読んでみてください。

ざっくり言うと以下二つの情報を取得することができます。
①国土数値情報の概要情報取得
②国土数値情報取得のURL情報取得

これだけだとよくわからないですね。。要するに・・・

①は国土数値情報では「海岸線」「ダム」などの情報を提供しているという概要的な情報を取得することができます。
f:id:sanvarie:20190212132519p:plain

②は「海岸線」などの各データの詳細情報を取得することができます。
f:id:sanvarie:20190212132559p:plain

各々の結果はXMLで取得でき、このようになります。

①国土数値情報の概要情報
f:id:sanvarie:20190212132938p:plain

②国土数値情報取得のURL情報
f:id:sanvarie:20190212133201p:plain

今回は「②国土数値情報取得のURL情報」を使用して、ファイルをダウンロードしようと思います。

環境

Python2.7.3
Windows10

サンプルコード

「国土数値情報 海岸線データ」をダウンロードするサンプルコードです。

# -*- coding: utf-8 -*-
import zipfile
import os.path
import requests
from lxml import etree

class WebAPI(object):

    def __init__(self):
        self.appId = "ksjapibeta1" #固定
        self.lang = "J"            #固定
        self.dataformat = 1        #固定
        self.identifier = "C23"    #各データ固有の番号
        self.prefCode = "1-47"     #都道府県
        self.fiscalyear = 2006     #年度
        self.dir = ur"D:\gis-py\WebAPI" #データ保存先

    def get_url(self):
        """
        国土数値情報のWEBAPIを使用して各情報をXMLで取得します。
        """

        url = "http://nlftp.mlit.go.jp/ksj/api/1.0b/index.php/app/getKSJURL.xml" \
              "?appId=%s&lang=%s&dataformat=%s&identifier=%s&prefCode=%s&" \
              "fiscalyear=%s" %(self.appId, self.lang, self.dataformat, \
                                self.identifier, self.prefCode, self.fiscalyear)

        resp = requests.get(url, timeout=10)
        tree = etree.fromstring(resp.content)

        for t in tree.iter():
            if t.tag == "zipFileUrl":
                self.download_zip(t.text)

    def download_zip(self,text):
        """
        指定したURLからZIPをダウンロードします。
        """

        filename = text.split('/')[-1]
        r = requests.get(text, stream=True)
        with open(filename, 'wb') as f:
            for chunk in r.iter_content(chunk_size=1024):
                if chunk:
                    f.write(chunk)
                    f.flush()

            self.uncompress_zip(filename)

    def uncompress_zip(self,filename):
        """
        ZIP ファイルを指定したディレクトリに展開します。
        """
        zfile = zipfile.ZipFile(filename)
        zfile.extractall(self.dir)

if __name__ == '__main__':
    w = WebAPI()
    w.get_url()

このようにファイルがダウンロードされました。
f:id:sanvarie:20190212150004p:plain

ファイルを読み込んだ結果です。
f:id:sanvarie:20190212150126p:plain

今回はデータの一部をダウンロードしてみましたが、全データをダウンロードして各データや年度ごとにフォルダ分けしてダウンロードファイルを格納する、みたいなプログラムも作ってみました。興味のある方は以下エントリーを読んでみてください。
www.gis-py.com

ArcPyレシピ集⑤ ~重複したフィーチャを抽出~

さて、本日は久しぶりにArcPyについて書いてみようと思います。

ArcPyとは

言わずもがなですが、ArcGISのデータを扱うためのPythonライブラリです。ArcPyに関しては当ブログで色々紹介していますので、興味のある方は以下エントリーを読んでみてください。

www.gis-py.com

www.gis-py.com

今回のゴール

GISデータを扱っているとたまに空間的に重複したデータに出くわすことがあります。こういった無駄なデータがあるとデータが重くなるだけでなく、システムの不具合の原因になる可能性もあります。今回はこういったデータを一括で抽出する方法を紹介しようと思います。

環境

Windows10 64bit
ArcGIS10.4.1
Python2.7.3

サンプルデータ

以下のようにポイントフィーチャクラス(平面直角座標系9系(JGD2000))を作成しました。
f:id:sanvarie:20190131100741p:plain

また、以下のようにポイントを作図しました(何点か重複している状態です)。

※「Key」カラムはキーとなるカラムと見立てて作りました。

f:id:sanvarie:20190131101841p:plain

サンプルコード

上記データの中にある重複したポイントデータを抽出するサンプルコードです。
(「Key」カラムの中身を出力します。)

# -*- coding: utf-8 -*-
import arcpy
from collections import Counter

class DuplicateGeometry(object):

    def  __init__(self):
        arcpy.env.workspace = ur"D:\gis-py\arcpy\data\SampleData.gdb"
        self.feature_class = u"SamplePoint"

    def check_dupulicate_feature(self):
        """重複しているフィーチャのFEAT_NUMをリストで返します

        arguments:
        none

        returns:
        duplicated_list -- FEAT_NUMのリスト

        """

        try:
            feature_list = []
            features_num_geo = []
            duplicated_list = []

            features = (row for row in arcpy.da.SearchCursor(self.feature_class, "SHAPE@XY"))
            feature_list = self.extract_data(features)
            features_num_geo = (row for row in arcpy.da.SearchCursor(self.feature_class, ["Key","SHAPE@XY"]))

            for ps in features_num_geo:
                for pl in feature_list:
                    if pl[0] == ps[1][0] and pl[1] == ps[1][1]:

                        duplicated_list.append(ps[0]) #重心XYが一致した場合FEAT_NUMを格納

        except Exception:
            return None
        else:
            return duplicated_list

    def extract_data(self,features):
        #重心XYが重複しているデータを抽出
        c = Counter(tuple(items) for items in features for items in items)
        return [items for items, count in c.most_common() if count > 1]

if __name__ == '__main__':
    d = DuplicateGeometry()
    duplicate_geometry = d.check_dupulicate_feature()
    print duplicate_geometry

実行すると以下のような結果になりました。
f:id:sanvarie:20190131102140p:plain

データを確認してみると、以下画像のようにポイントが重複していることがわかります。
f:id:sanvarie:20190131102317p:plain

今回抽出した他のデータも同様の状態ですが、同じような画像を貼り付けるだけになってしまうので、今回は割愛します。

このソースだとデータが多い時に少し時間がかかってしまうのが難ですね。もう少し効率の良い方法があればそれも今後紹介できればと思います。

今回はArcPyを使いましたが、ArcMapでも同様のことはできるかと思います。

例えば・・・
①対象のフィーチャクラスに「X」「Y」というカラムを追加
②ジオメトリ演算でそれぞれにXとYの座標を付与(以下画像参照)
③「②」のデータをEXCELか何かで重複チェック

f:id:sanvarie:20190131104119p:plain

これで同じ結果が得られるかと思います。

ArcGISユーザーだけど、ArcPyは敷居が高くてなかなか手が出せないという方もいらっしゃるかと思います。ただ、とても便利なライブラリなので、ぜひ使ってみてください。何かお困りのことがあれば、それに関するコメントを残していただければできる限り回答いたします(私にわかる範囲でですが)

簡単にではありますが、今回は以上です。

PythonでGISデータを扱う際に便利なライブラリのまとめ

本日は PythonでGISデータを扱う際に便利なライブラリについてまとめてみようと思います。ほとんどが過去のエントリーで紹介しているものですが、これからGISデータを扱おうと思っている方やたくさんライブラリがある中でどういったものを選択していいのか迷っている方はぜひ参考にしてみてください。

メジャーなライブラリ

まずはこれを使っておけば間違いないというライブラリを紹介します。

GDAL

PythonのGIS系ライブラリと言えばまずGDALを思い浮かべるくらいメジャーなライブラリです。ベクターからラスター、また、データの読み込みから書き込みまで基本的なことは何でもできると考えて大丈夫です。過去のエントリーでもGDALに関していくつかの機能を紹介しています。

www.gis-py.com

www.gis-py.com

www.gis-py.com

www.gis-py.com

GeoPandas

上記で紹介したGDALも素晴らしいライブラリなのですが、最近はGeoPandasの方が人気かもしれません。こちらも基本的なことは何でもできますし、ますます機能が強化されていくと思います。もしかしたら将来はGeoPandasが超メジャーなGIS系のライブラリになるかもしれません。

www.gis-py.com

www.gis-py.com

www.gis-py.com

pyproj

座標変換用ライブラリです。座標変換をするならこれ以外に考えられません。

www.gis-py.com

geopy

ジオコーディング用ライブラリです。

www.gis-py.com

piexif

画像にジオタグを追加したり、画像に付与されているジオタグを抽出するためのライブラリです。

www.gis-py.com

Arcpy

ArcGISに限定されてしまうのですが、ArcGISを使っている方はぜひ使ってみてください。こちらは商用のエンジンに付随するライブラリですので、他のライブラリよりも機能は強力だと思います。以下のエントリー以外もArcpyに関しては本ブログで色々紹介していますので、ぜひ読んでみてください。

www.gis-py.com

www.gis-py.com

www.gis-py.com

玄人向けなライブラリ

少しマニアックなのですが、使いこなせると便利なライブラリです。

Fiona

データの読み書きに特化したシンプルなライブラリです。GeoPandasは一部機能をFionaに依存しています。Fionaを深く知ることによってGeoPandasについても深く知ることができるかと思います。

www.gis-py.com

Shapely

地理空間データを計算するためのライブラリです。こちらもFionaと同じでGeoPandasは一部機能をShapelyに依存しています。

www.gis-py.com

Folium

WEBマップを作成するためのライブラリです。WEBサイトにインタラクティブなマップを組み込みたいという場合はかなり使えると思います。

www.gis-py.com

Rasterio

本ブログでは未紹介なのですが、ラスタを扱うためのライブラリです。GDALでもラスタを扱うことはできますが、Rasterioはラスタに特化したライブラリなので、機能はこちらの方が充実しているかもしれません。

pypi.org


ざっと紹介しましたが、ほとんどのことはGeoPandas、pyproj、geopyあたりで事足りるような気がします。GIS系ライブラリについて勉強したいけど、何から手を付けていいかわからないという方はこの辺から始めるのがいいかもしれません。

また、ライブラリではないのですが、GISについて色々勉強してみたいという方は以下エントリーを参考にしてみてください。

www.gis-py.com

geopandasの使い方をマスターしよう ~Shapeファイルの読込・作成、GeoDataFrameの扱い方まで~

本日はgeopandasについて書いてみようと思います。geopandasについては以下のエントリーでも紹介しているのですが、今回はShapeファイルの読込・作成、GeoDataFrameの扱い方まで、一通りの使い方を紹介できればと思っています。

www.gis-py.com

www.gis-py.com

geopandasとは

geopandasはpandasの拡張で、地理データを含むデータをpandasのように表形式で扱うことができるGIS系のPythonライブラリです。

pandasについては以下エントリーでも紹介していますので、ぜひ読んでみてください。

www.gis-py.com

環境

Windows10
Python2.7.4
Jupyter notebook

サンプルデータ

以下エントリーで使用した日本地図を使用します。
www.gis-py.com

サンプルコード

Shapeファイル読込

import geopandas as gpd

fp = r"D:\gis-py\library\geopandas\data\japan_ver80.shp"
data = gpd.read_file(fp)

これでOKです。それでは”data”のデータタイプを確認してみようと思います。

type(data)


f:id:sanvarie:20181101153703p:plain

GeoDataFrameというデータタイプですね。pandasにDataFrame というデータタイプがあるのですが、それプラス地理データが格納されているものをGeoDataFrameと認識していただいて大丈夫だと思います。

.head()を使用することでGeoDataFrameの最初の5行を取得することができます。

data.head()


f:id:sanvarie:20181101153955p:plain

こんな風にやっても同じ結果になります。

data[0:5]

また、カラム指定してデータを取得することもできます。

data['geometry'].head()


f:id:sanvarie:20181101155250p:plain

.plot() を使用することで読み込んだShapeを視覚化することができます。

data.plot();

f:id:sanvarie:20181101154134p:plain

.crsを使用することで読み込んだShapeの座標系などを確認することができます。

data.crs


f:id:sanvarie:20181101154306p:plain

Shapeファイル作成

.to_file()でShapeファイルの作成ができます。

out = r"D:\gis-py\library\geopandas\data\japan_ver80_2.shp"
selection = data[0:50]
selection.to_file(out)

日本地図の最初の50レコードを別のShapeファイルとして作成しました。

fp_50 = r"D:\gis-py\library\geopandas\data\japan_ver80_2.shp"
data_50 = gpd.read_file(fp_50)
data_50.plot();


f:id:sanvarie:20181101154844p:plain

GeoDataFrameの扱い方

GeoDataFrameをループさせる場合はこのようにします。最初の5行を抽出して、それをループさせてみました。

selection = data[0:5]
for index, row in selection.iterrows():
    poly_area = row['geometry'].area
    print("Polygon area at index {0} is: {1:.3f}".format(index, poly_area))


f:id:sanvarie:20181101155712p:plain

新しいカラムを追加場合はこのようにします。

data['area'] = None


5行目までを確認してみました。
f:id:sanvarie:20181101155918p:plain

新しく作ったカラムにデータを格納します。geometryカラムのareaの値をareaカラムに格納します。

for index, row in data.iterrows():
    data.loc[index, 'area'] = row['geometry'].area


f:id:sanvarie:20181101160141p:plain

指定したカラムの最大値、最大値を取得することができます。

max_area = data['area'].max()
min_area = data['area'].mean()
print("Max area: %s\nMin area: %s" % (round(max_area, 2), round(min_area, 2)))


f:id:sanvarie:20181101160357p:plain

ジオメトリをGeoDataFrameに対して作ってみる

shapelyとfionaを使用します。geopandas扱う上で非常に重要なライブラリです。これらに関しては以下エントリーで紹介していますので、ぜひ読んでみてください。

www.gis-py.com

www.gis-py.com


空のGeoDataFrameを作る場合は、以下のようにします。

import pandas as pd
import geopandas as gpd
from shapely.geometry import Point, Polygon
import fiona

newdata = gpd.GeoDataFrame()

f:id:sanvarie:20181101160740p:plain


geometryカラムを作成してみます。

newdata['geometry'] = None


f:id:sanvarie:20181101161014p:plain


Shapelyを使用してポリゴンを作成します。

coordinates = [(24.950899, 60.169158), (24.953492, 60.169158), (24.953510, 60.170104), (24.950958, 60.169990)]
poly = Polygon(coordinates)
poly

f:id:sanvarie:20181101161155p:plain

作成したポリゴンをgeometryカラムに格納します。

newdata.loc[0, 'geometry'] = poly

f:id:sanvarie:20181101161257p:plain


Fionaを使用して測地系の設定を行います。

from fiona.crs import from_epsg
newdata.crs = from_epsg(6668)
newdata.crs

f:id:sanvarie:20181101161537p:plain

Shapeファイルを作成します。

outfp = r"D:\gis-py\library\geopandas\data\test.shp"
newdata.to_file(outfp)

f:id:sanvarie:20181101161842p:plain

いかがでしたでしょうか?これ一つで大抵のことはできてしまう超便利ライブラリです。GISとPythonを使用される方は必修といっても過言ではないかもしれません。今後はgeopandasの応用について書ければと思っています。皆さんもぜひ使ってみてください。