1

I want to convert a shapefile (polygons in WGS84) which contains 3 variables (float values) for each municipality of an area. In fact it contains 1035 municipalities. I want to convert it to a raster file and distribute the values of municipalities to a specific grid with LCC projection .

After the suggestions described below in the answers/comments (to my initial question which I have edited) I reprojected the original shapefile with the municipalities to the final desired grid , I made the calculations with the area and then I rasterised the file as following:

ogr2ogr -t_srs "$outproj" emissions_municip_reproj.shp emissions_municip.shp

ogrinfo -dialect SQLite -sql "SELECT Esye_id, CO,NOx,SOx ST_Area(geometry) as AREA FROM emissions_municip_reproj" emissions_municip_reproj.shp

INFO: Open of $Database_outdir/emissions_municip_reproj.shp

# Add field, calculate and rasterize the shapefile
ogrinfo emissions_municip_reproj.shp -sql "ALTER TABLE emissions_gr_Cday_municip_reproj ADD COLUMN AREA float"

ogrinfo emissions_municip_reproj.shp -dialect SQLite -sql "UPDATE emissions_municip_reproj SET AREA =ST_Area(geometry)"
ogrinfo emissions_municip_reproj.shp -dialect SQLite -sql "UPDATE emissions_municip_reproj SET COarea = CO/ST_Area(geometry) * $RESX * $RESY"

gdal_rasterize -a COarea -ot Float32 -tr $RESX $RESY -te $xmin $ymin $xmax $ymax -a_nodata 0 -l emissions_municip_reproj emissions_municip_reproj.shp emissionsCO.tif

The result is almost correct. The problem is that it seems that values have not been distributed to all cells. In the picture below the municipality is shown with red line. It seems that over some part of the municipality there are no grid cells (no data).Any idea for this? enter image description here

3

You would need to calculate the appropriate value before you rasterized, by dividing your values by the polygon area then multiplying by grid cell area.

Steps:

  • add a new field,
  • calculate it to be CO/municipality area * 1000 * 1000.
  • rasterize

For example: if Municipality X had a total of 100kg CO and and area of 10000000, your COarea field would get a value of 100 / 10000000 * 1000 * 1000 = 10kg, so your 10 output grid cells would each get a value of 10.

A worked example with a shapefile called /tmp/test.shp:

# Original shapefile
ogrinfo -dialect SQLite -sql "SELECT ID, CO, ST_Area(geometry) as AREA FROM test" /tmp/test.shp
INFO: Open of `/tmp/test.shp'
      using driver `ESRI Shapefile' successful.

Layer name: SELECT
Geometry: None
Feature Count: 1
Layer SRS WKT:
(unknown)
id: Integer64 (0.0)
CO: Real (0.0)
AREA: Real (0.0)
OGRFeature(SELECT):0
  id (Integer64) = 1
  CO (Real) = 250
  AREA (Real) = 1010187004.92576

# Add field, calculate and rasterize 
ogrinfo /tmp/test.shp -sql "ALTER TABLE test ADD COLUMN COarea float"
ogrinfo /tmp/test.shp -dialect SQLite -sql "UPDATE test SET COarea = CO/ST_Area(geometry) * 1000 * 1000"
gdal_rasterize -a COarea -ot Float32 -ts 1000 1000 -a_nodata 0 -l test /tmp/test.shp /tmp/test.tif

# Output TIFF
gdalinfo -stats /tmp/test.tif
Driver: GTiff/GeoTIFF
Files: /tmp/test.tif
Size is 1000, 1000
<snip>
Band 1 Block=1000x2 Type=Float32, ColorInterp=Gray
  Minimum=0.247, Maximum=0.247, Mean=0.247, StdDev=0.000
  NoData Value=0
  Metadata:
    STATISTICS_MAXIMUM=0.24747893214226
    STATISTICS_MEAN=0.24747893214226
    STATISTICS_MINIMUM=0.24747893214226
    STATISTICS_STDDEV=0

# Modifed shapefile
ogrinfo -fields=YES -geom=SUMMARY /tmp/test.shp test
INFO: Open of `/tmp/test.shp'
      using driver `ESRI Shapefile' successful.

Layer name: test
Metadata:
  DBF_DATE_LAST_UPDATE=2019-02-02
Geometry: Polygon
Feature Count: 1
Extent: (1253387.735842, -4737341.287382) - (1298147.948308, -4702942.793010)
<snip>
id: Integer64 (18.0)
CO: Real (24.15)
COarea: Real (24.15)
OGRFeature(test):0
  id (Integer64) = 1
  CO (Real) = 250.000000000000000
  COarea (Real) = 0.247478930911780
  POLYGON : 7 points
  • thanks for the useful answer. I supposed that you mean that I had to divide with 1000*1000, right? Because in the code above it seems that you multiply with 1000*1000. Also, is it possible to rasterise in a specific grid with another projection and resolution than that of the origin shapefile? In fact, I want to convert it to a raster with 2000m resolution in a lcc projection. Is it possible to define the projection in the gdal_rasterize ? – Nat Feb 3 at 13:08
  • Regarding the second part of my question, I created a new question here gis.stackexchange.com/questions/310817/…. I think that working only with Shapefiles instead of converting a raster would be more easy to solve the problem with the reprojection and distribution of the variable values. – Nat Feb 3 at 14:46
  • @Nat - No I did not mean divide by 1000x1000 (the grid cell area). You divide by the area of the polygon then multiply by the area of a grid cell. – user2856 Feb 3 at 19:03
  • in this way , the final value of each grid cell takes very large value (much larger than the expected). Maybe the problem is the unit of the polygon area? Do I have to define it in someway? The shapefile is in WGS84. For example, for one municipality I see in GIS that the corresponding area is 20000000m2 and the parameter ST_Area(geometry) in gdal takes the value of ~0.002 – Nat Feb 3 at 19:32
  • Yes, your coordinate units need to be in the same units as your raster. You can either reproject the shapefile or create a field that stores the area in m2 and divide by that instead of ST_area – user2856 Feb 3 at 19:39

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.