Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

GeoArrays produces positively y spaced rasters by default #53

Open
kongdd opened this issue Oct 2, 2020 · 4 comments
Open

GeoArrays produces positively y spaced rasters by default #53

kongdd opened this issue Oct 2, 2020 · 4 comments

Comments

@kongdd
Copy link

kongdd commented Oct 2, 2020

I am writing tiff using the following scripts:

begin
    ga = GeoArray(rand(360, 150))
    bbox!(ga, (min_x=-180.0, min_y=-60.0, max_x=180.0, max_y=90.0))
    epsg!(ga, 4326)  # in WGS84
    GeoArrays.write!("test.tif", ga)    
end

When reading in R, the crsTransform is error. The latitude range is changed to [-59, 91].

library(raster)
b <- brick("test.tif")
Warning message:
In .rasterFromGDAL(x, band = band, objecttype, ...) :
  data seems flipped. Consider using: flip(x, direction='y')

b
class      : RasterBrick 
dimensions : 150, 360, 54000, 1  (nrow, ncol, ncell, nlayers)
resolution : 1, 1  (x, y)
extent     : -180, 180, -59, 91  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +datum=WGS84 +no_defs 
source     : N:/Research/GEE_repos/GEE-latest/test.tif 
names      : layer 
@kongdd
Copy link
Author

kongdd commented Oct 2, 2020

This problem can be solve by fixing bbox_to_affine.
jl-pkgs@2cc1c09

@evetion
Copy link
Owner

evetion commented Oct 2, 2020

Thanks for reporting! I agree that the data is flipped in the GDAL ecosystem sense, as it expects a negative y spacing.

But it shouldn't be wrong, especially not in terms of things like latitude:

❯ gdalinfo test.tif
Driver: GTiff/GeoTIFF
Files: test.tif
Size is 360, 150
Coordinate System is:
GEOGCRS["WGS 84",
    DATUM["World Geodetic System 1984",
        ELLIPSOID["WGS 84",6378137,298.257223563,
            LENGTHUNIT["metre",1]]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433]],
    CS[ellipsoidal,2],
        AXIS["geodetic latitude (Lat)",north,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433]],
        AXIS["geodetic longitude (Lon)",east,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433]],
    ID["EPSG",4326]]
Data axis to CRS axis mapping: 2,1
Origin = (-180.000000000000000,-60.000000000000000)
Pixel Size = (1.000000000000000,1.000000000000000)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  COMPRESSION=DEFLATE
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (-180.0000000, -60.0000000) (180d 0' 0.00"W, 60d 0' 0.00"S)
Lower Left  (-180.0000000,  90.0000000) (180d 0' 0.00"W, 90d 0' 0.00"N)
Upper Right ( 180.0000000, -60.0000000) (180d 0' 0.00"E, 60d 0' 0.00"S)
Lower Right ( 180.0000000,  90.0000000) (180d 0' 0.00"E, 90d 0' 0.00"N)
Center      (   0.0000000,  15.0000000) (  0d 0' 0.01"E, 15d 0' 0.00"N)
Band 1 Block=256x256 Type=Float64, ColorInterp=Gray

It boils down to how you would interpret an Array in Julia. Is it positively spaced? For your fix to work, the data needs to be flipped upside down too! Note there's a GeoArrays.flipud! utility exactly for this. I see a few possibilities:

  • We give a warning only when positive y spacing is encountered, hinting to use flipud!.
  • Always use flipud!, to arrive at the negative y spacing.

@visr Any thoughts on this?

@evetion evetion changed the title crsTransform error GeoArrays produces positively y spaced rasters by default Oct 2, 2020
@visr
Copy link
Collaborator

visr commented Oct 4, 2020

GDAL seems to support both positive and negative y spacing, so I think the current behaviour is fine? Although negative y spacing seems to be the norm. Are there many applications that only work with negative cell spacing? Does QGIS show it well?

@evetion
Copy link
Owner

evetion commented Oct 4, 2020

Showing is often no problem, but some GDAL command line tools refuse to work with them.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

3 participants