#===================raster sampling===============================
#get raster information at specific points

from qgis.core import *
from qgis import processing
import os.path, glob

##input variables--------------------------------------------------
pname="Stations_COSMO"
path = 'Z:/burc/COSMO/20200803_10/*.nc' #Change path to infiles
outfolder="Z:/burc/GIS/COSMO/CO2_vertical/"
var="CO2_RA"

##get raster layer and point file----------------------------------
#point file from map:
pfile = QgsProject.instance().mapLayersByName(pname)[0]
if not pfile.isValid():
    print("ERROR: Point file not found!")
#raster layer in folder:
fnames=[]
for file in glob.glob(path): # path to files
    uri = 'NETCDF:"' + file +'":'+ var
    #filename = str(int(''.join(filter(str.isdigit, file))))
    fnames.append(str(int(''.join(filter(str.isdigit, file)))))
fnames.pop(0)
for f in fnames:
    if len(f) > 10:
            f = f[10:]
    rlayer = QgsProject.instance().mapLayersByName(str(f))[0]
    if not rlayer.isValid():
        print("ERROR: Raster layer not found!")

##run rastersampling algorithm:------------------------------------
    out=outfolder+str(f)+'_vert.shp'
    processing.run("qgis:rastersampling", {'INPUT':pfile,'RASTERCOPY':rlayer,'COLUMN_PREFIX':var,'OUTPUT':out})
    reslayer=QgsVectorLayer(out,str(f)+"_vert","ogr")
    if not reslayer.isValid():
        print("ERROR: Failed to load result!")
    else:
        QgsProject.instance().addMapLayer(reslayer)

##set crs for results----------------------------------------------
    crs = QgsCoordinateReferenceSystem()
    crs.createFromProj("+proj=ob_tran +o_proj=longlat +o_lon_p=0 +o_lat_p=43 +lon_0=10")
    if not crs.isValid():
        print("CRS is not valid!")
    else:
        reslayer.setCrs(crs)
    QgsProject.instance().addMapLayer(reslayer)
    print(f)