0

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:

  1. 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).

  1. 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 at 23:42
  • Your code is almost identical to gis.stackexchange.com/questions/212795/…, 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) gdal.org/python/osgeo.gdal.Band-class.html before you RasterizeLayer. – Michael Stimson Jan 24 at 2: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 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 at 20:20
  • Since the masked raster is returned, you can transform the pixel values however you want. – Jon Jan 25 at 22:10
1

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.

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.