38
votes
Accepted
numpy array to GTiff using rasterio without source raster
This ended up being more straightforward than I thought, with all of the capabilities lying in the rasterio.open function.
Here is an example using a proj4 string instead of wkt.
import rasterio
...
- 1,342
28
votes
Accepted
Reading, modifying and writing a geotiff with GDAL in python
Your script is missing the ds.FlushCache method, that saves to disk what you have in memory at the end of the modifications. See below a corrected version of your example. Notice that I also added two ...
- 745
18
votes
Accepted
How does QGIS open so large raster datasets (about 40GB)?
If QGIS is runnig in a 1000x1000 pixel sized window on your screen there is no need to read all 32000x32000 pixels for showing the map. GDAL tries to read data from the source image so that no data at ...
- 69k
18
votes
Accepted
Creating an in memory rasterio Dataset from numpy array
Note: You can use rasterio.features.geometry_mask to mask your numpy array without writing a dataset (example).
Otherwise if you want to use rasterio.mask.mask, you can create a DatasetReader manually ...
- 69.5k
14
votes
gdal_calc raster calculator syntax for logical operators and other functions
Following on from Benjamin's answer, you can use logical_or() or logical_and(). See http://docs.scipy.org/doc/numpy/reference/routines.logic.html.
The following example worked nicely for me. This ...
- 191
14
votes
Accepted
How to split multiband image into image tiles using Rasterio?
Below is a simple example (rasterio 1.0.0 or later, won't work in 0.3.6). There might be better/simpler ways (and there is an easier way if your raster is internally tiled and the tile block sizes ...
- 69.5k
14
votes
Creating an in memory rasterio Dataset from numpy array
I wrote a wrapper for rasterio.mask.mask that accepts numpy arrays as inputs.
def mask_raster_with_geometry(raster, transform, shapes, **kwargs):
"""Wrapper for rasterio.mask.mask ...
- 241
13
votes
Shapely deprecation warning message when plotting GeoPandas geodataframe
( Answer for users who ended up here based on the title and just want to hide the error )
If you have acknowledged the error and do not want it to keep appearing, you can always ignore it using :
...
- 341
10
votes
How to use gdal_calc.py for multi-band images?
To calculate a grey-scale from the same input file using different bands u can open the file multiple times and define the band which you want to use with --A_band=n.
See my example for calculating ...
- 727
9
votes
Writing numpy array to raster file
An alternative to the approach suggested in the other answers is to use the rasterio package. I had issues generating these using gdal and found this site to be useful.
Assuming you have another tif ...
- 211
9
votes
Writing numpy array to raster file
There is also a nice solution in the official GDAL/OGR Cookbook for Python.
This recipe creates a raster from an array
import gdal, ogr, os, osr
import numpy as np
def array2raster(newRasterfn,...
- 301
8
votes
Accepted
Numpy is not setting properly nan values: arr[arr== 0] = np.nan
Your arrays are not all nan.
It prints nan because that is how np.mean and np.std work, if the array contains any nans, the result will be nan.
You can use nanmean and nanstd instead:
import numpy as ...
- 77.7k
7
votes
Accepted
Looping through all raster cell values using GDAL via Python?
You may read it as array, using numpy:
from osgeo import gdal
import sys
import numpy as np
src_ds = gdal.Open( "INPUT.tif" )
print "[ RASTER BAND COUNT ]: ", src_ds.RasterCount
for band in range( ...
- 16.3k
7
votes
Accepted
create a maximum raster using gdal_calc
you must use maximum instead of max
--calc="maximum(A,B)"
- 49.9k
7
votes
Accepted
Clip raster with another raster (by extent) in Python
Found the solution
#IN ORDER TO CLIP BY EXTENT EVERY IMAGE
gt1=IMG1.GetGeoTransform()
gt2=IMG2.GetGeoTransform()
if gt1[0] < gt2[0]: #CONDITIONAL TO SELECT THE CORRECT ORIGIN
...
7
votes
Accepted
How to create a TIFF file using GDAL from a numpy array and specifying NoData value
The two functions from the code snippet below, create_raster and numpy_array_to_raster should do the trick. In terms of maintaining the NoData value from the array in the output raster, that is set on ...
- 3,266
6
votes
Accepted
NumPyArrayToTable - RuntimeError: unsupported time units. use M8[us]
Playing around with dates using numpy, pandas, and arcpy tools for writing numpy arrays into geodatabase tables can be challenging.
I usually use arcpy.da.InsertCursor in favor of arcpy.da....
- 30k
6
votes
Buffering around raster using gdal and numpy?
I was able to build my own algorithm for this and it's working like a charm.
Since I didn't find this anywhere when I googled it, I'm posting my code here in case someone needs it.
from osgeo ...
6
votes
Accepted
converting a list of shapely geometry to numpy array
I created a shapefile with your point coordinates and the following code produces those numpy arrays:
import fiona
from shapely.geometry import shape
import numpy as np
path = '/home/zeito/...
- 30.2k
6
votes
Accepted
How to filter no data value with GDAL?
If you want to filter no data and get raw values you need following code:
import numpy as np
from osgeo import gdal, gdal_array
dataset = gdal.Open("path/to/file.tif")
array = dataset.ReadAsArray()
...
- 30.2k
6
votes
NumPy array to Raster file (GeoTIFF)
By using a raster with integer values (1, 100) and one equivalent condition (myarray >= 35, myarray <= 7), following code would work as expected:
from osgeo import gdal, osr
import numpy
ds = ...
- 30.2k
6
votes
Accepted
How to change a numpy array's dimensions
numpy.transpose is one way of doing this.
import numpy as np
zyx = np.ones((1, 2, 3)) # 1 band, 2 rows, 3 cols
yxz = np.transpose(zyx, (1,2,0))
print(yxz.shape)
# (2, 3, 1)
- 69.5k
6
votes
Accepted
Create pandas DataFrame from raster image - one row per pixel with bands as columns
Quick solution
pd.DataFrame(array.reshape([3,-1]).T)
Explanation
Take array of shape (3, x, y) and flatten out the 2nd and 3rd dimension. From the numpy docs: One shape dimension can be -1. In this ...
- 3,771
6
votes
Accepted
Issue in calculating NDVI using Rasterio Python
In order to solve your problem, you need to ensure the grids cover the same area and have the same dimensions. One method to achieve this is with the reproject_match method in rioxarray (geospatial ...
- 7,527
6
votes
Strange warning/error when working with polygons
This is not an error on user side. It is the result of a recent change in Numpy 1.21 coupled with a way shapely implements (rather not implements) __array_interface__. There is nothing the user should ...
- 2,447
6
votes
Accepted
Adding CSV file by Python to QGIS throws #Line20 got ...columns instead of... error
The problem is caused by the comma in lat/long values. In your previous question, lat/long contains dot instead of comma for decimal (this is the main reason why you get the error) and lat/long ...
- 78k
5
votes
Fully load raster into a numpy array?
My solution using gdal looks like this. I think it is very reusable.
import gdal
import osgeo.gdalnumeric as gdn
def img_to_array(input_file, dim_ordering="channels_last", dtype='float32'):
file ...
- 81
5
votes
Reconstructing MODIS time-series applying Savitzky-Golay Filter with Python/Numpy
Based on the SG filter from scipy.signal I built the NDVI timeseries smoothing algorithm proposed in:
A simple method for reconstructing a high quality NDVI time-series data set based on the ...
- 265
5
votes
Accepted
Inserting LiDAR points (from laspy) in GeoDataFrame without using a numpy array?
So, is there a better way, more efficient, more pythonesque way of injecting the .las file points (from laspy) to the GeoPandas dataframe without passing through a numpy array?
No, and I'm not sure ...
- 4,106
5
votes
Accepted
Classifying LiDAR ground points using laspy?
I was able to come up with a solution. The following code works, but it feels like there is a much simpler way to do this operation. If anyone has a cleaner or more efficiect way to filter and ...
- 371
Only top scored, non community-wiki answers of a minimum length are eligible