I'm now stuck with writing a piece of code to look up individual pixel values on a raster image.   This piece of code mainly executes two tasks: (1) converting world coordinates to raster coordinates (I use a package named "Affine" to do that); (2) then computing pixel offsets and figuring out which column and row the pixel is located at.  The major problem I'm facing now is that the code didn't return the correct column and row that the pixel is located at.   import sys, gdal, numpy from affine import Affine  # coordinates to get pixel values for lons = [-122.265373, -122.429139, -68.123095] lats = [37.873090, 37.783640, 44.981935]  # open the raster file ds = gdal.Open('data/nlcd_compressed.tif') if ds is None:     print 'Could not open the raster file'     sys.exit(1)  # get raster size cols = ds.RasterXSize #161190  rows = ds.RasterYSize #104424 bands = ds.RasterCount #1  # get georeference info transform = ds.GetGeoTransform() # (-2493045.0, 30.0, 0.0, 3310005.0, 0.0, -30.0) xOrigin = transform[0] # -2493045.0 yOrigin = transform[3] # 3310005.0 pixelWidth = transform[1] # 30.0 pixelHeight = transform[5] # -30.0  # convert global coords to raster coords aff = Affine.from_gdal(-2493045.0, 30.0, 0.0, 3310005.0, 0.0, -30.0) x_coords, y_coords = aff*(numpy.array(lons), numpy.array(lats)) # converted coordinates xValues = [-2496712.96119, -2496717.87417, -2495088.69285] yValues = [3308868.8073, 3308871.4908, 3308655.54195]   band = ds.GetRasterBand(1) # 1-based index #bandtype = gdal.GetDataTypeName(band.DataType) #Byte #data = band.ReadAsArray(0, 0, cols, rows) #value = data [122, 37] #127  data = band.ReadAsArray(xOffset)  # loop through the coordinates for i in range(len(lons)):     x = xValues[i]     y = yValues[i]     # compute pixel offset     xOffset = int((x - xOrigin) / pixelWidth)     yOffset = int((y - yOrigin) / pixelHeight)     # get individual pixel values     data = band.ReadAsArray(xOffset, yOffset, 1, 1)     value = data[0, 0]     print value
41081
I am working on a computational model of the abundance of wild pollinators across a landscape.  The model itself is complete, and I am now struggling with a post-processing step.  I have my GDAL pollinator supply raster that looks something like this (lighter colors mean higher pollinator visitation to a pixel):    And I have an OGR shapefile of points representing sample locations on the landscape:    I'm trying to run some analysis on the pixels under these points, but to do so, I need to be able to extract the value of a pixel under a point.  Is it possible to extract the value of a pixel under a point using only OGR and GDAL through Python?  I would prefer to avoid reading the entire raster into memory through ReadAsArray(), as my output rasters are very, very large (too large to fit into memory).  I noticed , which is similar, but requires a command-line call.
44745
I have been attempting to use the Spatial Analyst Curvature tool with mixed results.  I originally tried to use it with a raster that was projected using a 1983 Teale Albers projection and the results were black, indicating the lowest values.  I then transformed the raster into UTM projection and ended up with grid-like lines on the raster which have values that are obviously not accurate.  Does anyone know why these lines have appeared or what I'm doing wrong.  Thanks.
26515
I have written some code that calculates the world space coordinate I want to move the bone to but i am having trouble figuring out how to get that back into the bones coordinates.  from bpy import context, data from mathutils import Vector  def test_empty(matrix, name):     # view matrix by applying it to an object     obj_empty = data.objects.new(name, None)     context.scene.objects.link(obj_empty)     obj_empty.matrix_world = matrix   armature = data.objects["Armature"] pbones = armature.pose.bones root = pbones["root"] hips = pbones["hips"] test = pbones["test"]  # calculate root root_matrix = armature.matrix_world * root.matrix test_empty(root_matrix, "Root")  # calculate hips hips_matrix = armature.matrix_world * hips.matrix test_empty(hips_matrix, "Hips")  # calculate bone test_matrix = armature.matrix_world * test.matrix test_empty(test_matrix, "Test")  # get offset of hips from root bone offset = root_matrix.translation - hips_matrix.translation print(offset)  # calculate new location in local space new_position = test_matrix.translation + Vector((0,offset.y,0))  obj_empty = data.objects.new("TEST", None) context.scene.objects.link(obj_empty) obj_empty.location = new_position  # obviously something wrong here new_position = armature.matrix_world.inverted() * new_position  # move bone to new location test.location = new_position
114747
I have a Python function which converts Eastings and Northings to Latitude and Longitude if line_count == 0:                     keys = row                     print('keys', keys)                     line_count += 1                 else:                     inProj = Proj(init='epsg:27700') # British National Grid                     outProj = Proj(init='epsg:4326') # WGS84                     x1, y1 = row[1], row[2]  # easting, northing                     x2, y2 = transform(inProj, outProj, x1, y1)                     row[1], row[2] = y2, x2                                          json_doc = {}  At the moment it converts from epsg:27700 (British National Grid) to WGS84, however there is a possibility I will get TM75 Irish Grid too. Is there any way I can get Python to recognised which CRS is applicable for the conversion?
97331
I am trying to burn a shapefile to a raster using GDAL's RasterizeLayer. I pre-create an area of interest raster from a different shapefile, given a specific pixel size. This AOI then serves as a base for all following rasterizations (same number of collumns and rows, same projection and geotransform).  The issue occurs, however, when I go to burn the shapes to their own raster, based on the same pixel size and projections. The link below (don't have enough rep to post the image), shows the original shapefile in tan, and the dark pink where RasterizeLayer has burned data. The light pink is nodata values for the dark pink raster data. The grey is the AOI based on which the shapefile burn was completed.  Given the extents of the shapefile polygons, I would expect to see burn values in the bottom two corners, as well as the two pixels underneath the data that does show. Obviously, however, this is not the case.    As follows is the code that I used to generate these. All of the shapes were created using QGIS, and were all created in the same projection.  (It should be noted that the gridding in the picture shown was just to give an idea of the pixel size I was using.)  from osgeo import ogr from osgeo import gdal  aoi_uri = 'AOI_Raster.tif' aoi_raster = gdal.Open(aoi_uri)  def new_raster_from_base(base, outputURI, format, nodata, datatype):      cols = base.RasterXSize     rows = base.RasterYSize     projection = base.GetProjection()     geotransform = base.GetGeoTransform()     bands = base.RasterCount      driver = gdal.GetDriverByName(format)      new_raster = driver.Create(str(outputURI), cols, rows, bands, datatype)     new_raster.SetProjection(projection)     new_raster.SetGeoTransform(geotransform)      for i in range(bands):         new_raster.GetRasterBand(i + 1).SetNoDataValue(nodata)         new_raster.GetRasterBand(i + 1).Fill(nodata)      return new_raster  shape_uri = 'activity_3.shp' shape_datasource = ogr.Open(shape_uri) shape_layer = shape_datasource.GetLayer()  raster_out = 'new_raster.tif'  raster_dataset = new_raster_from_base(aoi_raster, raster_out, 'GTiff',                                 -1, gdal.GDT_Int32) band = raster_dataset.GetRasterBand(1) nodata = band.GetNoDataValue()  band.Fill(nodata)  gdal.RasterizeLayer(raster_dataset, [1], shape_layer, burn_values=[1])   Is this a bug in GDAL, or is RasterizeLayer burning data based on something other than just the presence or lack of a polygon within a specified pixel area?  The files that I was using can be found .
43376
I'm trying to white a Python plugin in QGIS (2.2) to get, for each pixel of a raster image, its X Y coordinates, as well as, its RGB (3 bands) values.  At first I opened the raster file using: rlayer = QgsRasterLayer(rfileName, rbaseName)  Now I don't know how to get, for example, for pixel (1,1) its coordinates (X,Y) and its RGB color values.  Anybody could help me?  I know that I'll need to implement a DO WHILE, but I don't know the commands to extract this information of each pixel.  Thanks in advance, Mateus
113293
I'm trying to write a script that can automatically use GDAL to crop the informational border from the FAA's published sectional charts (GeoTIFF) found :  I know that I can use gdal_translate -projwin ulx uly brx bry to crop it to a coordinate range, but I'm not sure where to get the values for the map area aside from by hand, or if there is possibly a better way to crop it.  Is the information necessary to crop the chart extractable from the files they provide?
81996
Trying to load and plot a raster (it's a DEM) in R, but having problems.    Using R I read in a raster file:  library(raster)  dem &lt;- raster("DEM_Victoria.tif")  # &gt;70 MB print(dem)  class       : RasterLayer  band        : 1  (of  3  bands) dimensions  : 5444, 6007, 32702108  (nrow, ncol, ncell) resolution  : 100, 100  (x, y) extent      : 176615.6, 777315.6, 5650397, 6194797  (xmin, xmax, ymin, ymax) coord. ref. : +proj=longlat +ellps=WGS84 +no_defs  data source : DEM_Victoria.tif names       : DEM_Victoria values      : 0, 255  (min, max)   My confusion here is the extent vs the CRS. The extent shows coordinates from a UTM Zone 55S system, while the CRS reports it as longlat.   Similarly, gdalinfo reports:  Driver: GTiff/GeoTIFF Files:  DEM_Victoria.tif Size is 6007, 5444 Coordinate System is: GEOGCS["GDA94",     DATUM["unknown",         SPHEROID["unretrievable - using WGS84",6378137,298.257223563]],     PRIMEM["Greenwich",0],     UNIT["degree",0.0174532925199433]] Origin = (176615.587999997427687,6194797.138000000268221) Pixel Size = (100.000000000000000,-100.000000000000000) Metadata:   AREA_OR_POINT=Area Image Structure Metadata:   COMPRESSION=LZW   INTERLEAVE=BAND Corner Coordinates: Upper Left  (  176615.588, 6194797.138) (Invalid angle,Invalid angle) Lower Left  (  176615.588, 5650397.138) (Invalid angle,Invalid angle) Upper Right (  777315.588, 6194797.138) (Invalid angle,Invalid angle) Lower Right (  777315.588, 5650397.138) (Invalid angle,Invalid angle) Center      (  476965.588, 5922597.138) (Invalid angle,Invalid angle) Band 1 Block=6007x1 Type=Byte, ColorInterp=Red Band 2 Block=6007x1 Type=Byte, ColorInterp=Green Band 3 Block=6007x1 Type=Byte, ColorInterp=Blue   Is this real? A problem? (well it is to me, but perhaps because I don't understand something) Fixable?  I tried 'Assign projection' in QGIS, which uses this command:  gdalwarp -t_srs EPSG:28355 -srcnodata -32768  DEM_Victoria.tif DEM_Victoria.tif.tmp   but it failed with:     ERROR 1: Too many points (12100 out of 12100) failed to transform,   unable to compute output bounds.
75830
Alright, I'm toying around with converting a PIL image object back and forth to a numpy array so I can do some faster pixel by pixel transformations than PIL's PixelAccess object would allow.  I've figured out how to place the pixel information in a useful 3D numpy array by way of:  pic = Image.open("foo.jpg") pix = numpy.array(pic.getdata()).reshape(pic.size[0], pic.size[1], 3)   But I can't seem to figure out how to load it back into the PIL object after I've done all my awesome transforms.  I'm aware of the putdata() method, but can't quite seem to get it to behave.
109234
I am attempting to fully automate the process of georeferencing a load of .png files - i have xml file containing the co-ordinates - but how can i use them?  Does anyone have any advice?
54956
I am trying to export existing raster layers as tif files. I just want to avoid saving each raster manually. this is the code I am trying to use, but it doesn't seem to work:  from PyQt4.QtCore import * from qgis.core import * from qgis.gui import * import os import sys  layers = iface.layerTreeView().selectedLayers() for layer in layers:     output='Desktop'%layer.name%'.tif'   As specified in a comment below, I modified my code to this:  from PyQt4.QtCore import * from qgis.core import * from qgis.gui import * import os import sys   layers = iface.layerTreeView().selectedLayers()  for layer in layers:      file_name = 'D:\tif' + layer.name() + '.tif'      file_writer = QgsRasterFileWriter(file_name)      pipe = QgsRasterPipe()      provider = layer.dataProvider()       if not pipe.set(provider.clone()):         print "Cannot set pipe provider"      continue       file_writer.writeRaster(         pipe,         provider.xSize(),         provider.ySize(),         provider.extent(),         provider.crs())   I can't seem to get it to work and I am not getting any error.  I already selected the layers and it still does not work.
120719
In Python, using GDAL, I've extracted a raster's projection as a WKT string as follows:  wkt = dataset.GetProjection() # wkt is 'PROJCS["GDA_1994_Transverse_Mercator",GEOGCS["GDA_1994",DATUM["GDA_1994",SPHEROID["GRS_1980",6378137,298.2572221010002],TOWGS84[0,0,0,0,0,0,0]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",117],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",10000000],UNIT["Meter",1]]'   Using the WKT string, I can create a SpatialReference instance as follows:  src = osr.SpatialReference() src.ImportFromWkt(wkt)   This is easy-peasy. I can extract several parameters of the projection from src such as the UTM zone, etc., easily enough. But I can't figure out how to extract something like the name of the projection, i.e. "GDA_1994_Transverse_Mercator". This has surely got to be possible, but the Python API documentation may as well be nonexistent for all the use it is.  How do I extract the names of the projection and the geographic coordinate system?
89230
I have 1-band GeoTIFF file. The color interpretation for band is gray. I want a color table in order to transform my Col.tiff to Col.jpeg or Col.png.  I tried putting a color table in this tiff via a .vrt:  gdal_translate -of VRT  Col.tiff Col.vrt   then adding in the .vrt:  &lt;ColorInterp&gt;Palette&lt;/ColorInterp&gt; &lt;ColorTable&gt;    &lt;Entry c1="0" c2="0" c3="255" c4="240"/&gt;    &lt;Entry c1="0" c2="0" c3="254" c4="240"/&gt;    &lt;Entry c1="0" c2="0" c3="253" c4="240"/&gt;    &lt;Entry c1="0" c2="0" c3="252" c4="240"/&gt; &lt;/ColorTable&gt;   And translate again :  gdal_translate Col.vrt Col.tiff   But it answered:  Warning 1: Unable to export color table to GeoTIFF file.  Color tables can only be written to 1 band or 2 bands Byte or UInt16 GeoTIFF files.   How can I do it?  For information, there is gdalinfo for the Band:  Band 1 Block=55x18 Type=Float64, ColorInterp=Gray   Description = Something   Metadata:     GRIB_COMMENT=Something[p]     GRIB_ELEMENT=TMP     GRIB_FORECAST_SECONDS=0 sec     GRIB_REF_TIME=1401926400 sec UTC     GRIB_SHORT_NAME=SMTHG     GRIB_UNIT=[p]     GRIB_VALID_TIME=1401926400 sec UTC
16280
I was wondering how to check the value of a pixel at a mouse click in QGIS (as of 2.4).  In ArcGIS, the Identify tool does the job for both rasters and vectors. However, in QGIS, the Identify tool seems to be greyed out when only one raster layer is loaded.
6180
I have a GeoTIFF published by Geoserver 2.2 and being used by Leaflet via WMS. How can I get Geoserver to return the raster value at a lat/long point?   I'm guessing it has to do with using WCS, but no examples seem to exist for doing this!  This will be similar to querying PostGIS with ST_Value(raster, ST_SetSRID(ST_Point(lon,lat),4326)). I dont really want to import the same raster into PostGIS just to find the raster value at a point. Or is this the recommended solution?
83700
Is it possible to create buffers for pixels of raster files? Actually I need to expand the pixels classified as 1 (in white in the figure) with a spatial range of 1 or 2 pixels, in order to perform some generalization.  I'm using QGIS, is this possible with gdal or something?
20066
I have two raster files that are the same cell size but I cannot get them to align properly. (I want the black to line up with the pink)    When creating the black raster file (from using polyline to raster), I imported the pink raster in the Snap Raster section in the Environments setting. I thought that Snap Raster was supposed to make it line up?  Am I doing something wrong or is it a bug in Arc (I remember reading that this was a problem about two years ago), and if it is a bug is there a work around?  I have also read there is a Raster Shift, but I'm not sure what to enter in the X Coordinate and Y Coordinate. Since I only want the black raster to shift down, I know that the X Coor = 0 and the Y Coor should be negative something, but what are the units?  Thanks.