pyhton 离散点数据转栅格img
# -*-coding:utf-8 -*-
import arcpy
import os
from dbcommon.views import OperateDB
from arcpy import env
import shutil
import numpy as np
class DataToImg(object):
"""
:站点和格点数据转img图层公用算法
"""
def __init__(self):
#img文件临时存储路径
self.OUTPATH = 'F:/arcgis_temp/raster/'
#img文件存储路径
self.IMGPATH = 'F:/arcgis_temp/img/'
#环境变量设置 上下左右经纬度
self.x_min = 104.250
self.y_min = 35.200
self.x_max = 107.700
self.y_max = 39.400
def station_to_img(self,data,out_path,file_name):
"""
:站点数据转换为img图层
:param data :数据 二维数据[[经度，纬度，数据],………]
out_path :输出路径 例如：FGFP/
file_name:文件名称
"""
try:
env.overwriteOutput = True
#工作空间
env.workspace = self.OUTPATH + out_path
#环境变量设置 上下左右经纬度
env.extent = arcpy.Extent( self.x_min, self.y_min,
self.x_max, self.y_max
)
#环境变量设置 上下左右经纬度
img_path = self.discretePoint2raster(data, out_path, file_name)
#如果路径不存在，新建目录
if not os.path.exists(self.IMGPATH+out_path):
os.makedirs(self.IMGPATH+out_path)
file_name = self.IMGPATH+out_path+"/"+file_name+".img"
#将文件永久保存
shutil.copyfile(img_path,file_name)
#删除所有的临时文件
self.remove_all_file(self.OUTPATH+out_path)
return file_name
except:
raise
def grid_to_img(self,data,out_path,file_name,usecols=None):
"""
:站点数据转换为img图层
:param data :数据 二维数据[[经度，纬度，数据],………]
out_path :输出路径 例如：FGFP/
file_name:文件名称
usecols :文件参与绘图的列序号，判断传入数据是否为文件 例如（0,1,2）
"""
try:
env.overwriteOutput = True
#工作空间
env.workspace = self.OUTPATH + out_path
#环境变量设置 上下左右经纬度
env.extent = arcpy.Extent(
self.x_min, self.y_min, self.x_max, self.y_max)
#如果传入数据为文件路径则读取文件数据
if usecols:
context = np.genfromtxt(data, skip_header=1, usecols=usecols)
#如果传入数据为数据，则删除第一行（文件头）
else:
context = data[1:]
#将格点数据转化为img图层
img_path = self.txt2raster(context, out_path, file_name)
#如果路径不存在，新建目录
if not os.path.exists(self.IMGPATH+out_path):
os.makedirs(self.IMGPATH+out_path)
file_name = self.IMGPATH+out_path+"/"+file_name+".img"
#将文件永久保存
shutil.copyfile(img_path,file_name)
#删除所有的临时文件
self.remove_all_file(self.OUTPATH+out_path)
return file_name
except:
raise
def discretePoint2raster(self,data,out_path,file_name,mark=None):
"""
:离散点到img栅格图层转换 公共函数
:params data :待插值的离散点，列表类型
out_path :输出路径 例如：FGFP/
file_name:文件名称
mark :标识，判断是否重采样 1是 0否
:return :生成的img图层
"""
try:
#img文件输出路径
filepath = self.OUTPATH + out_path
#如果路径不存在，新建目录
if not os.path.exists(filepath):
os.makedirs(filepath)
#输出1X1km的img栅格数据
outRaster = filepath+"/"+file_name+".img"
#输出重采样后100m*100m的img栅格数据
outRasterRes=filepath+"/"+file_name+"_100.img"
newfc= file_name+".shp"
has_m = "DISABLED"
has_z = "ENABLED"
#指定空间参考坐标系
sr = arcpy.SpatialReference("WGS 1984")
#创建要素类(工作空间，要素类名称，要素类的几何类型，‘’，‘’，……
arcpy.CreateFeatureclass_management(filepath, newfc, "POINT", "", has_m, has_z, sr)
#添加字段，（要素类名称，字段名称，字段类型，字段位数，字段小数位数）
arcpy.AddField_management(newfc, 'z', 'double', '12', '8')
#向要素类种插入行
cur = arcpy.InsertCursor(newfc)
min_value= 9999999.99
max_value= -9999999.99
#遍历数据，插入到要素类种
for item in data:
x = item[0]
y = item[1]
z = float(item[2])
#插入及纬度
pnt = arcpy.Point(float(x), float(y))
if min_value>z:
min_value=z
if max_value<z:
max_value=z
#向要素类插入Z字段值（数据值）
row = cur.newRow()
row.Shape = pnt
row.z = z
cur.insertRow(row)
del row
#清除
del cur
zField = "z"
cellSize = "0.01"
power=2
# Set variables for search neighborhood
majSemiaxis = 300000
minSemiaxis = 300000
angle = 0
smoothFactor = 0.5
#平滑插值
searchNeighbourhood = arcpy.SearchNeighborhoodSmooth(
majSemiaxis, minSemiaxis,
angle, smoothFactor
)
#加权平均法插值（输入要素类，字段，‘’，输出要素类，像素大小，功率，平滑）
arcpy.IDW_ga(
newfc, zField, "", outRaster,
cellSize, power,
searchNeighbourhood
)
#判断是否进行重采样
if mark:
inRaster = arcpy.Raster(outRaster)
# 重新采样，提供分辨率，100m分辨率
arcpy.Resample_management(inRaster, outRasterRes, "0.001 0.001", "CUBIC")
return outRasterRes
else:
return outRaster
except:
raise
def remove_all_file(self,path):
"""
:删除文件夹下多有文件
"""
try:
for i in os.listdir(path):
path_file = os.path.join(path,i) #取文件绝对路径
os.remove(path_file)
except:
raise
def txt2raster(self,data,file_path,file_name,mark=None):
"""
:格点数据转换为img栅格图层 公共函数
:param: data :格点数据 <type 'numpy.ndarray'> [[经度，纬度，数据],……]
file_path:输出路径 例如：FGFP/
file_name:文件名称
mark :标识，判断是否重采样 1是 0否
:return out_raster_path :生成的img图层
"""
# 定义ASCII文件头信息
try:
no_data = '-9999.00'
#文件保存临时路径
out_path = self.OUTPATH + file_path + '/' + file_name + '.asc'
out_raster_path = self.OUTPATH + file_path + '/' + file_name + '.img'
if not os.path.exists(self.OUTPATH + file_path):
os.makedirs(self.OUTPATH + file_path)
context = data
# 计算经纬度递增步长
if context is not None and context[1][1] != context[0][1]:
cellsize = context[1][1] - context[0][1]
elif context is not None and context[1][0] != context[0][0]:
cellsize = context[1][1] - context[0][1]
cellsize = round(cellsize, 2)
ncols = int(round(round((self.x_max - self.x_min),2) / cellsize,2)) + 1
nrows = int(round(round((self.y_max - self.y_min),2) / cellsize,2)) + 1
# 初始化二维数组，默认值为-9999.00,存储ASCII数值
data = np.full((nrows, ncols), -9999.00, np.float64)
for item in context:
x = int(round(round((item[0] - self.x_min),2) / cellsize,2))
y = int(round(round((item[1] - self.y_min),2) / cellsize,2))
data[nrows - y][x] = item[2]
# 设置文件头样式
header =( "ncols\t{ncols}\nnrows\t{nrows}\nxllcorner\t{xllcorner}\nyllcorner\t{yllcorner}\n"
"cellsize\t{cellsize}\nNODATA_value\t{NODATA_value}")
header = header.format(ncols=ncols, nrows=nrows, xllcorner=self.x_min, yllcorner=self.y_min,
cellsize=cellsize, NODATA_value=no_data)
# 写入ASCII文件
np.savetxt(out_path, data, fmt='%.2f', delimiter='\t', header=header, comments='')
# 栅格数据的 ASCII文件转换为栅格数据集
# 指定坐标系为1984
sr = arcpy.SpatialReference(4326)
in_raster = arcpy.Raster(out_path)
arcpy.DefineProjection_management(in_raster, sr)
# 重新采样，提供分辨率，100m分辨率
arcpy.Resample_management(in_raster,out_raster_path, "0.001 0.001", "CUBIC")
return out_raster_path
except:
raise
if __name__=='__main__':
di = DataToImg()
"""u站点数据转img示例"""
sql = "select dss.longitude,dss.latitude ,ds.swwc from dbcommon_soilday ds ,dbcommon_stationsite dss\
where ds.station_id_c=dss.station_no\
and date_time =to_date('20170813','yyyymmdd')\
and ds.swwc !='None' "
op = OperateDB()
data = op.select(sql)
raster_path = di.station_to_img(data, 'test', 'test')
"""格点数据转img示例——输入参数为数据"""
path = r'F:/dataStorage/BUNP/CDCF/20170406.txt'
data = np.genfromtxt(path)
di.grid_to_img(data, 'test','test')
"""格点数据文件转img示例——输入参数为格点文件的文件路径"""
path = r'F:/dataStorage/BUNP/CDCF/20170406.txt'
di.grid_to_img(path, 'test','test')
转载自：https://blog.csdn.net/qq_16645423/article/details/80359207
