def get_modis_tile_list(ds):
"""Helper function to identify MODIS tiles that intersect input geometry
modis_gird.py contains dictionary of tile boundaries (tile name and WKT polygon ring from bbox)
See: https://modis-land.gsfc.nasa.gov/MODLAND_grid.html
"""
from demcoreg import modis_grid
modis_dict = modis_grid.modis_dict
for key in modis_dict:
modis_dict[key] = ogr.CreateGeometryFromWkt(modis_dict[key])
geom = geolib.ds_geom(ds)
geom_dup = geolib.geom_dup(geom)
ct = osr.CoordinateTransformation(geom_dup.GetSpatialReference(), geolib.wgs_srs)
geom_dup.Transform(ct)
tile_list = []
for key, val in list(modis_dict.items()):
if geom_dup.Intersects(val):
tile_list.append(key)
return tile_list
python类CreateGeometryFromWkt()的实例源码
def shp2geom(shp_fn):
"""Extract geometries from input shapefile
Need to handle multi-part geom: http://osgeo-org.1560.x6.nabble.com/Multipart-to-singlepart-td3746767.html
"""
ds = ogr.Open(shp_fn)
lyr = ds.GetLayer()
srs = lyr.GetSpatialRef()
lyr.ResetReading()
geom_list = []
for feat in lyr:
geom = feat.GetGeometryRef()
geom.AssignSpatialReference(srs)
#Duplicate the geometry, or segfault
#See: http://trac.osgeo.org/gdal/wiki/PythonGotchas
#g = ogr.CreateGeometryFromWkt(geom.ExportToWkt())
#g.AssignSpatialReference(srs)
g = geom_dup(geom)
geom_list.append(g)
#geom = ogr.ForceToPolygon(' '.join(geom_list))
#Dissolve should convert multipolygon to single polygon
#return geom_list[0]
ds = None
return geom_list
def ds_geom(ds, t_srs=None):
"""Return dataset bbox envelope as geom
"""
gt = ds.GetGeoTransform()
ds_srs = get_ds_srs(ds)
if t_srs is None:
t_srs = ds_srs
ns = ds.RasterXSize
nl = ds.RasterYSize
x = np.array([0, ns, ns, 0, 0], dtype=float)
y = np.array([0, 0, nl, nl, 0], dtype=float)
#Note: pixelToMap adds 0.5 to input coords, need to account for this here
x -= 0.5
y -= 0.5
mx, my = pixelToMap(x, y, gt)
geom_wkt = 'POLYGON(({0}))'.format(', '.join(['{0} {1}'.format(*a) for a in zip(mx,my)]))
geom = ogr.CreateGeometryFromWkt(geom_wkt)
geom.AssignSpatialReference(ds_srs)
if not ds_srs.IsSame(t_srs):
geom_transform(geom, t_srs)
return geom
def get_features_as_geojson(layer, bbox=None, union=False):
""" Get features in this layer and return as GeoJSON """
if bbox is not None:
layer.SetSpatialFilterRect(bbox[0], bbox[3], bbox[2], bbox[1])
poly = ogr.Geometry(ogr.wkbPolygon)
if union:
for feature in layer:
geom = feature.GetGeometryRef()
# required for ogr2
if hasattr(geom, 'GetLinearGeometry'):
geom = geom.GetLinearGeometry()
poly = poly.Union(geom)
if bbox is not None:
wkt = "POLYGON ((%s %s, %s %s, %s %s, %s %s, %s %s))" % \
(bbox[0], bbox[1], bbox[2], bbox[1], bbox[2], bbox[3], bbox[0], bbox[3], bbox[0], bbox[1])
bbox_wkt = ogr.CreateGeometryFromWkt(wkt)
poly = poly.Intersection(bbox_wkt)
return json.loads(poly.ExportToJson())
def get_idx_as_shp(self, path, gt, wkt):
'''
Exports a Shapefile containing the locations of the extracted
endmembers. Assumes the coordinates are in decimal degrees.
'''
coords = pixel_to_xy(self.get_idx(), gt=gt, wkt=wkt, dd=True)
driver = ogr.GetDriverByName('ESRI Shapefile')
ds = driver.CreateDataSource(path)
srs = osr.SpatialReference()
srs.ImportFromEPSG(4326)
layer = ds.CreateLayer(path.split('.')[0], srs, ogr.wkbPoint)
for pair in coords:
feature = ogr.Feature(layer.GetLayerDefn())
# Create the point from the Well Known Text
point = ogr.CreateGeometryFromWkt('POINT(%f %f)' % pair)
feature.SetGeometry(point) # Set the feature geometry
layer.CreateFeature(feature) # Create the feature in the layer
feature.Destroy() # Destroy the feature to free resources
# Destroy the data source to free resources
ds.Destroy()
def readwktcsv(csv_path,removeNoBuildings=True, groundTruthFile=True):
#
# csv Format Expected = ['ImageId', 'BuildingId', 'PolygonWKT_Pix', 'PolygonWKT_Geo']
# returns list of Dictionaries {'ImageId': image_id, 'BuildingId': building_id, 'poly': poly}
# image_id is a string,
# BuildingId is an integer,
# poly is a ogr.Geometry Polygon
buildinglist = []
with open(csv_path, 'rb') as csvfile:
building_reader = csv.reader(csvfile, delimiter=',', quotechar='"')
next(building_reader, None) # skip the headers
for row in building_reader:
if removeNoBuildings:
if int(row[1]) != -1:
polyPix = ogr.CreateGeometryFromWkt(row[2])
if groundTruthFile:
polyGeo = ogr.CreateGeometryFromWkt(row[3])
else:
polyGeo = []
buildinglist.append({'ImageId': row[0], 'BuildingId': int(row[1]), 'polyPix': polyPix,
'polyGeo': polyGeo,
})
else:
polyPix = ogr.CreateGeometryFromWkt(row[2])
if groundTruthFile:
polyGeo = ogr.CreateGeometryFromWkt(row[3])
else:
polyGeo = []
buildinglist.append({'ImageId': row[0], 'BuildingId': int(row[1]), 'polyPix': polyPix,
'polyGeo': polyGeo,
})
return buildinglist
def pixelWKTToGeoWKT(geomWKT, inputRaster, targetSR='', geomTransform='', breakMultiPolygonPix=False):
# Returns GeoCoordinateList from PixelCoordinateList
geomPix = ogr.CreateGeometryFromWkt(geomWKT)
geomGeoList = pixelGeomToGeoGeom(geomPix, inputRaster, targetSR=targetSR,
geomTransform=geomTransform, breakMultiPolygonPix=breakMultiPolygonPix)
return geomGeoList
def geom_dup(geom):
"""Create duplicate geometry
Needed to avoid segfault when passing geom around. See: http://trac.osgeo.org/gdal/wiki/PythonGotchas
"""
g = ogr.CreateGeometryFromWkt(geom.ExportToWkt())
g.AssignSpatialReference(geom.GetSpatialReference())
return g
#This should be a function of a new geom class
#Assumes geom has srs defined
#Modifies geom in place
def bbox2geom(bbox, t_srs=None):
#Check bbox
#bbox = numpy.array([-180, 180, 60, 90])
x = [bbox[0], bbox[0], bbox[1], bbox[1], bbox[0]]
y = [bbox[2], bbox[3], bbox[3], bbox[2], bbox[2]]
geom_wkt = 'POLYGON(({0}))'.format(', '.join(['{0} {1}'.format(*a) for a in zip(x,y)]))
geom = ogr.CreateGeometryFromWkt(geom_wkt)
if t_srs is not None and not wgs_srs.IsSame(t_srs):
ct = osr.CoordinateTransformation(t_srs, wgs_srs)
geom.Transform(ct)
geom.AssignSpatialReference(t_srs)
return geom
def intersection(self,filename):
""" Return coordinates of intersection between this image and another.
:param filename : path to the second image (or another GeoRaster instance)
:type filename: str, georaster.__Raster
:returns: extent of the intersection between the 2 images \
(xmin, xmax, ymin, ymax)
:rtype: tuple
"""
# Create a polygon of the envelope of the first image
xmin, xmax, ymin, ymax = self.extent
wkt = "POLYGON ((%f %f, %f %f, %f %f, %f %f, %f %f))" \
%(xmin,ymin,xmin,ymax,xmax,ymax,xmax,ymin,xmin,ymin)
poly1 = ogr.CreateGeometryFromWkt(wkt)
# Create a polygon of the envelope of the second image
img = SingleBandRaster(filename, load_data=False)
xmin, xmax, ymin, ymax = img.extent
wkt = "POLYGON ((%f %f, %f %f, %f %f, %f %f, %f %f))" \
%(xmin,ymin,xmin,ymax,xmax,ymax,xmax,ymin,xmin,ymin)
poly2 = ogr.CreateGeometryFromWkt(wkt)
# Compute intersection envelope
intersect = poly1.Intersection(poly2)
extent = intersect.GetEnvelope()
# check that intersection is not void
if intersect.GetArea()==0:
print('Warning: Intersection is void')
return 0
else:
return extent
def has_geos():
pnt1 = ogr.CreateGeometryFromWkt('POINT(10 20)')
pnt2 = ogr.CreateGeometryFromWkt('POINT(30 20)')
ogrex = ogr.GetUseExceptions()
ogr.DontUseExceptions()
hasgeos = pnt1.Union(pnt2) is not None
if ogrex:
ogr.UseExceptions()
return hasgeos
def copy_features(src_layer, dst_layer, fix_geometry, simplify_geometry, start_index, total):
index = 0
batch_size = 200
index_batch = 0
for feat in src_layer:
if index < start_index:
index = index + 1
continue
try:
geom = shapely.wkt.loads(feat.GetGeometryRef().ExportToWkt())
except Exception as e:
print('Error({0}), skipping geometry.'.format(e))
continue
if fix_geometry and not geom.is_valid:
geom = geom.buffer(0.0)
if simplify_geometry:
geom = geom.simplify(0.004)
f = ogr.Feature(dst_layer.GetLayerDefn())
# set field values
for i in range(feat.GetFieldCount()):
fd = feat.GetFieldDefnRef(i)
f.SetField(fd.GetName(), feat.GetField(fd.GetName()))
# set geometry
f.SetGeometry(ogr.CreateGeometryFromWkt(geom.to_wkt()))
if index_batch == 0:
dst_layer.StartTransaction()
# create feature
feature = dst_layer.CreateFeature(f)
f.Destroy()
index_batch = index_batch + 1
if index_batch >= batch_size or index == total - 1:
dst_layer.CommitTransaction()
count = dst_layer.GetFeatureCount() # update number of inserted features
print('Inserted {0} of {1} features ({2:.2f}%)'.format(count, total, 100. * float(count) / total))
index_batch = 0
if index == total - 1:
break
index = index + 1
def shp2array(shp_fn, r_ds=None, res=None, extent=None, t_srs=None):
"""Rasterize input shapefile to match existing raster Dataset (or specified res/extent/t_srs)
"""
shp_ds = ogr.Open(shp_fn)
lyr = shp_ds.GetLayer()
#This returns xmin, ymin, xmax, ymax
shp_extent = lyr_extent(lyr)
shp_srs = lyr.GetSpatialRef()
# dst_dt = gdal.GDT_Byte
ndv = 0
if r_ds is not None:
r_extent = ds_extent(r_ds)
res = get_res(r_ds, square=True)[0]
if extent is None:
extent = r_extent
r_srs = get_ds_srs(r_ds)
r_geom = ds_geom(r_ds)
# dst_ns = r_ds.RasterXSize
# dst_nl = r_ds.RasterYSize
#Convert raster extent to shp_srs
cT = osr.CoordinateTransformation(r_srs, shp_srs)
r_geom_reproj = geom_dup(r_geom)
r_geom_reproj.Transform(cT)
r_geom_reproj.AssignSpatialReference(t_srs)
lyr.SetSpatialFilter(r_geom_reproj)
#lyr.SetSpatialFilter(ogr.CreateGeometryFromWkt(wkt))
else:
#TODO: clean this up
if res is None:
sys.exit("Must specify input res")
if extent is None:
print("Using input shp extent")
extent = shp_extent
if t_srs is None:
t_srs = r_srs
if not shp_srs.IsSame(t_srs):
print("Input shp srs: %s" % shp_srs.ExportToProj4())
print("Specified output srs: %s" % t_srs.ExportToProj4())
out_ds = lyr_proj(lyr, t_srs)
outlyr = out_ds.GetLayer()
else:
outlyr = lyr
#outlyr.SetSpatialFilter(r_geom)
m_ds = mem_ds(res, extent, srs=t_srs, dtype=gdal.GDT_Byte)
b = m_ds.GetRasterBand(1)
b.SetNoDataValue(ndv)
gdal.RasterizeLayer(m_ds, [1], outlyr, burn_values=[1])
a = b.ReadAsArray()
a = ~(a.astype('Bool'))
return a
