1

I'm asking here, because, I've got great results last time. Here it is the next of my adventure with gdal_rasterize, I want to make a raster picture where each pixel correspond to a square of 100 by 100 meters from a table of points of a Postgis DB (EPSG:3857).

From the manual : http://www.gdal.org/gdal_rasterize.html I found the option -tr normally is made to do it. I'm executing the following command :

gdal_rasterize -a value -tr 100 100 PG:"host=localhost dbname=myDB user=postgres" -sql "SELECT the_geom, value FROM table" out.tif

The resulting tiff is an irregular raster : sometime pixel size correspond to 50m sometime 150m (I calculate that thanks to the rule of Google earth)

What I did wrong ?

[Edit] As asked in comment, here it is the content of gdalinfo of the resulted tiff.

gdalinfo content for the file generated with the option "-tr 50 50" :

Driver: GTiff/GeoTIFF
Files: out_50_50.tif
Size is 22018, 15699
Coordinate System is:
PROJCS["WGS 84 / Pseudo-Mercator",
    GEOGCS["WGS 84",
        DATUM["WGS_1984",
            SPHEROID["WGS 84",6378137,298.257223563,
                AUTHORITY["EPSG","7030"]],
            AUTHORITY["EPSG","6326"]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433],
        AUTHORITY["EPSG","4326"]],
    PROJECTION["Mercator_1SP"],
    PARAMETER["central_meridian",0],
    PARAMETER["scale_factor",1],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    EXTENSION["PROJ4","+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +
x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext  +no_defs"],
    AUTHORITY["EPSG","3857"]]
Origin = (-1045006.900697760100000,5457301.723357450200000)
Pixel Size = (50.000000000000000,-50.000000000000000)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (-1045006.901, 5457301.723) (  9d23'14.84"W, 43d56'50.54"N)
Lower Left  (-1045006.901, 4672351.723) (  9d23'14.84"W, 38d39'17.46"N)
Upper Right (   55893.099, 5457301.723) (  0d30' 7.55"E, 43d56'50.54"N)
Lower Right (   55893.099, 4672351.723) (  0d30' 7.55"E, 38d39'17.46"N)
Center      ( -494556.901, 5064826.723) (  4d26'33.65"W, 41d21'17.53"N)
Band 1 Block=22018x1 Type=Float64, ColorInterp=Gray

gdalinfo content for the file generated with the option "-tr 150 150" :

Driver: GTiff/GeoTIFF
Files: out_150_150.tif
Size is 7339, 5233
Coordinate System is:
PROJCS["WGS 84 / Pseudo-Mercator",
    GEOGCS["WGS 84",
        DATUM["WGS_1984",
            SPHEROID["WGS 84",6378137,298.257223563,
                AUTHORITY["EPSG","7030"]],
            AUTHORITY["EPSG","6326"]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433],
        AUTHORITY["EPSG","4326"]],
    PROJECTION["Mercator_1SP"],
    PARAMETER["central_meridian",0],
    PARAMETER["scale_factor",1],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    EXTENSION["PROJ4","+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +
x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext  +no_defs"],
    AUTHORITY["EPSG","3857"]]
Origin = (-1045006.900697760100000,5457301.723357450200000)
Pixel Size = (150.000000000000000,-150.000000000000000)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (-1045006.901, 5457301.723) (  9d23'14.84"W, 43d56'50.54"N)
Lower Left  (-1045006.901, 4672351.723) (  9d23'14.84"W, 38d39'17.46"N)
Upper Right (   55843.099, 5457301.723) (  0d30' 5.93"E, 43d56'50.54"N)
Lower Right (   55843.099, 4672351.723) (  0d30' 5.93"E, 38d39'17.46"N)
Center      ( -494581.901, 5064826.723) (  4d26'34.46"W, 41d21'17.53"N)
Band 1 Block=7339x1 Type=Float64, ColorInterp=Gray

Is my projection well suited ?

  • 1
    just a guess, but I would write the -sql " " before the input name PG:" " . Maybe this is confusing gdal. – radouxju Apr 11 '14 at 14:24
  • 1
    I'm not understand, the following command 'gdal_rasterize -a value -tr 100 100 -sql PG:"host=localhost dbname=myDB user=postgres" "SELECT the_geom, value FROM table" out.tif' is not working – Jérôme B Apr 11 '14 at 14:34
  • 2
    I meant gdal_rasterize -a value -tr 100 100 -sql "SELECT the_geom, value FROM table" PG:"host=localhost dbname=myDB user=postgres" out.tif – radouxju Apr 11 '14 at 14:37
  • 2
    Could you edit your question, showing the output of gdalinfo out.tif for both the 150m and 50m rasters? – Simbamangu Apr 12 '14 at 10:37
  • 2
    Are some of your points far from the equator? If you look at how the World looks in EPSG:3857 openstreetmap.org/#map=1/51/7 you can see that pixels near the poles are much bigger in real world units than near the equator. If you want the size of the pixel to correspond better with reality you must use some other projection, for example UTM zone that fits the point location. – user30184 May 15 '14 at 7:23

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.