search

boku-ilen / geodot-plugin

Godot plugin for loading geospatial data
GNU General Public License v3.0
108 stars 18 forks source link

Can I render non-square geotiff? #62

Closed VRichardJP closed 2 years ago

VRichardJP commented 3 years ago

I could successfully run the demo. I also had some success with a custom geotiff file. However the area I want to render is not a square (1242x648) and if I modify the demo code to try to display the full image I get an error like this:

$ godot .
Godot Engine v3.3.stable.arch_linux - https://godotengine.org
OpenGL ES 3.0 Renderer: AMD RENOIR (DRM 3.40.0, 5.11.16-zen1-1-zen, LLVM 11.1.0)
OpenGL ES Batching: ON

ERROR 5: out1.tif, band 1: Access window out of range in RasterIO().  Requested
(0,-594) of size 647x647 on raster of 1242x648.
ERROR 5: out1.tif, band 2: Access window out of range in RasterIO().  Requested
(0,-594) of size 647x647 on raster of 1242x648.
ERROR 5: out1.tif, band 3: Access window out of range in RasterIO().  Requested
(0,-594) of size 647x647 on raster of 1242x648.
ERROR 5: out1.tif, band 4: Access window out of range in RasterIO().  Requested
(0,-594) of size 647x647 on raster of 1242x648.
ERROR: create_from_image: Invalid image
   At: scene/resources/texture.cpp:199.
^C

Is there any support for non-square image?

info about the geotiff:

gdalinfo out1.tif 
Driver: GTiff/GeoTIFF
Files: out1.tif
Size is 1242, 648
Coordinate System is:
PROJCRS["WGS 84 / Pseudo-Mercator",
    BASEGEOGCRS["WGS 84",
        ENSEMBLE["World Geodetic System 1984 ensemble",
            MEMBER["World Geodetic System 1984 (Transit)"],
            MEMBER["World Geodetic System 1984 (G730)"],
            MEMBER["World Geodetic System 1984 (G873)"],
            MEMBER["World Geodetic System 1984 (G1150)"],
            MEMBER["World Geodetic System 1984 (G1674)"],
            MEMBER["World Geodetic System 1984 (G1762)"],
            ELLIPSOID["WGS 84",6378137,298.257223563,
                LENGTHUNIT["metre",1]],
            ENSEMBLEACCURACY[2.0]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",4326]],
    CONVERSION["Popular Visualisation Pseudo-Mercator",
        METHOD["Popular Visualisation Pseudo Mercator",
            ID["EPSG",1024]],
        PARAMETER["Latitude of natural origin",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["False easting",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["easting (X)",east,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["northing (Y)",north,
            ORDER[2],
            LENGTHUNIT["metre",1]],
    USAGE[
        SCOPE["Web mapping and visualisation."],
        AREA["World between 85.06°S and 85.06°N."],
        BBOX[-85.06,-180,85.06,180]],
    ID["EPSG",3857]]
Data axis to CRS axis mapping: 1,2
Origin = (15208294.336400000378489,4159389.058000000193715)
Pixel Size = (1.000020450885421,-0.999814506173191)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  COMPRESSION=PACKBITS
  INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left  (15208294.336, 4159389.058) (136d37' 6.36"E, 34d58' 0.98"N)
Lower Left  (15208294.336, 4158741.178) (136d37' 6.36"E, 34d57'43.81"N)
Upper Right (15209536.362, 4159389.058) (136d37'46.52"E, 34d58' 0.98"N)
Lower Right (15209536.362, 4158741.178) (136d37'46.52"E, 34d57'43.81"N)
Center      (15208915.349, 4159065.118) (136d37'26.44"E, 34d57'52.39"N)
Band 1 Block=1242x1 Type=Byte, ColorInterp=Red
  Mask Flags: PER_DATASET ALPHA 
Band 2 Block=1242x1 Type=Byte, ColorInterp=Green
  Mask Flags: PER_DATASET ALPHA 
Band 3 Block=1242x1 Type=Byte, ColorInterp=Blue
  Mask Flags: PER_DATASET ALPHA 
Band 4 Block=1242x1 Type=Byte, ColorInterp=Alpha
kb173 commented 3 years ago

This is related to a GDAL limitation - its RasterIO function, which we use for loading data, doesn't support loading sections which contain pixels for which no data is available. I'm currently working on a workaround for this, so that Geodot returns an image with black pixels for such sections. With that, it should be easy to just load a square which encompasses all data, and then only use the part in which data is available. Directly fetching non-square segments as GeoImages is not something we're planning to support, but this workaround should solve most problems which arise from this limitation, so I'll keep this open until that's done.

However, it should generally be no problem to load a square window from a non-square dataset, as long as this square window is fully covered with available data. It says that you're requesting '(0,-594) of size 647x647 on raster of 1242x648.'. I'm assuming the y coordinate of the origin point (-594) is the problem here - are you sure the projected coordinates you put into the get_image call correspond to the top-left corner of the image?

VRichardJP commented 3 years ago

In the error case I just tried to set the tile size to 1242, which I guess expects a 1242x1242 image :^) https://github.com/boku-ilen/geodot-plugin/blob/f16b7c7c437a0410e0580297bccb2c0551ba8b3d/demo/RasterDemo.gd#L6

kb173 commented 3 years ago

Actually this size is independent of the source resolution, the data is scaled up or down depending on the resolution of the data vs. the requested resolution. You could theoretically request a 10000x10000 image from your data, it'll just be very blurry. The meter size of the window (in the line above the one you posted) is important here, as that needs to be within the available data until the workaround I'm building is done. But according to the error message you got, your meter size was fine, just the top-left y-coordinate was outside of the available data. Maybe try adding 594 to your start_position_y?

VRichardJP commented 3 years ago

Sorry these parameters are working:

export(float) var start_position_x = 15208294.336400000378489
export(float) var start_position_y = 4159389.058000000193715
export(float) var tile_size_meters = 648
export(int) var tile_size_pixels = 648

I have set these parameters based on Size, Origin and Pixel Size (1px = 1m) output from gdalinfo (see first comment) The non working case is with the following parameters:

export(float) var tile_size_meters = 1242
export(int) var tile_size_pixels = 1242

In such case, moving start_position by (0, 594) fix the demo:

export(float) var start_position_x = 15208294.336400000378489
export(float) var start_position_y = 4159389.058000000193715 + 594
export(float) var tile_size_meters = 1242
export(int) var tile_size_pixels = 1242
kb173 commented 2 years ago

The problem described in my first comment (the inability to load sections not fully covered by data) has been fixed in the last few commits: it's now possible to load any section, even if there is no data at all, and the parts without coverage will simply be 0.

In addition, https://github.com/boku-ilen/geodot-plugin/commit/817a65b2e46cb3769e53b28ccb626597afde43f0 added the function get_center() to the GeoRasterLayer - this makes it easier to place oneself within the data without having to use external tools like gdalinfo.