How to perform shapefile/raster integrated functions

Question

Description:

I'm having trouble with some Raster/Shapefile manipulation. My problem is that I have a shapefile that I would like to use to update my rastervalues with. The shapefile has multiple polygons (Countries), that I would like to to loop through and change the value of the cells that are masked by that shapefile. Most other questions are about clipping a raster with shapefile, but I just want the cells on the original raster to be masked. Here is an example of what I'm looking for:

for shapes in shapefile:
    masked_country=maskfunction(shapes,originalraster)
    masked_country=masked_country + arbitraryfunction(shape['layer'])

Do you have anything that works like that? I work with python and gdal.

Here are some of the other options I have tried:

Trying to rasterize shapefile using python gdal.RasterizeLayer:

I only get a black tiff using this method. I have even implemented the target_dst=None at the end, and I'm still getting a 0 value raster. Here is my code sample for this:

from osgeo import gdal #geo library
from osgeo import ogr

pathtofile1="rPath\to\my\shapefile"
pathtofile2="rPath\to\my\raster"

shapefile = ogr.Open(pathtofile1)
rasterfile=gdal.Open(pathtofile2)

mb_l = shapefile.GetLayer()
geo_transform = rasterfile.GetGeoTransform()
source_layer = rasterfile.GetLayer()

x_min = geo_transform[0]
y_max = geo_transform[3]
x_max = x_min + geo_transform[1] * rasterfile.RasterXSize
y_min = y_max + geo_transform[5] * rasterfile.RasterYSize
x_res = lur250.RasterXSize
y_res = lur250.RasterYSize
pixel_width = geo_transform[1]

target_ds = gdal.GetDriverByName('GTiff').Create('my.tiff', x_res, y_res, 1, gdal.GDT_Float32)
target_ds.SetGeoTransform((x_min, pixel_width, 0, y_max, 0, 0-pixel_width))

band = target_ds.GetRasterBand(1)
NoData_value = -999
band.SetNoDataValue(NoData_value)
band.FlushCache()

gdal.RasterizeLayer(target_ds, [1], mb_l,options=["ATTRIBUTE=OBJECTID"])

target_ds=None

There is a possibility that the "OBJECTID" is not Int and it's string, but I don't know how to check that.

Also I have played with the arguments in SetGeoTransform, but to no use (Setting Ymax or Ymin for different raaster definition).

Trying to burn the values of the shapefile into a raster with similar projection and resolution using gdal_rasterize from bash:

I need to create an empty raster with the same coordinate as my rasterfile (which can be achieved by copying that raster), but I don't know how to set the value for the entire raster to "0" or NoDataValue. All other solutions for these kind of problems are using "gdal_calc.py" and our HPC gdal installation is without such a tool.

Extra information

Here is the gdal_info results for my raster:

Driver: GTiff/GeoTIFF
Files: final4.tif
Size is 141969, 60829
Coordinate System is:
PROJCS["World_Mollweide",
    GEOGCS["GCS_WGS_1984",
        DATUM["D_WGS_1984",
            SPHEROID["WGS_1984",6378137.0,298.257223563]],
        PRIMEM["Greenwich",0.0],
        UNIT["Degree",0.017453292519943295]],
    PROJECTION["Mollweide"],
    PARAMETER["False_Easting",0.0],
    PARAMETER["False_Northing",0.0],
    PARAMETER["Central_Meridian",0.0],
    UNIT["Meter",1.0]]
Origin = (-17619594.546999998390675,8722279.461999999359250)
Pixel Size = (250.000000000000000,-250.000000000000000)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (-17619594.547, 8722279.462) ( 30d 5'16.03"E, 83d 8'14.83"N)
Lower Left  (-17619594.547,-6484970.538) (107d 3'58.73"E, 55d59' 2.97"S)
Upper Right (17872655.453, 8722279.462) ( 20d10'44.20"W, 83d 8'14.83"N)
Lower Right (17872655.453,-6484970.538) (103d26' 0.82"W, 55d59' 2.97"S)
Center      (  126530.453, 1118654.462) (  1d16'20.33"E,  9d 3'42.49"N)
Band 1 Block=141969x1 Type=Byte, ColorInterp=Gray
  NoData Value=-128

Are you opposed to using the rasterstats package? If not, you can feed it your shapefile and raster and it can return a masked array containing only the pixels within the shapefile polygons. Since that array will be the same shape as your input raster, you can update your raster with a one-liner. Then you have to know how to write that array to disk, but you can just copy all the georeferencing info from the original geotiff. — Jon, Jan 23 '19 at 23:42
Your code is almost identical to https://gis.stackexchange.com/questions/212795/rasterizing-shapefiles-with-gdal-and-python, is your vector in World_Mollweide projection? If not you will need to project the vector to the same CRS as the raster. To fill a raster with a single value (your NoData value) use band.Fill(-999) https://gdal.org/python/osgeo.gdal.Band-class.html before you RasterizeLayer. — Michael Stimson, Jan 24 '19 at 02:00
@Jon No I'm not opposed to them, but they usually have some basic operations such as Avg stats, but I need to apply a nonlinear function to the values in the masked raster. — Dr proctor, Jan 24 '19 at 20:19
Yes the code for rasterizing comes from there. I'll try both your answers to see if they work, and I will update the post. Thanks — Dr proctor, Jan 24 '19 at 20:20
Since the masked raster is returned, you can transform the pixel values however you want. — Jon, Jan 25 '19 at 22:10

score 1 · Answer 1 · answered Jan 28 '19 at 20:03

Thanks @jon and @michael-stimson for your answers. I found the solution. The problem was that I did not set the projection for my new file. I replicated the projection of the source file using these two lines:

proj=rasterfile.GetProjection()
target_ds.SetProjection(proj)

and then I filled the entire band with No-data value using Band.Fill() before doing the rasterization. I checked the results using gdalinfo -hist and it makes sense. Thanks a lot.

How to perform shapefile/raster integrated functions

1 Answers1