技术标签： 原型模式
# -*- 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()
文章浏览阅读86次。首发地址：https://yq.aliyun.com/articles/221708谈到机器学习，相信很多除学者都是通过斯坦福大学吴恩达老师的公开课《Machine Learning》开始具体的接触机器学习这个领域，但是学完之后又不知道自己的掌握情况，缺少一些实际的项目操作。对于机器学习的相关竞赛挑战，有些项目的门槛有些高，参加后难以具体的实现，因此造..._scrath五子棋下载
文章浏览阅读83次。原标题：Oracle 12c新特性系列专题-安徽Oracle授权认证中心 随着Oracle database 12c的普及，数据库管理员 (DBA) 的角色也随之发生了转变。 Oracle 12c数据库对 DBA 而言是下一代数据管理。它让 DBA 可以摆脱单调的日常管理任务，能够专注于如何从数据中获取更多价值。未来我们会推出基于Oracle12c的技术文章，帮助DBA尽快掌握新一代数据库的新特性..._ilm add policy row store compress advanced row after
文章浏览阅读150次。问题及代码：*Copyright(c)2016,烟台大学计算机与控制工程学院 *All right reserved. *文件名称：负数把正数赶出队列.cpp *作者：张冰 *完成日期；2016年10月09日 *版本号；v1.0 * *问题描述： 设从键盘输入一整数序列a1，a2，…an，试编程实现： 当ai>0时，ai进队，当ai<0时，将队首元素出队，当ai
文章浏览阅读376次。贪心+构造
文章浏览阅读150次。本文讲的是Linux命名空间学习教程（二） IPC，【编者的话】Docker核心解决的问题是利用LXC来实现类似VM的功能，从而利用更加节省的硬件资源提供给用户更多的计算资源。而 LXC所实现的隔离性主要是来自内核的命名空间, 其中pid、net、ipc、mnt、uts 等命名空间将容器的进程、网络、消息、文件系统和hostname 隔离开。本文是Li..._主机的 ipc 命名空间
文章浏览阅读2w次，点赞5次，收藏7次。在设备上强制安装apk。在app已有的情况下使用-r参数在app版本低于现有版本使用-d参数命令adb install -r -d xxx.apk_adb绕过安装程序强制安装app
文章浏览阅读290次。如果是越界进入硬件错误中断，MSP 或者 PSP 保存错误地址，跳转前会保存上一次执行的地址，lr 寄存器会保存子函数的地址，所以如果在 HardFault_CallBack 中直接调用 C 语言函数接口会间接修改了 lr，为了解决这个问题，直接绕过 lr 的 C 语言代码，用汇编语言提取 lr 寄存器再决定后面的操作。由于 STM32 加入了 FreeRTOS 操作系统，可能导致无法准确定位，仅供参考（日常编程需要考虑程序的健壮性，特别是对数组的访问，非常容易出现越界的情况）。_stm32flash地址越界怎么解决
文章浏览阅读1.8k次。学到了一种操作，说实话，我从来没想过还能这样正常情况下，为了管理方便，许多管理员都会开放MySQL数据库的secure_file_priv，这时就可以导入或者导出数据当我如图输入时，就会在D盘创建一个名为123456.php，内容为<?php phpinfo();?>的文件我们可以利用这一点运用到SQL注入中，从拿下数据库到拿下目标的服务器比如我们在使用联合查询注入，正常是这样的语句http://xxx?id=-1 union select 1,'你想知道的字段的内容或查询语句',
文章浏览阅读2.9w次，点赞12次，收藏63次。感谢原文：https://blog.csdn.net/abc5382334/article/details/24260817感谢原文：https://blog.csdn.net/jiaqingge/article/details/52564348Html CSS的三种链接方式css文本的链接方式有三种：分别是内联定义、链入内部css、和链入外部css1.代码为：<html>..._html链接css代码
文章浏览阅读625次。近几年，蓝牙耳机市场发展迅速，越来越多的消费者希望抛弃线缆，更自由地听音乐，对于运动人士来说，蓝牙耳机的便携性显得尤为重要。但目前市面上的大多数蓝牙耳机实际上都是“有线”的，运动过程中产生的听诊器效应会严重影响听歌的感受。而在“真无线”耳机领域，除了苹果的AirPods外，可供选择的产品并不多，而AirPods又不是为运动场景打造的，防水能力非常差。那么对于喜欢运动又想要“自由”的朋友来说，有没有一款产品能够满足他们的需求呢？下面这十款小编专门为大家搜罗的蓝牙耳机或许就能找到适合的！网红击音F1_适合游戏与运动的高音质蓝牙耳机
文章浏览阅读1k次，点赞6次，收藏7次。在本篇博文中，我们在 iOS 17 beta 4（SwiftUI 5.0）测试版中发现了 SwiftUI 视图首次显示时状态的改变会导致动画“副作用”的问题，并提出多种解决方案。
文章浏览阅读1.9k次。 在 上篇文章–Flutter 实现支持上拉加载和下拉刷新的 ListView 中，我们最终实现的效果是在 listView 上面留下了一段空白，本意是用来加载轮播图的，于是今天就开发了一下，希望能给各位灵感。一 、效果如下说一下大体思路 其实图片展示是用的 PageView ，然后，下面的指示器 是用的 TabPageSelector ，当然整体是用 Stack 包裹起来的。1、..._flutter pageview轮播图 site:csdn.net