Tag Info

Hot answers tagged

7

You could use gdal_rasterize either from the command line or QGIS to generate your raster. To make sure your points sit within a cell, you need to do two things. First, set the target resolution to 5m, and set the extents to be 2.5m bigger all around than the source data. So, assuming your dataset goes from [1000 2000] [2500 3250], giving you your 75,000 ...


6

OK so a second attempt to answer your question with a pure GDAL solution. Firstly, GDAL (Geospatial Data Abstraction Library) was originally just a library for working with raster geo-spatial data, while the separate OGR library was intended to work with vector data. However, the two libraries are now partially merged, and are generally downloaded and ...


5

If your data contains xyz data (where z is the raster value) and your points are on a regular grid (no need for interpolation). library("raster") r <- rasterFromXYZ(as.data.frame(travel)[, c("x", "y", "z")]) If you need interpolation, you can use akima library : library("raster") library("akima") steps <- 100 isu <- with(travel@data, interp(x, ...


5

I've had some luck reading from and writing to layers. Specifically, I have code that will read a shapefile layer containing polylines and output the geometry of each feature to text files (used as input for an old model). name = layer.name() provider = layer.dataProvider() feat = QgsFeature() # Now we can loop through all the defined features ...


5

I've been playing with GDALRasterizeLayers this week and have a pretty good idea of what it is doing. By default, it will rasterise a pixel if the pixel centre is within the polygon. If there is nothing in the centre, it won't be rasterised, even if there are parts of a polygon within the pixel limits. To allow the rasterising to work the way you intend, try ...


5

Following a recent question, you may want to make use of the functionalities offered by the rgeos package to solve your problem. For reasons of reproducibility, I downloaded a shapefile of Tanzanian roads from DIVA-GIS and put it in my current working directory. For the upcoming tasks, you will need three packages: rgdal for general spatial data handling ...


4

In addition to @Etiennebr's answer, I'd go for an apply style loop (which is more R-ish, and uses less code for the same thing): library("raster") filenames <- list.files(path="", pattern="XYhectareTravelTimes_ez+.*shp") raster_cat = lapply(filenames, function(x) { travel <- as.data.frame(readShapePoints(x)) r <- rasterFromXYZ(travel[, c("x", ...


4

You could try to create a Voronoi diagram. It seems to work: If the Voronoi tool should fail on 75,000 points, try divide&conquer: Split the point layer into multiple smaller ones to make computation easier.


3

This is a kludge but it does work - haven't found a way to go directly from points to raster yet (but am hoping someone gives a solution here!). Starting with a point grid (random points in the Serengeti from the Vector|Research tools|Random points tool): Create a polygonal grid of the same extent and cell size as the raster you'd like to have (this from ...


3

If you want to get zonal statistics for several features in one shapefile, you have to loop over the zonal_stats function. You can write the results of the loop for example to a dictionary. Below is the modified zonal_stats function together with a loop, looping over the input shapefile. As an output you get a Dictionary containing for each Feature ID the ...


3

Two things you could check: First make sure that your resolution is set correctly. This you do with the "Edit Current Grass Region" button. You mentioned that the original csv file had points every 40 km. If you want to stay with that, then assuming your points are in a projected coordinate system you should set the resolution to 40,000. If the points are ...


3

I've not tried it, but GDAL's gdal_rasterize should do the trick with its -3d option: gdal_rasterize -3d -tr 10.0 10.0 -l streams streams.shp streams.tif


3

According to the GDAL Proximity page: -distunits PIXEL/GEO: Indicate whether distances generated should be in pixel or georeferenced coordinates (default PIXEL). So you can specify the distances in pixels or in the unit of the coordinates you used to georeference the raster.


2

Go Raster -> Conversion -> Rasterize. Set a vector layer to process, a field with values and desired raster size.


2

This particular function does require a regularly spaced set of points. The error message is reporting that there are gaps in your data and it's not able to interpolate between the existing points to provide the value for the missing points. You need to see if R has a different rasterizing function that will interpolate to fill in the gaps, or find some ...


2

Your error says as.integer(putvals). The R rasterize function can't work on Strings. You have to transform your data first. Something like this may work, but i would assign different ranks (ala 1,2,3,4,..) to your data. However i still get a (different) error for which i don't have an explanation. Maybe the size of your raster is incorrect... ...


2

gdal_rasterize is not a python module, but an exe file. You can not call it that easy inside the python console. See this blogpost on how to wrap the operating system around it: http://geoinformaticstutorial.blogspot.de/2012/11/convert-shapefile-to-raster-with-gdal.html import os os.system('gdal_rasterize -a ICE_TYPE -where \"ICE_TYPE=\'Open Water\'\" ...


2

You could add a field to the lake shapefile called "value", and populate it with the value 1. In your Rasterize tool, set the value field to "value", and create your raster. Then just clip the raster using your min/max x & y values using the Clipper tool? Raster -> Extraction -> Clipper The output will be a clipped raster with 1 values for ...


2

This can be done in a two step process... Use Merge to combine all your polygons into one polygon (tool found in the Data Management / General toolbox) Use Feature to Raster to change the polygon into a raster (tool found in the Conversion / To Raster toolbox). These are both straightforward tools. The most complicated part with Feature to Raster is ...


1

You do not need a for loop. Just intersect everything at once and then add line lengths to the new line segments using the "SpatialLinesLengths" function in sp. Then, using the raster package rasterize function with the fun=sum argument you can create a raster with the sum of the line length(s) intersecting each cell. Using the above answer and associated ...


1

I was faced with the exact same problem some days ago. The only solution is to make a copy of the original grid and replace it. in that case you can use -TARGET 1 -GRID_GRID result.sgrd The result.sgrd file will be overwritten with your result, but will keep the boundaries that were present in the file you used. I do agree that a better solution should ...


1

With multiprocessing, for fastness! Has a little different output-formatting. #!/usr/bin/python import gdal, ogr, osr, numpy, sys from multiprocessing import Pool # Raster dataset input_value_raster = sys.argv[1] # Vector dataset(zones) input_zone_polygon = sys.argv[2] # Open data rast = gdal.Open(input_value_raster) shp = ogr.Open(input_zone_polygon) ...


1

Easiest and most straight-forward way: Assume you have a column with a unique point identifier (the species name) Split your Point layer by this Attribute (QGIS -> Data Management -> Split) Rasterize each individual point layer for instance with the GDAL Rasterize Tool, or the SAGA or GRASS tools available in the Processing Toolbox. Make sure that you use ...


1

I assume that you are using the GDAL Rasterize (vector to raster) tool. Using QGIS, you can edit the gdal rasterize command and had the "xmin ymin xmax ymax" parameter: Notice that you will have to adapt the coordinates to fit your case. You can use this same command directly on the console or terminal to bash rasterize all your vector in one go. You ...


1

I have noticed the same error when clipping to Pseudo mercator. It might have something to do with the special treatment of mercator projection on a sphere. Try any other CRS when clipping to avoid the bug. You can reproject afterwards back to EPSG:3857 if necessary. It might help to set -s_srs and -t_srs explicitely to EPSG:3857 in both command lines, but ...


1

Actually I found an alternative solution to my problem using R. Here is my sample code. In the sample code, xmn, xmx, ymn, ymx are the coordinates of the pre-defined extent. "LAKE_ID" is the attribute with value equal to 1 for all polygons. library(rgdal) library(raster) library(sp) ymn <- 6758000 ymx <- 6766000 xmn <- 422000 xmx <- 429000 ...


1

Assuming you know how to create the hectare sized grid (using the Vector->Research tools->Vector grid) then you can import both the point shapefile and the grid polygons into SpatiaLite and with an aggregate query get the totals as an attribute column in each grid sqare. For example if you have an id column in the grid polygon layer, and a ...


1

Two methods come to mind: You can use the gdal utility, gdal_rasterize, with the -a option to specify the attributie column, and the -te and -tr options to set the output extents and resolution the same for each raster. See details on the GDAL website. Alternatively, if you work in GRASS, then you first set the desired extents ("region settings" in GRASS ...



Only top voted, non community-wiki answers of a minimum length are eligible