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
