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 ...
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 ...
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:
0.1 0.2 0.3 0.4
0.2 0.3 0.4 0.5
0.3 0.4 0.5 0.6
You can use rasterio to interface with NumPy arrays. To read a raster to an array:
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, ...
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')
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 ...
Check out the Describe method.
Something like the following should work.
myArray = arcpy....
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 ...
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.
Create a reference to compare to:
gdalwarp -t_srs EPSG:4326 byte.tif warp_ref.tif
Run the ...
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.
Open a 32-bit Anaconda command prompt and type:
conda create -n esri101 python=2.7 numpy=1.6 ...
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 ...
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] == ...
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
file = "path+filename"
ds = gdal.Open(file)
band = ds.GetRasterBand(1)
arr = band....
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 numpy as np
from affine import Affine
from pyproj import Proj, transform
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)
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)
Or, if you want to ...
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 ...
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 ...
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()
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
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."
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 ...
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<=...
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)
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:
Use the script below to decypher WKBinary into a python image format. I prefer opencv, because it handles arbitrary number of image bands, ...
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()
for j in range(file.RasterXSize):
for i in range(file.RasterYSize):
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)