GIS奮闘記

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

スポンサーリンク

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.encode("shift-jis"))

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仕様

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

今回はデータの一部をダウンロードしてみましたが、全データをダウンロードして各データや年度ごとにフォルダ分けしてダウンロードファイルを格納する、みたいなプログラムを作ったらすごく面白そうですね。

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

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の応用について書ければと思っています。皆さんもぜひ使ってみてください。

PythonでGDALを使い倒そう!

さて、本日はGDALについて書いてみようと思います。GDLについては以下のように以前に何個か記事を書いているのですが、非常に奥が深いライブラリですので紹介していない機能がまだまだたくさん残っています。今回はラスターの扱いについて紹介します。

www.gis-py.com

www.gis-py.com

www.gis-py.com

GDALとは

公式サイトによると、GDALはフリーのGISライブラリでラスターとベクターを扱うことができます。よくGDALはラスター用のライブラリと書かれているサイトを見るのですが、ベクターも扱うことができます。GDAL: GDAL - Geospatial Data Abstraction Library

サンプルデータ

今回はラスターの読み込みとラスターの作成を行います。以下二つのデータを用意してください。

ラスターの読み込み

以下サイトで航空写真がダウンロードできますので、今回はそのうちの一つを使ってみようと思います。
空中写真 オープンデータ

航空写真とは?という方は以下エントリーを参考にしてみてください。
www.gis-py.com

ラスターの作成

ShapeファイルをGEOTIFFに変換してみます。以下エントリーで使用した日本地図をGEOTIFFに変換します。

www.gis-py.com

サンプルコード

ラスターの読み込み

from osgeo import gdal
raster = gdal.Open( r'D:\gis-py\library\gdal\data\19640721c-wide-niigata-eq-0-1.tif')

print raster.GetProjection() #投影法
print raster.RasterXSize      #Xサイズ
print raster.RasterYSize      #Yサイズ
print raster.RasterCount    #バンド数
print raster.GetMetadata() #メタデータ
print raster.GetDriver().LongName

geotransform= raster.GetGeoTransform()
print geotransform[0] #左上X
print geotransform[3] #左上Y
print geotransform[1] #セルサイズX
print geotransform[5] #セルサイズY

このようにラスターのプロパティを取得することができます。
f:id:sanvarie:20181029161812p:plain

ラスターの作成
from osgeo import gdal, ogr

# ピクセルサイズを設定
pixel_size = 0.01
NoData_value = -9999

# Shapeファイルを指定
vector_fn = ur'D:\gis-py\library\gdal\data\japan_ver80.shp'

# 作成するTIFFを指定
raster_fn = r'D:\gis-py\library\gdal\outcome\japan_geo.tif'

# Shapeをオープン
source_ds = ogr.Open(vector_fn)
source_layer = source_ds.GetLayer()
x_min, x_max, y_min, y_max = source_layer.GetExtent()

# TIFF作成のための設定
x_res = int((x_max - x_min) / pixel_size)
y_res = int((y_max - y_min) / pixel_size)
target_ds = gdal.GetDriverByName('GTiff').Create(raster_fn, x_res, y_res, 1, gdal.GDT_Int16)
target_ds.SetGeoTransform((x_min, pixel_size, 0, y_max, 0, -pixel_size))
band = target_ds.GetRasterBand(1)
band.SetNoDataValue(NoData_value)

# ラスタ化
gdal.RasterizeLayer(target_ds, [1], source_layer, burn_values=[0])


ばっちり日本地図のShapeをGEOTIFFに変換することができました。
f:id:sanvarie:20181031165923p:plain


ラスターの読み込みで使用したソースで、作成したGEOTIFFのプロパティを確認しました。座標などが設定されていることがわかりました。
f:id:sanvarie:20181031164809p:plain

今回は以上ですが、まだまだたくさんの機能がありますので、それはまた別の機会に紹介できればと思います。本当に便利なライブラリなので、ぜひ使ってみてください。