External vector processing

spacenetutilities/geoTools.py

@@ -428,25 +428,44 @@ def createPolygonFromCorners(left,bottom,right, top):
     return poly
 
 
-def clipShapeFile(geoDF, outputFileName, polyToCut, minpartialPerc=0.0, shapeLabel='Geo', debug=False):
+def clipShapeFile(geoDF, outputFileName, polyToCut, minpartialPerc=0.0, geomType="Polygon", shapeLabel='Geo', debug=False):
     # check if geoDF has origAreaField
     outGeoJSon = os.path.splitext(outputFileName)[0] + '.geojson'
     if not os.path.exists(os.path.dirname(outGeoJSon)):
         os.makedirs(os.path.dirname(outGeoJSon))
     if debug:
         print(outGeoJSon)
 
-    if 'origarea' in geoDF.columns:
-        pass
-    else:
-        geoDF['origarea'] = geoDF.area
 
     #TODO must implement different case for lines and for spatialIndex
+    if geomType=="LineString":
+        if 'origarea' in geoDF.columns:
+            pass
+        else:
+            geoDF['origarea'] = 0
+        if 'origlen' in geoDF.columns:
+            pass
+        else:
+            geoDF['origlen'] = geoDF.length
+
 
+        cutGeoDF = geoDF.loc[geoDF.intersection(polyToCut).length / geoDF['origlen'] > minpartialPerc].copy()
+        cutGeoDF['partialDec'] = cutGeoDF.length / cutGeoDF['origlen']
+        cutGeoDF['truncated'] = (cutGeoDF['partialDec'] == 1.0).astype(int)
+
+    else:
+        if 'origarea' in geoDF.columns:
+            pass
+        else:
+            geoDF['origarea'] = geoDF.area
+        if 'origlen' in geoDF.columns:
+            pass
+        else:
+            geoDF['origlen'] = 0
 
-    cutGeoDF = geoDF.loc[geoDF.intersection(polyToCut).area/geoDF['origarea'] > minpartialPerc].copy()
-    cutGeoDF['partialDec'] = cutGeoDF.area/cutGeoDF['origarea']
-    cutGeoDF['truncated'] = (cutGeoDF['partialDec']==1.0).astype(int)
+        cutGeoDF = geoDF.loc[geoDF.intersection(polyToCut).area/geoDF['origarea'] > minpartialPerc].copy()
+        cutGeoDF['partialDec'] = cutGeoDF.area/cutGeoDF['origarea']
+        cutGeoDF['truncated'] = (cutGeoDF['partialDec']==1.0).astype(int)
 
     if cutGeoDF.empty:
         with open(outGeoJSon, 'a'):

spacenetutilities/scripts/externalVectorProcessing.py

@@ -1,64 +1,56 @@
 from spacenetutilities import geoTools as gT
 from spacenetutilities import labeltools as lT
-from osgeo import gdal, ogr, osr
+import fiona as fio
+import geopandas as gpd
+import rasterio
 import argparse
 import os
 import glob
 
 
 
 def buildTindex(rasterFolder, rasterExtention='.tif'):
+
     rasterList = glob.glob(os.path.join(rasterFolder, '*{}'.format(rasterExtention)))
     print(rasterList)
 
     print(os.path.join(rasterFolder, '*{}'.format(rasterExtention)))
 
-    memDriver = ogr.GetDriverByName('MEMORY')
-    gTindex = memDriver.CreateDataSource('gTindex')
-    srcImage = gdal.Open(rasterList[0])
-    spat_ref = osr.SpatialReference()
-    spat_ref.SetProjection(srcImage.GetProjection())
-    gTindexLayer = gTindex.CreateLayer("gtindexlayer", spat_ref, geom_type=ogr.wkbPolygon)
 
-    # Add an ID field
-    idField = ogr.FieldDefn("location", ogr.OFTString)
-    gTindexLayer.CreateField(idField)
 
-    # Create the feature and set values
-    featureDefn = gTindexLayer.GetLayerDefn()
+    featureList = []
 
+    for rasterFile in rasterList:
+        with rasterio.open(rasterFile) as srcImage:
 
 
-    for rasterFile in rasterList:
-        srcImage = gdal.Open(rasterFile)
+            geoTrans, polyToCut, ulX, ulY, lrX, lrY = gT.getRasterExtent(srcImage)
 
-        geoTrans, polyToCut, ulX, ulY, lrX, lrY = gT.getRasterExtent(srcImage)
+        feature = {'geometry': polyToCut,
+                   'location': rasterFile}
 
-        feature = ogr.Feature(featureDefn)
-        feature.SetGeometry(polyToCut)
-        feature.SetField("location", rasterFile)
-        gTindexLayer.CreateFeature(feature)
-        feature = None
+        featureList.append(feature)
 
+        geoDFTindex = gpd.GeoDataFrame(featureList)
 
-    return gTindex, gTindexLayer
+    return geoDFTindex
 
 
 def createTiledGeoJsonFromSrc(rasterFolderLocation, vectorSrcFile, geoJsonOutputDirectory, rasterTileIndex='',
                               geoJsonPrefix='GEO', rasterFileExtenstion='.tif',
                               rasterPrefixToReplace='PAN'
                               ):
     if rasterTileIndex == '':
-        gTindex, gTindexLayer = buildTindex(rasterFolderLocation, rasterExtention=rasterFileExtenstion)
+        geoDFTindex = buildTindex(rasterFolderLocation, rasterExtention=rasterFileExtenstion)
     else:
-        gTindex = ogr.Open(rasterTileIndex,0)
-        gTindexLayer = gTindex.GetLayer()
+        geoDFTindex = gpd.read_file(rasterTileIndex)
+
+    shapeSrc = gpd.read_file(vectorSrcFile)
 
-    shapeSrc = ogr.Open(vectorSrcFile,0)
     chipSummaryList = []
-    for feature in gTindexLayer:
-        featureGeom = feature.GetGeometryRef()
-        rasterFileName = feature.GetField('location')
+    for idx, feature in geoDFTindex.iterrows():
+        featureGeom = feature['geometry']
+        rasterFileName = feature['location']
         rasterFileBaseName = os.path.basename(rasterFileName)
         outGeoJson = rasterFileBaseName.replace(rasterPrefixToReplace, geoJsonPrefix)
         outGeoJson = outGeoJson.replace(rasterFileExtenstion, '.geojson')
@@ -127,6 +119,7 @@ def createTiledGeoJsonFromSrc(rasterFolderLocation, vectorSrcFile, geoJsonOutput
                               )
 
 
+
     outputCSVFileName = geoJsonOutputDirectory+"OSM_Proposal.csv"
     lT.createCSVSummaryFile(chipSummaryList, outputCSVFileName,
                                 replaceImageID=rasterPrefix+"_",