translates csv file of "National Seismic Hazard Maps" provided by NIED into GeoTIFF files

j-shis_csv2geotiff.py

#!/usr/bin/env python
# j-shis_csv2geotiff.py
# translates csv files of "National Seismic Hazard Maps" provided by NIED into GeoTIFF files
# Copyright (c) 2013, Minoru Akagi
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program.  If not, see <http://www.gnu.org/licenses/>.

import sys
import csv
import gdal
import numpy

filename = "P-Y2012-MAP-AVR-TTL_MTTL.csv"
column_number = 3   # 3: T30_I55_PS
# columns: CODE, T30_I45_PS, T30_I50_PS, T30_I55_PS, T30_I60_PS, T30_P03_SI, T30_P03_BV, T30_P03_SV, T30_P06_SI, T30_P06_BV, T30_P06_SV, T50_P02_SI, T50_P02_BV, T50_P02_SV, T50_P05_SI, T50_P05_BV, T50_P05_SV, T50_P10_SI, T50_P10_BV, T50_P10_SV, T50_P39_SI, T50_P39_BV, T50_P39_SV

nodata_value = -9999.
xsize = ysize = 320
driver = gdal.GetDriverByName("GTiff")

def writeGeoTIFF(meshcode1, grid):
    dest_file = "%s.tif" % meshcode1

    uly = float(int(meshcode1[0:2]) + 1) * 2 / 3
    ulx = float(meshcode1[2:4]) + 100
    psize_x = 1. / xsize
    psize_y = 2. / 3 / ysize

    geotransform = [ulx, psize_x, 0, uly, 0, -psize_y]
    dst_ds = driver.Create(dest_file, xsize, ysize, 1, gdal.GDT_Float32, [])
    if dst_ds is None:
        sys.exit("Creation failed.")
    dst_ds.SetGeoTransform(geotransform)
    dst_ds.SetProjection('GEOGCS["Tokyo",DATUM["Tokyo",SPHEROID["Bessel 1841",6377397.155,299.1528128,AUTHORITY["EPSG","7004"]],AUTHORITY["EPSG","6301"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.01745329251994328,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4301"]]')

    rband = dst_ds.GetRasterBand(1)
    rband.WriteArray(narray, 0, 0)
    rband.SetNoDataValue(nodata_value)

reader = csv.reader(open(filename, "rb"), delimiter=",")
last_meshcode1 = 0
narray = None
for row in reader:
    if row[0][0] == "#":
        continue
    meshcode = row[0][0:10]
    meshcode1 = meshcode[0:4]
    p = float(row[column_number])
    if last_meshcode1 != meshcode1:
        print meshcode1
        if narray is not None:
            writeGeoTIFF(last_meshcode1, narray)
        narray = numpy.empty((ysize, xsize), numpy.float32)
        narray.fill(nodata_value)
        last_meshcode1 = meshcode1

    c8 = int(meshcode[8]) - 1
    c9 = int(meshcode[9]) - 1
    x = int(meshcode[5]) * 40 + int(meshcode[7]) * 4 + (c8 % 2) * 2 + (c9 % 2)
    y = ysize - 1 - (int(meshcode[4]) * 40 + int(meshcode[6]) * 4 + (c8 // 2) * 2 + (c9 // 2))
    narray[y][x] = p

writeGeoTIFF(meshcode1, narray)

J-SHISサイト(http://www.j-shis.bosai.go.jp/ )からダウンロードできる確率論的地震動予測地図の地図データ(CSVファイル形式)をGeoTIFFファイルに変換するスクリプトです。1次メッシュ毎にファイルが生成されます。