I assume you want your irregular point data on a regular raster. In that case, rasterize should work, and the examples in ?rasterize show how. Here is something based on your data
s100 <- matrix(c(267573.9, 2633781, 213.29545, 262224.4, 2633781, 69.78261, 263742.7, 2633781, 51.21951, 259328.4, 2633781, 301.98413, 264109.8, 2633781, 141.72414, 255094.8, ...
here is an example. You need to extract the bounding box from a describe object.
desc = arcpy.Describe(fc)
arcpy.CreateFishnet_management(fc[:-4]+"_c200.shp",str(desc.extent.lowerLeft),str(desc.extent.XMin) + " " + str(desc.extent.YMax + 10),"200","200","0","0",str(desc.extent.upperRight),"NO_LABELS","#","POLYGON")
here is a solution using sf and sf::st_make_grid:
e <- as(raster::extent(-125, -65, 25, 49), "SpatialPolygons") %>%
grd_lrg <- st_make_grid(e, cellsize = c(5, 2)) %>%
use the MMQGIS plugin, after installing you find it in qgis 1.8 under plugins--mmqgis. select create --create grid layer and enter the appropriate values
for h spacing and v spacing select 1 or 0.5 (instead of 10 in my screenshot)
hope this helps
This is quite simple you can use gdaltindex to build this grid.
In QGIS go to Raster -> Misc. -> Tile index. Give the path to the raster-files, set the name and directory for the resulting shape. I prefer absolute paths for catalogue purposes.
Then hit ok. :)
The following script will do the job with GDAL and Python:
import os, sys
from math import ceil
# convert sys.argv to float
xmin = float(xmin)
xmax = float(xmax)
ymin = float(ymin)
ymax = float(ymax)
gridWidth = float(gridWidth)
gridHeight = float(...
The New York dataset provided in the question is no longer available for download. I use the nc dataset from sf package to demonstrate a solution using sf package:
# read nc polygon data and transform to UTM
nc <- st_read(system.file('shape/nc.shp', package = 'sf')) %>%
# random sample of 5 points
You need to use Split Vector Layer from processing toolbox. You can find the tool from Processing toolbox -> QGIS geoalgorithms -> Vector general tools -> Split Vector Layer
Select the input shapefile and use a unique ID field to save each grid block into separate file. You need to specify the output folder to save all the polygons in it.
Vector grid creation has been incorporated in QGIS for a long time:
Vector -> Research Tools -> Vector Grid
Note that the default grid spacing of 0.0001000000 units of the layer CRS is too small in most cases.
You might want to use Vector -> Geometry Tools -> Densify Geometries afterwards if you need to reproject the grid to another CRS where ...
Create a polygon grid using the Vector Grid tool instead of lines. Make sure to check the polygon output.
Once you have a polygon grid (also known as fishnet), you can use the Sum line length tool in the QGIS Vector analysis tools. This will result in a new field for each cell with the total road length inside it
Here's a simple example of a vector ...
You can do this with both the Create Fishnet tool or the Grid Index Features tool.
Create Fishnet is a bit more configurable, and can output either polygon or polyline. Both will allow you to set a width and height.
Instead of looking for the Create Fishnet (or any other) tool in ArcToolbox, I recommend using the Search window to search for and open it.
If you are curious to know which toolset it is located in for the ArcGIS version that you are using then after searching for it you can click on the green link returned.
It looks like you're using a geographic coordinate system instead of a projected coordinate system. Because NAD83 (HARN) uses degrees as the units, when you type in 500 in the grid spacing, it's using a spacing of 500 degrees, which is huge. You'll need to reproject the data to a more suitable projected coordinate system which uses metres such as UTM Zone ...
As Vince said, the Create Fishnet tool is what you want to use. At the top of that help page in the summary section is a link to the How Create Fishnet Works page. Note that the main help page for most tools has a chart outlining various parameters, but many of them have such a 'how it works' link the same place that explains in more detail, often with ...
This way it should work:
Reproject your points to a projected CRS with metres as units (if not already done)
Draw buffers with a radius of 50m (half of the grid size) around them with Vector -> Geoprocessing Tools -> buffer(s)
Create the vector grid with the extent of the buffer layer and 100m spacing
You can use the floor function in your expression, like this:
floor($id /24) IN (1,3,5,7,9,11,13,15,17,19,21)
On a sample grid I've obtained this result with my own numbers (my grid has 13 columns per row):
I've used: floor($id/13) IN (1,3,5,7,9)
So, in general:
floor($id / number_of_columns) IN (1, 3, 5, ..., 2*number_of_rows_to_be_selected - 1)
I just ended up creating my own tools to do this.
I used the Clipper library (http://www.angusj.com/delphi/clipper.php) along with OGR to divide my data set up. Something to note is performing intersections naively with this lib takes very long, so I instead used a quadtree approach... ie, divide into four grid cells, divide each of those into four more, ...
The ogr2ogr clip operation creates polygons but for some reason also linestrings to to South and East edges of the area. Pink lines below show those 100 linestings (one is selected). Because of mixed geometrytypes the result cannot be saved into shapefile and therefore the error. I am not sure if this is an intended behaviour of ogr2ogr clipping or a bug. ...
The first part of your question is easy. Go to Vector->Analysis Tools->Points in Polygon. This tool will iterate over all your polygons and add a column (called 'PNTCNT' by default) with the number of points in each polygon.
This tool gives you just a simple count of the number of points in each polygon as the name suggests. You can refine this by ...
Previous answer gives correct result, but I allowed myself to improve the code, to avoid many unnecessary lines, as well as iterating on indexes instead of values on the list.
import geopandas as gpd
from shapely.geometry import Polygon
import numpy as np
points = gpd.read_file('points.shp')
xmin,ymin,xmax,ymax = points.total_bounds
length = 1000
wide = ...
You can make a long/lat raster and transform
r <- raster(ext = extent(-125, -65, 25, 49), res=c(5,2))
values(r) <- 1:ncell(r)
rA <- projectRaster(r, crs="+proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=23 +lon_0=-96 +x_0=0 +y_0=0 +datum=NAD83 +units=m")
But what you want seems to be
p <- rasterToPolygons(r)
# or do: p ...
That was a fun problem! I solved it using QGIS 2.18 but I don't think the tools have drastically changed in 3.0.
1. Generate the lines
I won't give much details here since you already have the lines you need. I have written a script to generate these lines, but it works only up to QGIS 2.18. Needless to say, the more lines the better your estimation.
You have two options here. First one will create a new vector layer with a grid and the second one will just display a grid as overlay. What you are looking for is the first option. However, in case someone else is finding this question, the second option might be helpful as well.
Edit: of course 10m xy, not 100m as said below...
Run "Create ...
I've created 4 tools which might help you out.
1) The first 2 tools create a quarter section grid from an existing
section grid and provide custom labeling for each quarter section.
2) The other 2 tools create a quarter quarter section grid from an
existing section grid and provide custom labeling for each quarter
I say from ...