技术标签： 原型模式
# -*- coding: utf-8 -*-
import os
import numpy as np
from osgeo import gdal, gdalnumeric, ogr, osr, gdal_array
gdal.UseExceptions()
def world2Pixel(geoMatrix, x, y):
"""
Uses a gdal geomatrix (gdal.GetGeoTransform()) to calculate
the pixel location of a geospatial coordinate
"""
ulX = geoMatrix[0]
ulY = geoMatrix[3]
xDist = geoMatrix[1]
yDist = geoMatrix[5]
rtnX = geoMatrix[2]
rtnY = geoMatrix[4]
pixel = int((x - ulX) / xDist)
line = int((ulY - y) / xDist)
return (pixel, line)
#
# EDIT: this is basically an overloaded
# version of the gdal_array.OpenArray passing in xoff, yoff explicitly
# so we can pass these params off to CopyDatasetInfo
#
def OpenArray( array, prototype_ds = None, xoff=0, yoff=0 ):
# ds = gdal.Open( gdalnumeric.GetArrayFilename(array))
ds = gdal_array.OpenArray(array)
if ds is not None and prototype_ds is not None:
if type(prototype_ds).__name__ == 'str':
prototype_ds = gdal.Open( prototype_ds )
if prototype_ds is not None:
gdalnumeric.CopyDatasetInfo( prototype_ds, ds, xoff=xoff, yoff=yoff )
return ds
def write_img(filename,im_proj,im_geotrans,im_data):
if 'int8' in im_data.dtype.name:
datatype = gdal.GDT_Byte
elif 'int16' in im_data.dtype.name:
datatype = gdal.GDT_UInt16
else:
datatype = gdal.GDT_Float32
if len(im_data.shape) == 3:
im_bands, im_height, im_width = im_data.shape
else:
im_bands, (im_height, im_width) = 1,im_data.shape
driver = gdal.GetDriverByName("GTiff")
dataset = driver.Create(filename, im_width, im_height, im_bands, datatype)
dataset.SetGeoTransform(im_geotrans)
dataset.SetProjection(im_proj)
if im_bands == 1:
dataset.GetRasterBand(1).WriteArray(im_data)
else:
for i in range(im_bands):
dataset.GetRasterBand(i+1).WriteArray(im_data[i])
del dataset
def shpClipRaster(shapefile_path, raster_path, save_path):
# Load the source data as a gdalnumeric array
# srcArray = gdalnumeric.LoadFile(raster_path)
# Also load as a gdal image to get geotransform
# (world file) info
srcImage = gdal.Open(raster_path)
geoTrans = srcImage.GetGeoTransform()
geoProj = srcImage.GetProjection()
# Create an OGR layer from a boundary shapefile
shapef = ogr.Open(shapefile_path)
lyr = shapef.GetLayer( os.path.split( os.path.splitext( shapefile_path )[0] )[1] )
poly = lyr.GetNextFeature()
# Convert the layer extent to image pixel coordinates
minX, maxX, minY, maxY = lyr.GetExtent()
ulX, ulY = world2Pixel(geoTrans, minX, maxY)
lrX, lrY = world2Pixel(geoTrans, maxX, minY)
# Calculate the pixel size of the new image
pxWidth = int(lrX - ulX)
pxHeight = int(lrY - ulY)
# clip = srcArray[:, ulY:lrY, ulX:lrX]
clip = srcImage.ReadAsArray(ulX,ulY,pxWidth,pxHeight) #***只读要的那块***
#
# EDIT: create pixel offset to pass to new image Projection info
#
xoffset = ulX
yoffset = ulY
print ("Xoffset, Yoffset = ( %f, %f )" % ( xoffset, yoffset ))
# Create a new geomatrix for the image
geoTrans = list(geoTrans)
geoTrans[0] = minX
geoTrans[3] = maxY
write_img(save_path, geoProj, geoTrans, clip)
gdal.ErrorReset()
if __name__ == "__main__":
shp = "dataset/E22_Bound.shp"
img = "dataset/CGdomYRJ-114(CK0-17)_E_22.tif"
out = "dataset/E22.tif"
shpClipRaster(shp,img,out)
print(img)
import os
from osgeo import gdal, ogr
def ShapeClip(
baseFilePath,
maskFilePath,
saveFolderPath):
"""
矢量裁剪
:param baseFilePath: 要裁剪的矢量文件
:param maskFilePath: 掩膜矢量文件
:param saveFolderPath: 裁剪后的矢量文件保存目录
:return:
"""
ogr.RegisterAll()
gdal.SetConfigOption("GDAL_FILENAME_IS_UTF8", "YES")
# 载入要裁剪的矢量文件
baseData = ogr.Open(baseFilePath)
print(os.path.split( os.path.splitext( baseFilePath )[0] )[1])
baseLayer = baseData.GetLayer( os.path.split( os.path.splitext( baseFilePath )[0] )[1] )
spatial = baseLayer.GetSpatialRef()
geomType = baseLayer.GetGeomType()
baseLayerName = baseLayer.GetName()
# 载入掩膜矢量文件
maskData = ogr.Open(maskFilePath)
maskLayer = maskData.GetLayer()
maskLayerName = maskLayer.GetName()
# 生成裁剪后的矢量文件
outLayerName = maskLayerName + "_Clip_" + baseLayerName
outFilePath = saveFolderPath
gdal.SetConfigOption("SHAPE_ENCODING", "GBK")
driver = ogr.GetDriverByName("ESRI Shapefile")
outData = driver.CreateDataSource(outFilePath)
outLayer = outData.CreateLayer(outLayerName, spatial, geomType)
baseLayer.Clip(maskLayer, outLayer)
outData.Release()
baseData.Release()
maskData.Release()
return outFilePath
if __name__ == "__main__":
baseFilePath = 'dataset/veg_E_22.shp'
maskFilePath = 'dataset/E22_Bound.shp'
saveFolderPath = 'dataset/E22.shp'
outFilePath=ShapeClip(baseFilePath,maskFilePath,saveFolderPath)
print(outFilePath)
from osgeo import gdal, ogr, gdalconst
def shp2Raster(shp,templatePic,output,nodata):
"""
shp:字符串，一个矢量，从0开始计数，整数
templatePic:字符串，模板栅格，一个tif，地理变换信息从这里读，栅格大小与该栅格一致
output:字符串，输出栅格，一个tif
field:字符串，栅格值的字段
nodata:整型或浮点型，矢量空白区转换后的值
"""
ndsm = templatePic
data = gdal.Open(ndsm, gdalconst.GA_ReadOnly)
geo_transform = data.GetGeoTransform()
proj=data.GetProjection()
#source_layer = data.GetLayer()
x_min = geo_transform[0]
y_max = geo_transform[3]
x_max = x_min + geo_transform[1] * data.RasterXSize
y_min = y_max + geo_transform[5] * data.RasterYSize
x_res = data.RasterXSize
y_res = data.RasterYSize
mb_v = ogr.Open(shp)
mb_l = mb_v.GetLayer()
pixel_width = geo_transform[1]
#输出影像为24位整型
target_ds = gdal.GetDriverByName('GTiff').Create(output, x_res, y_res, 1, gdal.GPI_RGB)
target_ds.SetGeoTransform(geo_transform)
target_ds.SetProjection(proj)
band = target_ds.GetRasterBand(1)
NoData_value = nodata
band.SetNoDataValue(NoData_value)
band.FlushCache()
gdal.RasterizeLayer(target_ds, [1], mb_l, options=['ALL_TOUCHED=TRUE'])
target_ds = None
if __name__ == "__main__":
shp = "dataset/E22.shp"
templatePic= "dataset/E22.tif"
output = "dataset/E22_mask.tif"
nodata=0
shp2Raster(shp,templatePic,output,nodata)
import os
import os
import numpy as np
from osgeo import gdal, gdalnumeric, ogr, osr, gdal_array
gdal.UseExceptions()
def world2Pixel(geoMatrix, x, y):
"""
Uses a gdal geomatrix (gdal.GetGeoTransform()) to calculate
the pixel location of a geospatial coordinate
"""
ulX = geoMatrix[0]
ulY = geoMatrix[3]
xDist = geoMatrix[1]
yDist = geoMatrix[5]
rtnX = geoMatrix[2]
rtnY = geoMatrix[4]
pixel = int((x - ulX) / xDist)
line = int((ulY - y) / xDist)
return (pixel, line)
#
# EDIT: this is basically an overloaded
# version of the gdal_array.OpenArray passing in xoff, yoff explicitly
# so we can pass these params off to CopyDatasetInfo
#
def OpenArray( array, prototype_ds = None, xoff=0, yoff=0 ):
# ds = gdal.Open( gdalnumeric.GetArrayFilename(array))
ds = gdal_array.OpenArray(array)
if ds is not None and prototype_ds is not None:
if type(prototype_ds).__name__ == 'str':
prototype_ds = gdal.Open( prototype_ds )
if prototype_ds is not None:
gdalnumeric.CopyDatasetInfo( prototype_ds, ds, xoff=xoff, yoff=yoff )
return ds
def write_img(filename,im_proj,im_geotrans,im_data):
if 'int8' in im_data.dtype.name:
datatype = gdal.GDT_Byte
elif 'int16' in im_data.dtype.name:
datatype = gdal.GDT_UInt16
else:
datatype = gdal.GDT_Float32
if len(im_data.shape) == 3:
im_bands, im_height, im_width = im_data.shape
else:
im_bands, (im_height, im_width) = 1,im_data.shape
driver = gdal.GetDriverByName("GTiff")
dataset = driver.Create(filename, im_width, im_height, im_bands, datatype)
dataset.SetGeoTransform(im_geotrans)
dataset.SetProjection(im_proj)
if im_bands == 1:
dataset.GetRasterBand(1).WriteArray(im_data)
else:
for i in range(im_bands):
dataset.GetRasterBand(i+1).WriteArray(im_data[i])
del dataset
pre_path='dataset/pre/'
labellist = filter(lambda x: x.find('label')!=-1, os.listdir(pre_path))
list1 = list(map(lambda x: x[:], labellist))
label_name=pre_path + list1[0]
boundarylist = filter(lambda x: x.find('shp')!=-1, os.listdir(pre_path+'boundary/'))
list2 = list(map(lambda x: x[:], boundarylist))
imagelist = filter(lambda x: x.find('tif')!=-1, os.listdir(pre_path))
list3 = list(map(lambda x: x[:], imagelist))
img_path=pre_path + list3[0]
"""
矢量裁剪
:param label_name: 要裁剪的矢量文件
:param boundary_name: 掩膜矢量文件
img_path: 影像
:param saveFolderPath: 裁剪后的矢量文件保存目录
:return:
"""
ogr.RegisterAll()
gdal.SetConfigOption("GDAL_FILENAME_IS_UTF8", "YES")
# 载入要裁剪的矢量文件
labelData = ogr.Open(label_name)
labelLayer = labelData.GetLayer( os.path.split( os.path.splitext( label_name )[0] )[1] )
spatial = labelLayer.GetSpatialRef()
geomType = labelLayer.GetGeomType()
# 载入掩膜矢量文件
def new_func(outLayerName):
return outLayerName
for i in list2:
boundary_name=pre_path+'boundary/'+ i
maskData = ogr.Open(boundary_name)
maskLayer = maskData.GetLayer()
#裁剪shp
# 生成裁剪后的矢量文件
save_shp_dir='./dataset/pre/shp/'
if not os.path.exists(save_shp_dir):
os.mkdir(save_shp_dir)
outLayerName = (save_shp_dir+i)
gdal.SetConfigOption("SHAPE_ENCODING", "GBK")
driver = ogr.GetDriverByName("ESRI Shapefile")
outData = driver.CreateDataSource(outLayerName)
outLayer = outData.CreateLayer(new_func(outLayerName), spatial, geomType)
labelLayer.Clip(maskLayer, outLayer)
outData.Release()
maskData.Release()
#裁剪tif
shp = "dataset/E22_Bound.shp"
img = "dataset/CGdomYRJ-114(CK0-17)_E_22.tif"
out = "dataset/E22.tif"
# Load the source data as a gdalnumeric array
# srcArray = gdalnumeric.LoadFile(raster_path)
# Also load as a gdal image to get geotransform
# (world file) info
srcImage = gdal.Open(raster_path)
geoTrans = srcImage.GetGeoTransform()
geoProj = srcImage.GetProjection()
# Create an OGR layer from a boundary shapefile
shapef = ogr.Open(shapefile_path)
lyr = shapef.GetLayer( os.path.split( os.path.splitext( shapefile_path )[0] )[1] )
poly = lyr.GetNextFeature()
# Convert the layer extent to image pixel coordinates
minX, maxX, minY, maxY = lyr.GetExtent()
ulX, ulY = world2Pixel(geoTrans, minX, maxY)
lrX, lrY = world2Pixel(geoTrans, maxX, minY)
# Calculate the pixel size of the new image
pxWidth = int(lrX - ulX)
pxHeight = int(lrY - ulY)
# clip = srcArray[:, ulY:lrY, ulX:lrX]
clip = srcImage.ReadAsArray(ulX,ulY,pxWidth,pxHeight) #***只读要的那块***
#
# EDIT: create pixel offset to pass to new image Projection info
#
xoffset = ulX
yoffset = ulY
print ("Xoffset, Yoffset = ( %f, %f )" % ( xoffset, yoffset ))
# Create a new geomatrix for the image
geoTrans = list(geoTrans)
geoTrans[0] = minX
geoTrans[3] = maxY
write_img(save_path, geoProj, geoTrans, clip)
gdal.ErrorReset()
labelData.Release()#影像裁剪shp,转为栅格，为该影像标签
#输入：存放影像文件夹dataset/sat_train，存放标签矢量文件夹dataset/mask_shp
#输出：标签（栅格），存放在dataset/mask_train
from osgeo import gdal, ogr, osr, gdalconst
import fnmatch
import os
def ShapeClip(
baseFilePath,
maskFilePath,
saveFolderPath):
"""
矢量裁剪
:param baseFilePath: 要裁剪的矢量文件
:param maskFilePath: 掩膜矢量文件
:param saveFolderPath: 裁剪后的矢量文件保存目录
:return:
"""
ogr.RegisterAll()
gdal.SetConfigOption("GDAL_FILENAME_IS_UTF8", "YES")
# 载入要裁剪的矢量文件
baseData = ogr.Open(baseFilePath)
baseLayer = baseData.GetLayer( os.path.split( os.path.splitext( baseFilePath )[0] )[1] )
spatial = baseLayer.GetSpatialRef()
geomType = baseLayer.GetGeomType()
baseLayerName = baseLayer.GetName()
# 载入掩膜矢量文件
maskData = ogr.Open(maskFilePath)
maskLayer = maskData.GetLayer()
maskLayerName = maskLayer.GetName()
# 生成裁剪后的矢量文件
outLayerName = maskLayerName + "_Clip_" + baseLayerName
outFilePath = saveFolderPath
gdal.SetConfigOption("SHAPE_ENCODING", "GBK")
driver = ogr.GetDriverByName("ESRI Shapefile")
outData = driver.CreateDataSource(outFilePath)
outLayer = outData.CreateLayer(outLayerName, spatial, geomType)
baseLayer.Clip(maskLayer, outLayer)
outData.Release()
baseData.Release()
maskData.Release()
return outFilePath
def shp2Raster(shp,templatePic,output,nodata):
"""
shp:字符串，一个矢量，从0开始计数，整数
templatePic:字符串，模板栅格，一个tif，地理变换信息从这里读，栅格大小与该栅格一致
output:字符串，输出栅格，一个tif
field:字符串，栅格值的字段
nodata:整型或浮点型，矢量空白区转换后的值
"""
ndsm = templatePic
data = gdal.Open(ndsm, gdalconst.GA_ReadOnly)
geo_transform = data.GetGeoTransform()
proj=data.GetProjection()
#source_layer = data.GetLayer()
x_min = geo_transform[0]
y_max = geo_transform[3]
x_max = x_min + geo_transform[1] * data.RasterXSize
y_min = y_max + geo_transform[5] * data.RasterYSize
x_res = data.RasterXSize
y_res = data.RasterYSize
mb_v = ogr.Open(shp)
mb_l = mb_v.GetLayer()
pixel_width = geo_transform[1]
#输出影像为24位整型
target_ds = gdal.GetDriverByName('GTiff').Create(output, x_res, y_res, 1, gdal.GPI_RGB)
target_ds.SetGeoTransform(geo_transform)
target_ds.SetProjection(proj)
band = target_ds.GetRasterBand(1)
NoData_value = nodata
band.SetNoDataValue(NoData_value)
band.FlushCache()
gdal.RasterizeLayer(target_ds, [1], mb_l, options=['ALL_TOUCHED=TRUE'])
target_ds = None
print("开始制作标签")
ogr.RegisterAll()
img_path="dataset/sat_train/" #影像所在的文件夹
mask_shp_path="dataset/mask_shp/" #原始标签shp位置
shape_path="dataset/mask_boundary_shp/" #shape输出位置
mask_clip_path='dataset/mask_clip_train/'#裁剪后shp
mask_train_path='dataset/mask_train/'#最终输出标签文件夹
if not os.path.exists(shape_path):
os.mkdir(shape_path)
if not os.path.exists(mask_clip_path):
os.mkdir(mask_clip_path)
imagelist = filter(lambda x: x.find('shp')!=-1, os.listdir(mask_shp_path))
list = list(map(lambda x: x[:], imagelist))
mask_shp_name=mask_shp_path + list[0]
img_list = fnmatch.filter(os.listdir(img_path), '*.tif')
for img in img_list:
p_img=img_path+img
outfilename = shape_path+img[:-4]+".shp"
dataset = gdal.Open(p_img)
oDriver = ogr.GetDriverByName('ESRI Shapefile')
oDS = oDriver.CreateDataSource(outfilename)
srs = osr.SpatialReference(wkt=dataset.GetProjection())
geocd = dataset.GetGeoTransform()
oLayer = oDS.CreateLayer("polygon", srs, ogr.wkbPolygon)
oDefn = oLayer.GetLayerDefn()
row = dataset.RasterXSize
line = dataset.RasterYSize
geoxmin = geocd[0]
geoymin = geocd[3]
geoxmax = geocd[0] + (row) * geocd[1] + (line) * geocd[2]
geoymax = geocd[3] + (row) * geocd[4] + (line) * geocd[5]
ring = ogr.Geometry(ogr.wkbLinearRing)
ring.AddPoint(geoxmin, geoymin)
ring.AddPoint(geoxmax, geoymin)
ring.AddPoint(geoxmax, geoymax)
ring.AddPoint(geoxmin, geoymax)
ring.CloseRings()
poly = ogr.Geometry(ogr.wkbPolygon)
poly.AddGeometry(ring)
outfeat = ogr.Feature(oDefn)
outfeat.SetGeometry(poly)
oLayer.CreateFeature(outfeat)
outfeat = None
oDS.Destroy()
mask_train_name = mask_clip_path+img[:-4]+".shp"
#裁剪
outFilePath=ShapeClip(mask_shp_name,outfilename, mask_train_name)
#矢量转栅格
output = mask_train_path + img
nodata=0
shp2Raster(mask_train_name,p_img,output,nodata)
print(output)
print('标签制作完成')做mask
'根据多个给定范围shp，对画好的标签进行裁剪并转栅格，做为标签样本，对影像进行裁剪，作为影像样本'
'输入：'
'输出'
import os
import os
import numpy as np
from osgeo import gdal, gdalnumeric, ogr, osr, gdal_array
gdal.UseExceptions()
def world2Pixel(geoMatrix, x, y):
"""
Uses a gdal geomatrix (gdal.GetGeoTransform()) to calculate
the pixel location of a geospatial coordinate
"""
ulX = geoMatrix[0]
ulY = geoMatrix[3]
xDist = geoMatrix[1]
yDist = geoMatrix[5]
rtnX = geoMatrix[2]
rtnY = geoMatrix[4]
pixel = int((x - ulX) / xDist)
line = int((ulY - y) / xDist)
return (pixel, line)
#
# EDIT: this is basically an overloaded
# version of the gdal_array.OpenArray passing in xoff, yoff explicitly
# so we can pass these params off to CopyDatasetInfo
#
def OpenArray( array, prototype_ds = None, xoff=0, yoff=0 ):
# ds = gdal.Open( gdalnumeric.GetArrayFilename(array))
ds = gdal_array.OpenArray(array)
if ds is not None and prototype_ds is not None:
if type(prototype_ds).__name__ == 'str':
prototype_ds = gdal.Open( prototype_ds )
if prototype_ds is not None:
gdalnumeric.CopyDatasetInfo( prototype_ds, ds, xoff=xoff, yoff=yoff )
return ds
def write_img(filename,im_proj,im_geotrans,im_data):
if 'int8' in im_data.dtype.name:
datatype = gdal.GDT_Byte
elif 'int16' in im_data.dtype.name:
datatype = gdal.GDT_UInt16
else:
datatype = gdal.GDT_Float32
if len(im_data.shape) == 3:
im_bands, im_height, im_width = im_data.shape
else:
im_bands, (im_height, im_width) = 1,im_data.shape
driver = gdal.GetDriverByName("GTiff")
dataset = driver.Create(filename, im_width, im_height, im_bands, datatype)
dataset.SetGeoTransform(im_geotrans)
dataset.SetProjection(im_proj)
if im_bands == 1:
dataset.GetRasterBand(1).WriteArray(im_data)
else:
for i in range(im_bands):
dataset.GetRasterBand(i+1).WriteArray(im_data[i])
del dataset
def shp2Raster(shp,templatePic,output,nodata):
"""
shp:字符串，一个矢量，从0开始计数，整数
templatePic:字符串，模板栅格，一个tif，地理变换信息从这里读，栅格大小与该栅格一致
output:字符串，输出栅格，一个tif
field:字符串，栅格值的字段
nodata:整型或浮点型，矢量空白区转换后的值
"""
ndsm = templatePic
data = gdal.Open(ndsm, gdalconst.GA_ReadOnly)
geo_transform = data.GetGeoTransform()
proj=data.GetProjection()
#source_layer = data.GetLayer()
x_min = geo_transform[0]
y_max = geo_transform[3]
x_max = x_min + geo_transform[1] * data.RasterXSize
y_min = y_max + geo_transform[5] * data.RasterYSize
x_res = data.RasterXSize
y_res = data.RasterYSize
mb_v = ogr.Open(shp)
mb_l = mb_v.GetLayer()
pixel_width = geo_transform[1]
#输出影像为24位整型
target_ds = gdal.GetDriverByName('GTiff').Create(output, x_res, y_res, 1, gdal.GPI_RGB)
target_ds.SetGeoTransform(geo_transform)
target_ds.SetProjection(proj)
band = target_ds.GetRasterBand(1)
NoData_value = nodata
band.SetNoDataValue(NoData_value)
band.FlushCache()
gdal.RasterizeLayer(target_ds, [1], mb_l, options=['ALL_TOUCHED=TRUE'])
target_ds = None
pre_path='dataset/pre/'
mask_train_path='dataset/mask_train/'#最终输出标签文件夹
labellist = filter(lambda x: x.find('label')!=-1, os.listdir(pre_path))
list1 = list(map(lambda x: x[:], labellist))
label_name=pre_path + list1[0]
boundarylist = filter(lambda x: x.find('.shp')!=-1, os.listdir(pre_path+'boundary/'))
list2 = list(map(lambda x: x[:], boundarylist))
imagelist = filter(lambda x: x.find('tif')!=-1, os.listdir(pre_path+'img'))
list3 = list(map(lambda x: x[:], imagelist))
img_path=pre_path + list3[0]
"""
矢量裁剪
:param label_name: 要裁剪的矢量文件
:param boundary_name: 掩膜矢量文件
img_path: 影像
:param saveFolderPath: 裁剪后的矢量文件保存目录
:return:
"""
print('开始用矢量范围裁剪影像')
ogr.RegisterAll()
gdal.SetConfigOption("GDAL_FILENAME_IS_UTF8", "YES")
# 载入要裁剪的矢量文件
labelData = ogr.Open(label_name)
labelLayer = labelData.GetLayer( os.path.split( os.path.splitext( label_name )[0] )[1] )
spatial = labelLayer.GetSpatialRef()
geomType = labelLayer.GetGeomType()
# 载入掩膜矢量文件
def new_func(outLayerName):
return outLayerName
for i in list2:
boundary_name=pre_path+'boundary/'+ i
maskData = ogr.Open(boundary_name)
maskLayer = maskData.GetLayer()
#裁剪shp
# 生成裁剪后的矢量文件
save_shp_dir='./dataset/pre/shp/'
if not os.path.exists(save_shp_dir):
os.mkdir(save_shp_dir)
outLayerName = (save_shp_dir+i)
gdal.SetConfigOption("SHAPE_ENCODING", "GBK")
driver = ogr.GetDriverByName("ESRI Shapefile")
outData = driver.CreateDataSource(outLayerName)
outLayer = outData.CreateLayer(new_func(outLayerName), spatial, geomType)
labelLayer.Clip(maskLayer, outLayer)
lyr = maskData.GetLayer( os.path.split( os.path.splitext( boundary_name )[0] )[1] )
shpminX, shpmaxX, shpminY, shpmaxY = lyr.GetExtent()
#裁剪tif
flag=0
for j in list3:
raster_path = pre_path+'img/'+j
srcImage = gdal.Open(raster_path)
geocd = srcImage.GetGeoTransform()
geoProj = srcImage.GetProjection()
row = srcImage.RasterXSize
line = srcImage.RasterYSize
tifxmin = geocd[0]
tifymin = geocd[3]
tifxmax = geocd[0] + (row) * geocd[1] + (line) * geocd[2]
tifymax = geocd[3] + (row) * geocd[4] + (line) * geocd[5]
if shpminX>=tifxmin and shpmaxX<=tifxmax and shpminY<=tifymin and shpmaxY>=tifymax:
ulX, ulY = world2Pixel(geocd, shpminX, shpmaxY)
lrX, lrY = world2Pixel(geocd, shpmaxX, shpminY)
# Calculate the pixel size of the new image
pxWidth = int(lrX - ulX)
pxHeight = int(lrY - ulY)
clip = srcImage.ReadAsArray(ulX,ulY,pxWidth,pxHeight) #***只读要的那块***
xoffset = ulX
yoffset = ulY
geoTrans = list(geoTrans)
geoTrans[0] = shpminX
geoTrans[3] = shpmaxY
save_path='dataset/sat_train/'+i[:-4]+'.tif'
write_img(save_path, geoProj, geoTrans, clip)
gdal.ErrorReset()
outData.Release()
maskData.Release()
flag=1
output = mask_train_path + i[:-4] +'.tif'
nodata=0
shp2Raster(outLayerName,save_path,output,nodata)
if flag==0:
print(raster_path+"没有制作")
else:
print(raster_path)
labelData.Release()
做了一半的
import os
import os
import numpy as np
from osgeo import gdal, gdalnumeric, ogr, osr, gdal_array
gdal.UseExceptions()
def world2Pixel(geoMatrix, x, y):
"""
Uses a gdal geomatrix (gdal.GetGeoTransform()) to calculate
the pixel location of a geospatial coordinate
"""
ulX = geoMatrix[0]
ulY = geoMatrix[3]
xDist = geoMatrix[1]
yDist = geoMatrix[5]
rtnX = geoMatrix[2]
rtnY = geoMatrix[4]
pixel = int((x - ulX) / xDist)
line = int((ulY - y) / xDist)
return (pixel, line)
#
# EDIT: this is basically an overloaded
# version of the gdal_array.OpenArray passing in xoff, yoff explicitly
# so we can pass these params off to CopyDatasetInfo
#
def OpenArray( array, prototype_ds = None, xoff=0, yoff=0 ):
# ds = gdal.Open( gdalnumeric.GetArrayFilename(array))
ds = gdal_array.OpenArray(array)
if ds is not None and prototype_ds is not None:
if type(prototype_ds).__name__ == 'str':
prototype_ds = gdal.Open( prototype_ds )
if prototype_ds is not None:
gdalnumeric.CopyDatasetInfo( prototype_ds, ds, xoff=xoff, yoff=yoff )
return ds
def write_img(filename,im_proj,im_geotrans,im_data):
if 'int8' in im_data.dtype.name:
datatype = gdal.GDT_Byte
elif 'int16' in im_data.dtype.name:
datatype = gdal.GDT_UInt16
else:
datatype = gdal.GDT_Float32
if len(im_data.shape) == 3:
im_bands, im_height, im_width = im_data.shape
else:
im_bands, (im_height, im_width) = 1,im_data.shape
driver = gdal.GetDriverByName("GTiff")
dataset = driver.Create(filename, im_width, im_height, im_bands, datatype)
dataset.SetGeoTransform(im_geotrans)
dataset.SetProjection(im_proj)
if im_bands == 1:
dataset.GetRasterBand(1).WriteArray(im_data)
else:
for i in range(im_bands):
dataset.GetRasterBand(i+1).WriteArray(im_data[i])
del dataset
pre_path='dataset/pre/'
labellist = filter(lambda x: x.find('label')!=-1, os.listdir(pre_path))
list1 = list(map(lambda x: x[:], labellist))
label_name=pre_path + list1[0]
boundarylist = filter(lambda x: x.find('.shp')!=-1, os.listdir(pre_path+'boundary/'))
list2 = list(map(lambda x: x[:], boundarylist))
imagelist = filter(lambda x: x.find('tif')!=-1, os.listdir(pre_path+'img'))
list3 = list(map(lambda x: x[:], imagelist))
img_path=pre_path + list3[0]
"""
矢量裁剪
:param label_name: 要裁剪的矢量文件
:param boundary_name: 掩膜矢量文件
img_path: 影像
:param saveFolderPath: 裁剪后的矢量文件保存目录
:return:
"""
print('开始用矢量范围裁剪影像')
ogr.RegisterAll()
gdal.SetConfigOption("GDAL_FILENAME_IS_UTF8", "YES")
# 载入要裁剪的矢量文件
labelData = ogr.Open(label_name)
labelLayer = labelData.GetLayer( os.path.split( os.path.splitext( label_name )[0] )[1] )
spatial = labelLayer.GetSpatialRef()
geomType = labelLayer.GetGeomType()
# 载入掩膜矢量文件
def new_func(outLayerName):
return outLayerName
for i in list2:
boundary_name=pre_path+'boundary/'+ i
maskData = ogr.Open(boundary_name)
maskLayer = maskData.GetLayer()
#裁剪shp
# 生成裁剪后的矢量文件
save_shp_dir='./dataset/pre/shp/'
if not os.path.exists(save_shp_dir):
os.mkdir(save_shp_dir)
outLayerName = (save_shp_dir+i)
gdal.SetConfigOption("SHAPE_ENCODING", "GBK")
driver = ogr.GetDriverByName("ESRI Shapefile")
outData = driver.CreateDataSource(outLayerName)
outLayer = outData.CreateLayer(new_func(outLayerName), spatial, geomType)
labelLayer.Clip(maskLayer, outLayer)
lyr = maskData.GetLayer( os.path.split( os.path.splitext( boundary_name )[0] )[1] )
shpminX, shpmaxX, shpminY, shpmaxY = lyr.GetExtent()
#裁剪tif
flag=0
for j in list3:
raster_path = pre_path+'img/'+j
srcImage = gdal.Open(raster_path)
geocd = srcImage.GetGeoTransform()
geoProj = srcImage.GetProjection()
row = srcImage.RasterXSize
line = srcImage.RasterYSize
tifxmin = geocd[0]
tifymin = geocd[3]
tifxmax = geocd[0] + (row) * geocd[1] + (line) * geocd[2]
tifymax = geocd[3] + (row) * geocd[4] + (line) * geocd[5]
if shpminX>=tifxmin and shpmaxX<=tifxmax and shpminY<=tifymin and shpmaxY>=tifymax:
ulX, ulY = world2Pixel(geocd, shpminX, shpmaxY)
lrX, lrY = world2Pixel(geocd, shpmaxX, shpminY)
# Calculate the pixel size of the new image
pxWidth = int(lrX - ulX)
pxHeight = int(lrY - ulY)
clip = srcImage.ReadAsArray(ulX,ulY,pxWidth,pxHeight) #***只读要的那块***
xoffset = ulX
yoffset = ulY
geocd = list(geocd)
geocd[0] = shpminX
geocd[3] = shpmaxY
save_path='dataset/sat_train/'+i[:-4]+'.tif'
write_img(save_path, geoProj, geocd, clip)
gdal.ErrorReset()
outData.Release()
maskData.Release()
flag=1
if flag==0:
print(raster_path+"没有制作")
else:
print(raster_path)
labelData.Release()
文章浏览阅读3.8k次。1、将下载好的萤石js插件，添加到SoringBoot项目中。位置可参考下图所示。（容易出错的地方，在将js插件在html页面引入时，发生路径错误的问题）所以如果对页面中引入js的路径不清楚，可参考下图所示存放路径。2、将ezuikit.js引入到demo-live.html中。（可直接将如下代码复制到你创建的html页面中）<!DOCTYPE html><html lan..._ezuikit 测试的url
文章浏览阅读322次。第二步，在弹出的对话框选择，设备驱动—>PLC—>莫迪康—>ModbusRTU—>COM,根据配置软件选择的协议选期期，这里以此为例，然后点击“下一步”。第四步，把使用虚拟串口打勾(GPRS设备)，根据需要选择要生成虚拟口，这里以选择KVCOM1为例，然后点击“下一步”设备ID即Modbus地址(1-255) 使用DTU时，为下485接口上的设备地址。第六步，Modbus的从机地址，与配置软件相同，这里以1为例，点击“下一步“第五步，Modbus的从机地址，与配置软件相同，这里以1为例，点击“下一步“_组态王ua
文章浏览阅读9.4k次，点赞22次，收藏19次。安装npm相当于安装node.js,Node.js已自带npm，安装Node.js时会一起安装，npm的作用就是对Node.js依赖的包进行管理，也可以理解为用来安装/卸载Node.js需要装的东西_npm安装配置
文章浏览阅读748次，点赞21次，收藏26次。大家好，小编来为大家解答以下问题，python基础训练100题，python入门100例题，现在让我们一起来看看吧！宝子们还在新手村练级的时候，不单要吸入基础知识，夯实自己的理论基础，还要去实际操作练练手啊！由于文章篇幅限制，不可能将100道题全部呈现在此除了这些，下面还有我整理好的基础入门学习资料，视频和讲解文案都很齐全，用来入门绝对靠谱，需要的自提。保证100%免费这不，贴心的我爆肝给大家整理了这份今天给大家分享100道Python练习题。大家一定要给我三连啊~
文章浏览阅读1k次。 为了在 Linux （ Ubuntu） 上安装sublime，一般大家都会选择常见的教程或是 sublime 官网教程，然而在国内这种方法可能失效。为此，需要用安装包安装。以下就是使用官网安装包安装的教程。打开 sublime 官网后，点击右上角 download， 或是直接访问点击打开链接，即可看到各个平台上的安装包。选择 Linux 64 位版并下载。下载后，打开终端，进入安装..._ubuntu 安装sumlime text打不开
文章浏览阅读563次，点赞13次，收藏6次。CrossOver24是一款类虚拟机软件，专为macOS和Linux用户设计。它的核心技术是Wine，这是一种在Linux和macOS等非Windows操作系统上运行Windows应用程序的开源软件。通过CrossOver24，用户可以在不购买Windows授权或使用传统虚拟机的情况下，直接在Mac或Linux系统上运行Windows软件和游戏。该软件还提供了丰富的功能，如自动配置、无缝集成和实时传输等，以实现高效的跨平台操作体验。
文章浏览阅读1.7k次。一个用聊天的方式让ChatGPT帮我写的线程安全的环形List_为什么gpt一写list就卡
文章浏览阅读336次。我们在前面的文章里曾写过Web应用中乱码产生的原因和处理方式，旧文回顾：深度揭秘乱码问题背后的原因及解决方式其中我们提到可以通过Filter的方式来设置请求和响应的encoding，来解..._filterconfig selectencoding
文章浏览阅读651次。转自：http://www.jb51.net/article/36480.htmencodeURI和decodeURI是成对来使用的，因为浏览器的地址栏有中文字符的话，可以会出现不可预期的错误，所以可以encodeURI把非英文字符转化为英文编码，decodeURI可以用来把字符还原回来_js encodeur decodeurl
文章浏览阅读1.9w次，点赞6次，收藏3次。前言在日常的Android开发当中，我们肯定要打包apk。但是今天我打包的时候遇到一个很奇怪的问题Android The destination folder does not exist or is not writeable，大意是目标文件夹不存在或不可写。出现问题的原因以及解决办法上面有说报错的中文大意是：目标文件夹不存在或不可写。其实问题就在我们的打包界面当中图中标红的Desti..._the destination folder does not exist or is not writeable
文章浏览阅读94次。一、配置代码编辑区的样式 <1>打开Eclipse，Help —> Install NewSoftware，界面如下： <2>点击add...，按下图所示操作： name：随意填写，Location：http://eclipse-color-th..._ecplise高大上设置
文章浏览阅读2.8k次。一，下载mysql:http://dev.mysql.com/downloads/mysql/； 打开页面之后，在Select Platform:下选择linux Generic，如果没有出现Linux的选项，请换一个浏览器试试。我用的谷歌版本不可以，换一个别的浏览器就行了，如果还是不行，需要换一个翻墙的浏览器。 二，下载完后解压缩并放到安装文件夹下： 1、MySQL-client-5.6.2_linux mysql 安装 mysql-5.6.24-1.linux_glibc2.5.x86_64.rpm-bundle