biglimp / PhDCourseVT2021

This repo include teaching materials for the PhD-course "Automate your GIS - Scripting in Python" (NGEO306)
7 stars 1 forks source link

Help with nested loop #4

Open Damboia opened 3 years ago

Damboia commented 3 years ago

Hi..

Finally I could managed to calculate the area of my polygons ( I had to fix geometries but finally it runs kakaka) for one of the months of my data but now I would like to run the script for other months and plot the graphs with differences. But I am bit confusing of how can I do it. The data as same format (.shp files).

Loading a vector layer

vlayer = QgsVectorLayer ('c:/Cossa/Thesis_Data/FixedJulyScars_2020.shp', 'FixedJulyScars_2020','ogr') print(vlayer.isValid()) vlayer2 = QgsVectorLayer ('c:/Cossa/Thesis_Data/DugongScars_Aug20.shp', 'DugongScars_Aug20','ogr') print(vlayer2.isValid())

check attributes in vector layers

for features in vlayer.getFeatures(): attrs = features.attributes() print(attrs[0]) for features in vlayer2.getFeatures(): attrs = features.attributes() print(attrs[0])

geom = features.geometry()
if geom.isMultipart() :
    print (vlayer.name(), ': Warning: Multipart', features. id())
    print (vlayer2.name(), ': Warning: Multipart', features. id())

vlayer = iface.activeLayer()

caps = vlayer.dataProvider().capabilities() fields_name = [f.name() for f in vlayer.fields()]

we check if we can add an attribute to the layer

if caps & QgsVectorDataProvider.AddAttributes: if 'Area' not in fields_name: vlayer.dataProvider().addAttributes ( [ QgsField ('Area',QVariant.Double) ] ) vlayer.updateFields() fields_name = [f. name() for f in vlayer.fields()] fareaidx = fields_name.index ('Area')

calculate area of polygons

if caps & QgsVectorDataProvider.ChangeAttributeValues: for gFeat in vlayer.getFeatures(): d = QgsDistanceArea() d.setSourceCrs(vlayer.crs(), QgsProject.instance().transformContext()) d.setEllipsoid (QgsProject.instance().ellipsoid())

for features in vlayer.getFeatures():
geom = features.geometry() area = [0] landArea = gFeat ['Area'] if geom.isMultipart(): polyg = geom.asMultiPolygon() for polyg in polyg: area = d.measurePolygon(polyg[0]) landArea = area print(landArea)

biglimp commented 3 years ago

Start by listing the monthly shapefiles using e.g. import pathlib monthlyFiles = list(pathlib.Path('your_directory').glob('*.shp'))

then you can make a loop outside your code i.e. for i in monthlyFiles: load monthly vector layer calulate area

To load a monthly file in the loop, replace your shapefile name (e.g. c:/Cossa/Thesis_Data/FixedJulyScars_2020.shp) with i in your QgsVectorLayer-function

Damboia commented 3 years ago

Hi ..

Something like that: monthlyFiles = list(pathlib.Path('c:/Cossa/Thesis_Datafinal').glob('*.shp')) print (monthlyFiles)

for i in monthlyFiles: vLayer= QgsVectorLayer(monthlyFiles, 'i', 'org') list (vLayer)

biglimp commented 3 years ago

yes. @nilswallenberg , Can you help Damboia if she need more assistance today. I am teaching today.

nilswallenberg commented 3 years ago

Yes, I can help!

Damboia commented 3 years ago

Hi Nils, thank you..

Well, I am a bit stuck now because I am tying to load a monthly file to loop but it says

File "c:/Cossa/Projectscriptmonths_Damboia.py", line 66, in vlayer = QgsVectorLayer(monthlyFiles,'i','org') TypeError: QgsVectorLayer(): argument 1 has unexpected type 'list'

nilswallenberg commented 3 years ago

Hi Damboia,

If this is related to the code in your comment from two days ago I suggest you try this:

monthlyFiles = list(pathlib.Path('c:/Cossa/Thesis_Datafinal').glob('*.shp'))
print (monthlyFiles)
for i in monthlyFiles:
    vLayer= QgsVectorLayer(i, 'ogr')
    list (vLayer)

you can also add this in your for loop to see what i actually is: print(i)

Damboia commented 3 years ago

Hi Nils thanks,

It list the monthly files and recognize the windows path. Seems that the problem is when I tried to create Vectorlayer for months (load monthly data). I don't know if there is something related with my database.

I tried your suggestion and now says that : "TypeError: 'QgsVectorLayer' object is not iterable"

nilswallenberg commented 3 years ago

Can you show which row that you get this error? Maybe it is list(vLayer) that is causing this error.

Damboia commented 3 years ago

Hi..Nils, Thank you and this error its solved, but now says that my vLayer is false so, seems that I cannot run my script :(

Thesis_Datafinal/._DugongScars_Aug20.shp'), WindowsPath('c:/Cossa/Thesis_Datafinal/DugongScarsfinal_June20.shp'), WindowsPath('c:/Cossa/Thesis_Datafinal/DugongScars_Aug20.shp'), WindowsPath('c:/Cossa/Thesis_Datafinal/DugongScars_Dec20.shp'), WindowsPath('c:/Cossa/Thesis_Datafinal/DugongScars_July20.shp'), WindowsPath('c:/Cossa/Thesis_Datafinal/DugongScars_Oct20.shp'), WindowsPath('c:/Cossa/Thesis_Datafinal/DugongScars_Sept20.shp')] c:\Cossa\Thesis_Datafinal._DugongScarsfinal_June20.shp c:\Cossa\Thesis_Datafinal._DugongScars_Aug20.shp c:\Cossa\Thesis_Datafinal\DugongScarsfinal_June20.shp c:\Cossa\Thesis_Datafinal\DugongScars_Aug20.shp c:\Cossa\Thesis_Datafinal\DugongScars_Dec20.shp c:\Cossa\Thesis_Datafinal\DugongScars_July20.shp c:\Cossa\Thesis_Datafinal\DugongScars_Oct20.shp c:\Cossa\Thesis_Datafinal\DugongScars_Sept20.shp c:\Cossa\Thesis_Datafinal._DugongScarsfinal_June20.shp False c:\Cossa\Thesis_Datafinal._DugongScars_Aug20.shp False c:\Cossa\Thesis_Datafinal\DugongScarsfinal_June20.shp False c:\Cossa\Thesis_Datafinal\DugongScars_Aug20.shp False

nilswallenberg commented 3 years ago

You will be able to run your script. We just have to figure out why it says False. :)

Try this snippet to create a list with files:

file_directory = 'C:/Cossa/Thesis_Datafinal/'
monthlyFiles = [f for f in glob.glob(file_directory + '/*.shp')]

for file in file_directory:
    vlayer = QgsVectorLayer(file, 'ogr')
    print(vlayer.isValid())

Maybe there is something with the paths in your previous code.

Damboia commented 3 years ago

Thank you Nils rsrsrs.

Says again that is False.

I import glob to script, so I think its not the problem.

nilswallenberg commented 3 years ago

Okay, we can have a look on Zoom today if you want to.

Damboia commented 3 years ago

Hi.. I want, yes...Let me know when you are available :) and the link

Damboia commented 3 years ago

Hi Nils,

Seems that we still with an issue on vLayer. I did a double check and I found the area was calculated for just one of the months. I started debugging and and look what says when I run in vLayer row

At line:1 char:39

nilswallenberg commented 3 years ago

Okay, come to Zoom again please (same as before). :)