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 y_max = geo_transform x_max = x_min + geo_transform * rasterfile.RasterXSize y_min = y_max + geo_transform * rasterfile.RasterYSize x_res = lur250.RasterXSize y_res = lur250.RasterYSize pixel_width = geo_transform 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, , 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.
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