Hot answers tagged

41

if you have python-gdal bindings: import numpy as np from osgeo import gdal ds = gdal.Open("mypic.tif") myarray = np.array(ds.GetRasterBand(1).ReadAsArray()) And you're done: myarray.shape (2610,4583) myarray.size 11961630 myarray array([[ nan, nan, nan, ..., 0.38068664, 0.37952521, 0.14506227], [ nan, nan, ...


31

Below is an example that I wrote for a workshop that utilizes the numpy and gdal Python modules. It reads data from one .tif file into a numpy array, does a reclass of the values in the array and then writes it back out to a .tif. From your explanation, it sounds like you might have succeeded in writing out a valid file, but you just need to symbolize it ...


22

I've finally hit upon this solution, which I gained from this discussion (http://osgeo-org.1560.n6.nabble.com/gdal-dev-numpy-array-to-raster-td4354924.html). I like it because I can go straight from a numpy array to a tif raster file. I'd be very grateful for comments that could improve on the solution. I'll post it here in case anyone else searches for a ...


21

One possible solution to your problem: Convert it into a ASCII Raster, documention for which is here. This should be fairly easy to do with python. So with your example data above, you'd end up with the following in a .asc file: ncols 4 nrows 4 xllcorner 20 yllcorner 8.5 cellsize 0.5 nodata_value -9999 0.1 0.2 0.3 0.4 0.2 0.3 0.4 0.5 0.3 0.4 0.5 0.6 0.4 0....


20

You can use rasterio to interface with NumPy arrays. To read a raster to an array: import rasterio with rasterio.open('/path/to/raster.tif', 'r') as ds: arr = ds.read() # read all raster values print(arr.shape) # this is a 3D numpy array, with dimensions [band, row, col] This will read everything into a 3D numpy array arr, with dimensions [band, ...


17

Numpy is made for processing arrays and not for reading image files. You need other modules to read the raster and convert it to an array. If you do not want use GDAL or ArcPy: Numpy use Scipy for that: Image manipulation and processing using Numpy and Scipy from scipy import misc raster = misc.imread('image.tif') type(raster) <type 'numpy.ndarray'> ...


16

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 all is read outsize the bounding box, and if image has overviews the data come from the resolution level that is best suitable for the map resolution. There is ...


15

Check out the Describe method. Something like the following should work. #Using arcpy.env import arcpy import numpy inRaster='C:/workspace/test1.tif' outRaster='C:/workspace/test2.tif' dsc=arcpy.Describe(inRaster) arcpy.env.extent=dsc.Extent arcpy.env.outputCoordinateSystem=dsc.SpatialReference arcpy.env.cellSize=dsc.meanCellWidth myArray = arcpy....


15

Matplotlib knows nothing about georeferenced surfaces, it only knows x,y,z coordinates. You can also use Visvis or Mayavi. the original DEM you must first extract the x,y, z coordinates of the grid ( a raster is a grid of pixels and with a DEM, the value of the pixel is the elevation, z) with osgeo.gdal. No script here because it is possible to find the ...


15

I get the same results as gdalwarp from gdal.AutoCreateWarpedVRT if I set the error threshold to 0.125 to match the default (-et) in gdalwarp. Alternatively, you could set -et 0.0 in your call to gdalwarp to match the default in gdal.AutoCreateWarpedVRT. Example Create a reference to compare to: gdalwarp -t_srs EPSG:4326 byte.tif warp_ref.tif Run the ...


15

Instead of doing the reclassification as a double for loop described by dmh126, do it using np.where: # reclassification lista[np.where( lista < 200 )] = 1 lista[np.where((200 < lista) & (lista < 400)) ] = 2 lista[np.where((400 < lista) & (lista < 600)) ] = 3 lista[np.where((600 < lista) & (lista < 800)) ] = 4 lista[np....


14

I was able to get pandas working in ArcMap 10.1 using pandas packages already built for Anaconda. One of the nice things about Anaconda (besides being free) is that it's easy to create environments with specific versions of python and numpy. Try this: Open a 32-bit Anaconda command prompt and type: conda create -n esri101 python=2.7 numpy=1.6 ...


14

rt_raster_to_gdal: Could not load the output GDAL driver As for the first error with ST_AsTIFF, you need to enable your GDAL drivers, which by default are not enabled for PostGIS 2.1. See the manual on ways to do this. For instance, I have an environment variable set up on a Windows computer with: POSTGIS_GDAL_ENABLED_DRIVERS=GTiff PNG JPEG GIF XYZ DTED ...


13

Your PointsXYZIC is now a numpy array. Which means you can use numpy indexing to filter the data you're interested in. For example you can use an index of booleans to determine which points to grab. #the values we're classifying against unclassified = 1 ground = 2 #create an array of booleans filter_array = np.any( [ PointsXYZIC[:, 4] == ...


13

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 lines to set projection and geotransform as input import os import gdal file = "path+filename" ds = gdal.Open(file) band = ds.GetRasterBand(1) arr = band....


12

This should get you going. The raster values are read using rasterio, and pixel centre coordinates are converted to Eastings/Northings using affine, which are then converted to Latitude/Longitude using pyproj. Most arrays have the same shape as the input raster. import rasterio import numpy as np from affine import Affine from pyproj import Proj, transform ...


12

The Python API method that supports the rio-sample command is documented here: https://rasterio.readthedocs.io/en/latest/api/rasterio._io.html#rasterio._io.DatasetReaderBase.sample src.sample() takes an iterator over x, y tuples, so do: for val in src.sample([(x, y)]): print(val)


11

newarray = array * 2.0 performs the math on the entire array, not just on one element. It should instead be something like this: raster = arcpy.Raster(r"C:\test.jpg") array = arcpy.RasterToNumPyArray(raster) # modify cell array[0,0] *= 2.0 # save to a new raster newraster = NumPyArrayToRaster(array) newraster.save(r"C:\export.gdb\t") Or, if you want to ...


10

In the source for gdal_calc.py, the calculation is made directly using eval: myResult = eval(opts.calc, global_namespace, local_namespace) That would suggest that any well-formed expression that also evaluates on the command line will work. According to the documentation, you may use gdalnumeric syntax with +-/*, and/or numpy functions. You can test your ...


8

GDAL 1.10 added a few resampling methods which will help, see gdalwarp. In particular, the -r average method, documented as: average resampling, computes the average of all non-NODATA contributing pixels. This isn't tested, but should look something like: gdalwarp -t_srs EPSG:4326 -tr 0.5 0.66 -r average fine_one_sq_km.tif coarse_average.tif Then to ...


7

There is an option in GDAL to rasterize polygons based on their attribute. But as far as I know it can not be string. But you can just add an attribute to your features and then give each feature a unique id. Let's say we call this field ID. Open your shapefile source_ds = ogr.Open("Longhurst_world_v4_2010.shp") source_layer = source_ds.GetLayer() Create ...


7

Yes, some of the tools use matplotlib. For example (in my 10.1 install): Multi-Distance Spatial Cluster Analysis (Ripleys K Function) <ArcGIS install folder>\ArcToolbox\Scripts\KFunction.py Incremental Spatial Autocorrelation (Moran's I) <ArcGIS install folder>\ArcToolbox\Scripts\MoransI_Increment.py Ordinary Least Squares <ArcGIS ...


7

I have purposed @radouxju comment referencing THIS link as the answer to this question for future viewers: "One solution to this is to explicitly append the PYTHONPATH environment variable to reference the ArcGIS10.1 Python install’s site-packages directory."


7

would add as comment, but a bit long - in case you wanted to use gdal/ogr within python - something like this might work (hacked together from some other code i had - not tested!) This also assumes that rather than finding the nearest raster pixel to a polygon centroid, you simply query the raster at the xy of the centroid. i have no idea what the speed ...


7

The correct name for this function is TableToNumPyArray (arcpy.da) i.e. there should be a capital "T" on "to". Correct capitalization is very important to the Python programming language.


7

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 sets all values between 177 and 185 (inclusive) to 0, which is then treated as nodata. gdal_calc.py -A input.tif --outfile=output.tif --calc="A*logical_or(A<=...


6

Thanks for the help, Branco and om_henners. The answer to my problem appears to be to use numpy.ravel() to change the array produced by arcpy.RasterToNumPy() to a 1D array: import arcpy, pysal from pysal.esda.mapclassify import Natural_Breaks as nb # code to create greenIndex arcpy Raster object here greenArray = arcpy.RasterToNumPyArray(greenIndex) breaks ...


6

I think the question was whether you can read from postgis raster tables WITHOUT gdal drivers enabled. As all things Python, you can! Make sure you select your raster result as WKBinary: select St_AsBinary(rast)... Use the script below to decypher WKBinary into a python image format. I prefer opencv, because it handles arbitrary number of image bands, ...


6

Here you have a simple python script for reclassification, I wrote it and it works for me: from osgeo import gdal driver = gdal.GetDriverByName('GTiff') file = gdal.Open('/home/user/workspace/raster.tif') band = file.GetRasterBand(1) lista = band.ReadAsArray() # reclassification for j in range(file.RasterXSize): for i in range(file.RasterYSize): ...


6

Here's a basic example using rasterio and numpy: import rasterio as rio import numpy as np with rio.open('~/rasterio/tests/data/rgb1.tif') as src: # Read the raster into a (rows, cols, depth) array, # dstack this into a (depth, rows, cols) array, # the sum along the last axis (~= grayscale) grey = np.mean(np.dstack(src.read()), axis=2) ...


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