2

I am clipping raster by mask with python. My input raster owned Float32, however the output of clipping was Byte.

I don't know why? I suspect that transformation of tif to numpy.array is changing them.

I would like that output will have the same values from its origin image.

Here my code:

"""Clip a raster image using a shapefile"""

import gdal, gdalnumeric
import shapefile
import Image, ImageDraw

# Raster image to clip
raster = "/home/infra/PycharmProjects/map_ailton/dados/raster/resultado/custos_it_sema_2015xmunicipio_uc_multiplicacao.tif"

# Polygon shapefile used to clip
shp = "/home/infra/PycharmProjects/map_ailton/dados/vector/chapada/chapada_square"

# Name of clipped raster file(s)
output = "/home/infra/PycharmProjects/map_ailton/script/4558_06_Codes/output"

def imageToArray(i):
    """
    Converts a Python Imaging Library array to a gdalnumeric image.
    """
    a=gdalnumeric.numpy.fromstring(i.tostring(),'b')
    a.shape=i.im.size[1], i.im.size[0]
    return a

def world2Pixel(geoMatrix, x, y):
  """
  Uses a gdal geomatrix (gdal.GetGeoTransform()) to calculate
  the pixel location of a geospatial coordinate 
  """
  ulX = geoMatrix[0]
  ulY = geoMatrix[3]
  xDist = geoMatrix[1]
  yDist = geoMatrix[5]
  rtnX = geoMatrix[2]
  rtnY = geoMatrix[4]
  pixel = int((x - ulX) / xDist)
  line = int((ulY - y) / xDist)
  return (pixel, line) 

# Load the source data as a gdalnumeric array
srcArray = gdalnumeric.LoadFile(raster)

# Also load as a gdal image to get geotransform (world file) info
srcImage = gdal.Open(raster)
geoTrans = srcImage.GetGeoTransform()

# Use pyshp to open the shapefile
r = shapefile.Reader("%s.shp" % shp)

# Convert the layer extent to image pixel coordinates
minX, minY, maxX, maxY = r.bbox
ulX, ulY = world2Pixel(geoTrans, minX, maxY)
lrX, lrY = world2Pixel(geoTrans, maxX, minY)

# Calculate the pixel size of the new image
pxWidth = int(lrX - ulX)
pxHeight = int(lrY - ulY)

# Multi-band image?
#Check this modification in script in: http://karthur.org/2015/clipping-rasters-in-python.html
try:
    clip = srcArray[:, ulY:lrY, ulX:lrX]

# Nope: Must be single-band
except IndexError:
    clip = srcArray[ulY:lrY, ulX:lrX]


# Create a new geomatrix for the image
geoTrans = list(geoTrans)
geoTrans[0] = minX
geoTrans[3] = maxY

# Map points to pixels for drawing the county boundary 
# on a blank 8-bit, black and white, mask image.
pixels = []
for p in r.shape(0).points:
  pixels.append(world2Pixel(geoTrans, p[0], p[1]))
rasterPoly = Image.new("L", (pxWidth, pxHeight), 1)
# Create a blank image in PIL to draw the polygon.
rasterize = ImageDraw.Draw(rasterPoly)
rasterize.polygon(pixels, 0)
# Convert the PIL image to a NumPy array
mask = imageToArray(rasterPoly)   

# Clip the image using the mask
clip = gdalnumeric.numpy.choose(mask, (clip, 0)).astype(gdalnumeric.numpy.uint8)

# Save clipping as tiff
gdalnumeric.SaveArray(clip, "%s.tif" % output, format="GTiff", prototype=raster)
3

The reason your output image is 8-bit is because you explicitly ensure the image is 8-bit the step before you write it.

# Clip the image using the mask
clip = gdalnumeric.numpy.choose(mask, (clip,0)).astype(gdalnumeric.numpy.uint8)

# Save clipping as tiff
gdalnumeric.SaveArray(clip, "%s.tif" % output, format="GTiff", prototype=raster)

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.