pythonqgispyqgis

Upgrading Windows yeilds empty shape files


I have a script I have been running in Windows 10 for quite some time and it has been working flawlessly. What it does is irrelevant, you can probably figure out by reading the code, but to conclude it it would be to say that it extract certain slopes in certain aspects.

After a major power out the main drive crashed and I had to reinstall Windows, I thought upgrading to Windows 11 would be a good idea. However, now that I am running the script on Windows 11 the shapefiles produced are empty. There are no errors reported.

Script:

from os.path import join, normpath
import os, shutil
from datetime import datetime
import time
from pushbullet import Pushbullet

print("Start point:", datetime.now())

pb = Pushbullet("xx")

path_start = normpath(r'H:\\GIS\\worklayers0\\')
path_temp = normpath(r'H:\\GIS\\templayers0\\')
path_temp2 = normpath(r'H:\\GIS\\temp2layers0\\')
path_end = normpath(r'H:\\GIS\\resultlayers0\\')

slope_list = [4, 4.9, 5.9, 7.1, 8.2, 10, 11.8, 13.8, 15.7, 17.9, 20.6, 25]
aspect_list_1 = [[0, 15], [15,30], [30, 45], [45, 60], [60, 75], [75, 90], [90, 105], [105, 120], [120, 135], [135, 150], [150, 165], [165, 180]]
aspect_list_2 = [[180, 195], [195, 210], [210, 225], [225, 240], [240, 255], [255, 270], [270, 285], [285, 300], [300, 315], [315, 330], [330, 345], [345, 360]]

for filename in os.listdir(path_temp):
    file_path = os.path.join(path_temp, filename)
    try:
        if os.path.isfile(file_path) or os.path.islink(file_path):
            os.unlink(file_path)
        elif os.path.isdir(file_path):
            shutil.rmtree(file_path)
    except:
        pass

for filename in os.listdir(path_temp2):
    file_path = os.path.join(path_temp2, filename)
    try:
        if os.path.isfile(file_path) or os.path.islink(file_path):
            os.unlink(file_path)
        elif os.path.isdir(file_path):
            shutil.rmtree(file_path)
    except:
        pass
    
for filename in os.listdir(path_end):
    file_path = os.path.join(path_end, filename)
    try:
        if os.path.isfile(file_path) or os.path.islink(file_path):
            os.unlink(file_path)
        elif os.path.isdir(file_path):
            shutil.rmtree(file_path)
    except:
        pass

files = os.listdir(path_start)

for index, file in enumerate(files):
    os.rename(os.path.join(path_start, file), os.path.join(path_start, ''.join([str(index), '.tif'])))

for number in range(3):
    
    processing.run("native:slope", {'INPUT':path_start+'\\'+str(number)+'.tif','Z_FACTOR':1,'OUTPUT':path_temp+'\\'+'slope_'+str(number)+'.tif'})
    processing.run("native:aspect", {'INPUT':path_start+'\\'+str(number)+'.tif','Z_FACTOR':1,'OUTPUT':path_temp+'\\'+'aspect_'+str(number)+'.tif'})
    
    processing.run("gdal:merge", {'INPUT':[path_temp+'\\'+'aspect_'+str(number)+'.tif',path_temp+'\\'+'slope_'+str(number)+'.tif'],'PCT':True,'SEPARATE':True,'NODATA_INPUT':None,'NODATA_OUTPUT':None,'OPTIONS':'','EXTRA':'','DATA_TYPE':5,'OUTPUT':path_temp+'\\'+'merge_'+str(number)+'.tif'})
        
    processing.run("native:pixelstopolygons", {'INPUT_RASTER':path_temp+'\\'+'merge_'+str(number)+'.tif','RASTER_BAND':1,'FIELD_NAME':'ASPECT','OUTPUT':path_temp+'\\'+'aspect_'+str(number)+'.shp'})
    processing.run("native:createspatialindex", {'INPUT':path_temp+'\\'+'aspect_'+str(number)+'.shp'})
        
    processing.run("native:pixelstopolygons", {'INPUT_RASTER':path_temp+'\\'+'merge_'+str(number)+'.tif','RASTER_BAND':2,'FIELD_NAME':'SLOPE','OUTPUT':path_temp+'\\'+'slope_'+str(number)+'.shp'})
    processing.run("native:createspatialindex", {'INPUT':path_temp+'\\'+'slope_'+str(number)+'.shp'})
        
    processing.run("native:joinattributesbylocation", {'INPUT':path_temp+'\\'+'slope_'+str(number)+'.shp','PREDICATE':[5],'JOIN':path_temp+'\\'+'aspect_'+str(number)+'.shp','JOIN_FIELDS':[],'METHOD':0,'DISCARD_NONMATCHING':False,'PREFIX':'','OUTPUT':path_temp+'\\'+'aspect_slope_'+str(number)+'.shp'})
    aspect_slope_layer = iface.addVectorLayer(path_temp+'\\'+'aspect_slope_'+str(number)+'.shp', "", "ogr")
        
    
    for i in range(len(slope_list)):
        aspect_list_1[i].append(slope_list[i])
    
    for i in range(len(slope_list)):
        aspect_list_2[i].append(slope_list[::-1][i])
        
    for aspect_interval in aspect_list_1:
            start = aspect_interval[0] 
            end = aspect_interval[1]
            slope_loop = aspect_interval[2]
            aspect_slope_layer.selectByExpression('"ASPECT">'+str(start)+' and "ASPECT"<='+str(end)+' and "SLOPE">='+str(slope_loop))
            QgsVectorFileWriter.writeAsVectorFormat(aspect_slope_layer, str(path_temp)+'\\'+'aspect_slope_'+str(start)+'-'+str(end)+'_'+str(number)+'.shp', "UTF-8", aspect_slope_layer.crs(), "ESRI Shapefile", onlySelected=True)
            processing.run("native:createspatialindex", {'INPUT':str(path_temp)+'\\'+'aspect_slope_'+str(start)+'-'+str(end)+'_'+str(number)+'.shp'})
            
    for aspect_interval in aspect_list_2:
            start = aspect_interval[0] 
            end = aspect_interval[1]
            slope_loop = aspect_interval[2]
            aspect_slope_layer.selectByExpression('"ASPECT">'+str(start)+' and "ASPECT"<='+str(end)+' and "SLOPE">='+str(slope_loop))
            QgsVectorFileWriter.writeAsVectorFormat(aspect_slope_layer, str(path_temp)+'\\'+'aspect_slope_'+str(start)+'-'+str(end)+'_'+str(number)+'.shp', "UTF-8", aspect_slope_layer.crs(), "ESRI Shapefile", onlySelected=True)
            processing.run("native:createspatialindex", {'INPUT':str(path_temp)+'\\'+'aspect_slope_'+str(start)+'-'+str(end)+'_'+str(number)+'.shp'})
    
    processing.run("native:mergevectorlayers", {'LAYERS':[path_temp+'\\'+'aspect_slope_0-15_'+str(number)+'.shp',path_temp+'\\'+'aspect_slope_345-360_'+str(number)+'.shp'],'CRS':None,'OUTPUT':str(path_temp2)+'\\'+'1_aspect_slope_int_'+str(number)+'.shp'})
    processing.run("native:mergevectorlayers", {'LAYERS':[path_temp+'\\'+'aspect_slope_15-30_'+str(number)+'.shp',path_temp+'\\'+'aspect_slope_330-345_'+str(number)+'.shp'],'CRS':None,'OUTPUT':str(path_temp2)+'\\'+'2_aspect_slope_int_'+str(number)+'.shp'})
    processing.run("native:mergevectorlayers", {'LAYERS':[path_temp+'\\'+'aspect_slope_30-45_'+str(number)+'.shp',path_temp+'\\'+'aspect_slope_315-330_'+str(number)+'.shp'],'CRS':None,'OUTPUT':str(path_temp2)+'\\'+'3_aspect_slope_int_'+str(number)+'.shp'})
    processing.run("native:mergevectorlayers", {'LAYERS':[path_temp+'\\'+'aspect_slope_45-60_'+str(number)+'.shp',path_temp+'\\'+'aspect_slope_300-315_'+str(number)+'.shp'],'CRS':None,'OUTPUT':str(path_temp2)+'\\'+'4_aspect_slope_int_'+str(number)+'.shp'})
    processing.run("native:mergevectorlayers", {'LAYERS':[path_temp+'\\'+'aspect_slope_60-75_'+str(number)+'.shp',path_temp+'\\'+'aspect_slope_285-300_'+str(number)+'.shp'],'CRS':None,'OUTPUT':str(path_temp2)+'\\'+'5_aspect_slope_int_'+str(number)+'.shp'})
    processing.run("native:mergevectorlayers", {'LAYERS':[path_temp+'\\'+'aspect_slope_75-90_'+str(number)+'.shp',path_temp+'\\'+'aspect_slope_270-285_'+str(number)+'.shp'],'CRS':None,'OUTPUT':str(path_temp2)+'\\'+'6_aspect_slope_int_'+str(number)+'.shp'})
    processing.run("native:mergevectorlayers", {'LAYERS':[path_temp+'\\'+'aspect_slope_90-105_'+str(number)+'.shp',path_temp+'\\'+'aspect_slope_255-270_'+str(number)+'.shp'],'CRS':None,'OUTPUT':str(path_temp2)+'\\'+'7_aspect_slope_int_'+str(number)+'.shp'})
    processing.run("native:mergevectorlayers", {'LAYERS':[path_temp+'\\'+'aspect_slope_105-120_'+str(number)+'.shp',path_temp+'\\'+'aspect_slope_240-255_'+str(number)+'.shp'],'CRS':None,'OUTPUT':str(path_temp2)+'\\'+'8_aspect_slope_int_'+str(number)+'.shp'})
    processing.run("native:mergevectorlayers", {'LAYERS':[path_temp+'\\'+'aspect_slope_120-135_'+str(number)+'.shp',path_temp+'\\'+'aspect_slope_225-240_'+str(number)+'.shp'],'CRS':None,'OUTPUT':str(path_temp2)+'\\'+'9_aspect_slope_int_'+str(number)+'.shp'})
    processing.run("native:mergevectorlayers", {'LAYERS':[path_temp+'\\'+'aspect_slope_135-150_'+str(number)+'.shp',path_temp+'\\'+'aspect_slope_210-225_'+str(number)+'.shp'],'CRS':None,'OUTPUT':str(path_temp2)+'\\'+'10_aspect_slope_int_'+str(number)+'.shp'})
    processing.run("native:mergevectorlayers", {'LAYERS':[path_temp+'\\'+'aspect_slope_150-165_'+str(number)+'.shp',path_temp+'\\'+'aspect_slope_195-210_'+str(number)+'.shp'],'CRS':None,'OUTPUT':str(path_temp2)+'\\'+'11_aspect_slope_int_'+str(number)+'.shp'})
    processing.run("native:mergevectorlayers", {'LAYERS':[path_temp+'\\'+'aspect_slope_165-180_'+str(number)+'.shp',path_temp+'\\'+'aspect_slope_180-195_'+str(number)+'.shp'],'CRS':None,'OUTPUT':str(path_temp2)+'\\'+'12_aspect_slope_int_'+str(number)+'.shp'})


    for i in range(1,13):
        processing.run("native:createspatialindex", {'INPUT':str(path_temp2)+'\\'+str(i)+'_aspect_slope_int_'+str(number)+'.shp'})     
    
    QgsProject.instance().removeAllMapLayers()

merge_list = os.listdir(path_temp2)

for i in range(1,13):
    check = str(i)+'_'
    name_ext_matches = []
    name_matches = [idx for idx in merge_list if idx.lower().startswith(check.lower())]
    
    for file in name_matches:
        
        if file.endswith('.shp'):
            name_ext_matches.append(file)
            
        matches = [path_temp2 +'\\'+ x for x in name_ext_matches]
            
    processing.run("native:mergevectorlayers", {'LAYERS':matches,'CRS':None,'OUTPUT':str(path_temp)+'\\'+str(i)+'_aspect_slope_int_int.shp'})
    processing.run("native:createspatialindex", {'INPUT':str(path_temp)+'\\'+str(i)+'_aspect_slope_int_int.shp'})    
    processing.run("native:refactorfields", {'INPUT':str(path_temp)+'\\'+str(i)+'_aspect_slope_int_int.shp','FIELDS_MAPPING':[{'expression': '"SLOPE"','length': 20,'name': 'SLOPE','precision': 8,'sub_type': 0,'type': 2,'type_name': 'integer'},{'expression': '"ASPECT"','length': 20,'name': 'ASPECT','precision': 8,'sub_type': 0,'type': 4,'type_name': 'int8'}],'OUTPUT':str(path_end)+'\\'+str(i)+'_aspect_slope.shp'})
    processing.run("native:createspatialindex", {'INPUT':str(path_end)+'\\'+str(i)+'_aspect_slope.shp'}) 
  

phone = pb.devices[0]
push = pb.push_note("Wakie wakie!", "All done with extracting the slopes.")
print("Notifications sent!")
print("Done!")
print("End point:", datetime.now())

Solution

  • I finally found what was wrong. It was the version that didn't support my script. I was using 3.26.2 before upgrading to Win11. Running on Win11 I have been using 3.28. Could there be some change in the GDAL which shows in the scripts? In any case, it's weird since it raises no errors.