GeoJSONなどベクターデータをGeoTiffなどのラスターに変換する

from osgeo import gdal, ogr, osr
source = ogr.Open('src.geojson') # geojsonやShapefile
dst = 'dst.tiff' # 作成するファイル名
pixel_size = 10000 # ベクターの領域に合わせていい塩梅で決める
nodata = -9999 # テキトウに選択
value = 125 # intの範囲内でテキトウに選択

layer = source.GetLayer()
x_min, x_max, y_min, y_max = layer.GetExtent()
x_res = int((x_max - x_min) / pixel_size)
y_res = int((y_max - y_min) / pixel_size)

feature = layer.GetNextFeature()
geometry = feature.GetGeometryRef()
epsg = int(geometry.GetSpatialReference().GetAuthorityCode("PROJCS")) 
# データによってはうまくとれない可能性がある。だめなときはソースの投影座標系を手動で入れる。

target = gdal.GetDriverByName('GTiff').Create(dst, x_res, y_res, 1, gdal.GDT_Int16)
target.SetGeoTransform((x_min, pixel_size, 0, y_max, 0, -pixel_size))
band = target.GetRasterBand(1)
band.SetNoDataValue(nodata)
srs = osr.SpatialReference()
srs.ImportFromEPSG(epsg)
target.SetProjection(srs.ExportToWkt())
gdal.RasterizeLayer(target, [1], source, burn_values=[value])
target.FlushCache()
target = None

from PIL import Image
from PIL import ImageFile
import numpy as np
Image.MAX_IMAGE_PIXELS = 1000000000 # おまじない
ImageFile.LOAD_TRUNCATED_IMAGES = True # おまじない

dst = 'dst.tiff' # 作成したgeotiff
nodata = -9999 # 作成時に設定したnodata
dst_png = 'dst.png' # 作成するpng
rate = 1 # pngのpixel縮小率, 0 < rate <=1

im = Image.open(dst)
arr = np.array(im)
size = (int(im.size[0]*rate),int(im.size[1]*rate))

mask = np.copy(arr)
mask[arr!=nodata] = 255
mask[arr==nodata] = 0

outim = im.convert('RGBA')
outmask = Image.fromarray(mask).convert('L')

outim.putalpha(outmask)
outim = outim.resize(size)

outim.save(dst_png, 'png')

from osgeo import gdal, ogr
nodata = -9999 # tiffのNoDATAの値
layername = "tovector" # 出力名
ds = gdal.Open('source.tif')
drv = ogr.GetDriverByName("GeoJSON")
dst = drv.CreateDataSource( layername + ".geojson" )
layer = dst.CreateLayer(layername, srs = None )

drv2 = gdal.GetDriverByName( 'MEM' )
ds2 = drv2.CreateCopy('', ds)
srcband = ds2.GetRasterBand(1)
arr = srcband.ReadAsArray()
arr[arr[:, ]!=nodata] = 255
# 値ごとにポリゴンを分けてしまうのでNoDATA以外のピクセルはすべて同じ値にしてなるべくポリゴンが細かくならないようにする
srcband.WriteArray(arr)
srcband.SetNoDataValue(nodata)

gdal.Polygonize(srcband, srcband, layer, 0, [], callback=None )
dst.Destroy()
ds = None
ds2 = None

{
    "type": "FeatureCollection",
    "name": "tovector",
    "features": [{ "type": "Feature", "properties": { }, "geometry": { "type": "Polygon", "coordinates": [略]}]
}

def area(linearRing):
    ret = 0
    length = len(linearRing) - 1
    for i in range(length):
        ret = ret + linearRing[i][0] * linearRing[i + 1][1] - linearRing[i][1] * linearRing[i + 1][0]

    return abs(ret / 2);



def getCandidate(i, linearRing):
    if i == 0 or i == len(linearRing) - 1:
        return {
            'score': float('inf'),
            'index': i
        }

    return {
        'score': area([linearRing[i - 1], linearRing[i], linearRing[i + 1], linearRing[i - 1]]),
        'index': i
    }

def vw(linearRing, limit):
    candidates = [getCandidate(i, linearRing) for i in range(len(linearRing))]
    while (len(candidates) > 4):
        memory = float('inf')
        eliminate = -1
        for i in range(len(candidates)):
            if candidates[i]['score'] < memory:
                memory = candidates[i]['score']
                eliminate = i
         
        if len(candidates) > limit:
            if eliminate - 2 >= 0:
                p1 = candidates[eliminate - 2]['index']
                p2 = candidates[eliminate - 1]['index']
                p3 = candidates[eliminate + 1]['index']
                candidates[eliminate - 1]['score'] = area([
                    linearRing[p1], linearRing[p2], linearRing[p3], linearRing[p1]
                ])
            if eliminate + 2 <= len(candidates)-1:
                p1 = candidates[eliminate - 1]['index']
                p2 = candidates[eliminate + 1]['index']
                p3 = candidates[eliminate + 2]['index']
                candidates[eliminate + 1]['score'] = area([
                    linearRing[p1], linearRing[p2], linearRing[p3], linearRing[p1]
                ])
            candidates.pop(eliminate) 
        else:
            break;
    candidateIndexes = [v['index'] for v in candidates];
    filtered = []
    for i, v in enumerate(linearRing):
        if i in candidateIndexes:
            filtered.append(v)
    
    return filtered

import json
source = 'input.geojson'
dist = 'out.geojson'
limit = 20

with open(source, 'r') as f:
    json_load = json.load(f)
    print('before', len(json_load['features'][0]['geometry']['coordinates'][0]))
    decreased = vw(json_load['features'][0]['geometry']['coordinates'][0], limit)
    json_load['features'][0]['geometry']['coordinates'][0] = decreased
    print('after', len(json_load['features'][0]['geometry']['coordinates'][0]))

with open(dist, 'w') as f:
    json.dump(json_load, f)