export in wgs84

[bhbraswell]

I can't remember if there is a trick to this. I tried several incantations, and even searched the issues, read the manual, etc.

I know that using -r [raster template] will work but that's not really practical for my current use case.

Example:

gips_export cdl -p cdl -d 2008 --days 1,1 -s /export/tiles.shp -w "tileid='548_341'" -v4 --fetch --outdir /export/cdl_548_341 --notld --res 0.00005 0.00005

Since the shapefile is epsg:4326, I naively expected that this would produce an image 0.15x0.15 degrees in extent (the size of a tile) and 3000x3000 pixels (the right number of pixels at that resolution).

What came out had the correct resolution, but the extent was bonkers - it wasn't the extent of the tile, or of the shapefile (which I initially thought had happened). I assume the bonkersness has to do with gdal_warp doing something or not doing something. Then I tried some other stuff that didn't work, but remembered that I'd been through all this at some point.

Is this a thing / a thing that can be accomplished?

[ircwaves]

Not having your files, I just did the following (clipped for brevity):

(venv) icooke@north:~$ gips_export cdl -p cdl -d 2017 --days 1,1 -s src/gips/gips/data/sentinel2/tiles.shp -w "Name='18TYP'" --outdir ./cdl4326 --notld --res 0.00005 0.00005
GIPS Data Export (v0.13.0-dev)
  Dates: 1 dates (2017-01-01 - 2017-01-01)
  Products: cdl
    DATE     Coverage  Product  
2017        
    001     cdl  


1 files on 1 dates
(venv) icooke@north:~$ ogr2ogr  -where "Name='18TYP'" 18TYP_only.shp src/gips/gips/data/sentinel2/tiles.shp
(venv) icooke@north:~$ ogrinfo -so -al 18TYP_only.shp
INFO: Open of `18TYP_only.shp'
      using driver `ESRI Shapefile' successful.

Layer name: 18TYP_only
Metadata:
  DBF_DATE_LAST_UPDATE=2019-01-17
Geometry: Polygon
Feature Count: 1
Extent: (-72.537265, 43.201211) - (-71.124261, 44.225975)
Layer SRS WKT:
GEOGCS["GCS_WGS_1984",
    DATUM["WGS_1984",
        SPHEROID["WGS_84",6378137,298.257223563]],
    PRIMEM["Greenwich",0],
    UNIT["Degree",0.017453292519943295],
    AUTHORITY["EPSG","4326"]]
   ...
(venv) icooke@north:~$ gdalinfo -noct cdl4326/16982/2017001_cdl_cdl.tif 
Driver: GTiff/GeoTIFF
Files: cdl4326/16982/2017001_cdl_cdl.tif
Size is 48261, 40496
Coordinate System is:
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"]]
Origin = (-72.537264522300006,44.225975127500000)
Pixel Size = (0.000049999998737,-0.000049999998737)
...
Corner Coordinates:
Upper Left  ( -72.5372645,  44.2259751) ( 72d32'14.15"W, 44d13'33.51"N)
Lower Left  ( -72.5372645,  42.2011752) ( 72d32'14.15"W, 42d12' 4.23"N)
Upper Right ( -70.1242146,  44.2259751) ( 70d 7'27.17"W, 44d13'33.51"N)
Lower Right ( -70.1242146,  42.2011752) ( 70d 7'27.17"W, 42d12' 4.23"N)
Center      ( -71.3307396,  43.2135752) ( 71d19'50.66"W, 43d12'48.87"N)
...

I think somewhere, gippy is calculating the extent incorrectly by 1 your-shapefiles-units, which doesn't even register as an issue when the units are meters, but is a problem when the units are degrees:

[ircwaves]

Given that this is more of a gippy==0.3.x issue, I'm going to try the GH issue transfer feature and bump this over to there.

[bhbraswell]

Thanks for confirming that it is a thing, just maybe not a GIPS thing. I feel like a temporary work-around shouldn't be too hard e.g. a couple GDAL system calls. I'll look into it.

[ircwaves]

I've identified the issue, the extent width calculation:

            T width() const { return _p1.x()-_p0.x()+1; }

Seems like an assumption that the extent is ok a unit grid.

[ircwaves]

Looks like this has been addressed in the 1.0 branch. Really makes me want to get some time to do that port to gippy>=1.0.