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)