diff --git a/sample/geslotenheid/create_pnts.py b/sample/geslotenheid/create_pnts.py index 97328e62bd1971f59538a259121833a7a621784b..e49c543b3ee15c7dbb5f80f39b7ee260915ab493 100644 --- a/sample/geslotenheid/create_pnts.py +++ b/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')