clipping to ply in viewpoint script

sample/geslotenheid/create_pnts.py

@@ -15,21 +15,25 @@ import rasterio as rio
 mask = rio.open(r'w:\PROJECTS\Landschapsmonitor\cIndicatoren\zOverigeGeodata\provincies\raster\provincies_2018.tif')
 mask_arr = np.where(mask.read(1) > 0, True, False)
 
+
+mask = gp.read_file(r'c:\Users\roelo008\OneDrive - Wageningen University & Research\b_geodata\nl_grenzen\provincies_2018\singly_ply\nl_2018_prov_outline.shp')
+
 # Set extent of output shapefile in meters in RD-New
 xmin = 10000
 ymin = 300000
-xmax = 280000
+xmax = 280000   
 ymax = 625000
 
 # Spacing between points
-spacing = 100
+spacing = 500
 
 # How many points?
 ncols = int(((xmax - xmin)/spacing))
 nrows = int(((ymax - ymin)/spacing))
 
 # Sequence of x-coordinates
-x_coords = [x for x in range(int(xmin + (spacing/2)), int(xmax - (spacing/2)+spacing), spacing)]
+x_coords = [x for x in range(int(xmin + (spacing/2)),
+                             int(xmax - (spacing/2)+spacing), spacing)]
 
 # Sequence of y-coordinates
 y_coords = [y for y in range(int(ymin + (spacing/2)), int(ymax - (spacing/2)+spacing), spacing)]
@@ -38,11 +42,9 @@ y_coords = [y for y in range(int(ymin + (spacing/2)), int(ymax - (spacing/2)+spa
 xx, yy = np.meshgrid(x_coords, y_coords)
 
 # Zip x and y coordinates together with boolean inside/outside mask. Create shapely Point when mask is True
-points = [geometry.Point(x, y) for x, y, m in zip(np.flipud(xx).flatten(),
-                                                  np.flipud(yy).flatten(),
-                                                  mask_arr.flatten())
-          if m]
+points = gp.GeoDataFrame(data={'id': range(1, np.multiply(nrows, ncols)+1)},
+                         geometry=[geometry.Point(x, y) for x, y in zip(np.flipud(xx).flatten(),
+                                                                        np.flipud(yy).flatten())]).set_crs("EPSG:28992")
 
 # Write to file
-gp.GeoDataFrame(data={'id': [i for i,_ in enumerate(points)]},
-                geometry=points).set_crs("EPSG:28992").to_file(r'w:\PROJECTS\Landschapsmonitor\cIndicatoren\Openheid\c_viewscape_input\viewpoints\pnts_100m_v2021_fine.shp')
+gp.clip(points, mask).to_file(r'c:\apps\temp_geodata\viewpoints\viewpoints_500m.shp')