conda create -y -n tmp python=3.7
conda activate tmp
conda install -y -c conda-forge rioxarray
conda install -y -c conda-forge gdal
conda install -y -c conda-forge jupyter
import re
import glob
import pandas as pd
import numpy as np
import xarray as xr
import rioxarray
import os
def getcoords(word, xorg, yorg, dx, dy):
anum = lambda x : ord(x) - ord('a')
v = str(word)
if re.match('\d+', v):
iy, ix = int(v[0]), int(v[1])
else:
iy, ix = anum(v[0]), anum(v[1])
yc = yorg - iy*dy
xc = xorg + ix*dx
return xc, yc
# 4分割
def getcoords2(word, xorg, yorg, dx, dy):
v = int(word)
if v == 1:
iy, ix = 0, 0
elif v == 2:
iy, ix = 0, 1
elif v == 3:
iy, ix = 1, 0
elif v == 4:
iy, ix = 1, 1
else:
print('error')
yc = yorg - iy*dy
xc = xorg + ix*dx
return xc, yc
fs =glob.glob('*.txt')
for fp in fs:
f = os.path.basename(fp)
df = pd.read_csv(fp, header=None)
X = df[1].values
Y = df[2].values
Z = df[3].values
# float to int
Z = (100*Z).astype(np.int32)
w = f[:8]
nepsg = int(w[:2])
# 図郭の北西の座標を求める
# 地図情報レベル 50000
xc, yc = getcoords(w[2:4], xorg=-160000, yorg=300000, dx=40000, dy=30000)
# 地図情報レベル 5000
xc, yc = getcoords(w[4:6], xorg=xc, yorg=yc, dx=4000, dy=3000)
# 地図情報レベル 5000を4分割
xc, yc = getcoords2(w[6], xorg=xc, yorg=yc, dx=2000, dy=1500)
# # # さらに4分割
# xc, yc = getcoords2(w[7], xorg=xc, yorg=yc, dx=1000, dy=750)
delta = float(0.5)
dx=2000
dy=1500
ix = ((X - xc - 0.5*delta)/delta).astype(int)
iy = ((- Y + yc - 0.5*delta)/delta).astype(int)
zout= np.full((int(dx/delta),int(dy/delta)), int(-9999) )
# zout= np.full((int(dx/delta),int(dy/delta)), float(-9999) )
for ixp, iyp, zp in zip(ix, iy, Z) : zout[ixp,iyp] = zp
xcoord = np.arange(xc, xc+dx, delta) + 0.5*delta
ycoord = np.arange(yc, yc-dy, -delta) - 0.5*delta
epsg = str(6668 + nepsg)
ds = xr.Dataset({'z': (['y','x'], zout.T) }, coords={'x': xcoord, 'y': ycoord})
ds = ds.rio.write_crs('EPSG:' + epsg, inplace = True)
# ds.rio.reproject("epsg:****")
# export geotiff file
out = ds['z'].rio.to_raster( f[:-4] + '.tif')
del out
from osgeo import gdal
opt=gdal.BuildVRTOptions(VRTNodata=float(-9999), srcNodata=float(-9999))
my_vrt = gdal.BuildVRT('output.vrt', glob.glob( '*.tif'), options=opt)
del my_vrt