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.

Sign up for GitHub

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

Jump to bottom

gdal diffs #1

Open
mattijn opened this issue Dec 1, 2020 · 5 comments
Open

gdal diffs #1

mattijn opened this issue Dec 1, 2020 · 5 comments

Comments

@mattijn
Copy link
Owner

mattijn commented Dec 1, 2020

Created radian tif:

gdalinfo using 2.4.4

Driver: GTiff/GeoTIFF
Files: RAD_NL25_PCP_NA_201910281900_rad.tif
Size is 700, 765
Coordinate System is:
PROJCS["unnamed",
    GEOGCS["unnamed ellipse",
        DATUM["unknown",
            SPHEROID["unnamed",6378.137,298.2528407762547]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433]],
    PROJECTION["Polar_Stereographic"],
    PARAMETER["latitude_of_origin",60],
    PARAMETER["central_meridian",0],
    PARAMETER["scale_factor",1],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]]]
Origin = (0.000000000000000,-3650.000000000000000)
Pixel Size = (1.000000000000000,-1.000000000000000)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (       0.000,   -3650.000) (  0d 0' 0.01"E, 55d58'24.82"N)
Lower Left  (       0.000,   -4415.000) (  0d 0' 0.01"E, 49d21'43.40"N)
Upper Right (     700.000,   -3650.000) ( 10d51'23.09"E, 55d23'20.17"N)
Lower Right (     700.000,   -4415.000) (  9d 0'33.39"E, 48d53'43.07"N)
Center      (     350.000,   -4032.500) (  4d57'37.96"E, 52d30'20.09"N)
Band 1 Block=700x2 Type=Float32, ColorInterp=Gray

gdalinfo using 3.x

Driver: GTiff/GeoTIFF
Files: RAD_NL25_PCP_NA_201910281900_rad.tif
Size is 700, 765
Coordinate System is:
PROJCRS["unnamed",
    BASEGEOGCRS["unnamed ellipse",
        DATUM["unknown",
            ELLIPSOID["unnamed",6378.137,298.252840776255,
                LENGTHUNIT["metre",1,
                    ID["EPSG",9001]]]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433,
                ID["EPSG",9122]]]],
    CONVERSION["Polar Stereographic (variant B)",
        METHOD["Polar Stereographic (variant B)",
            ID["EPSG",9829]],
        PARAMETER["Latitude of standard parallel",60,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8832]],
        PARAMETER["Longitude of origin",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8833]],
        PARAMETER["False easting",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["(E)",south,
            MERIDIAN[90,
                ANGLEUNIT["degree",0.0174532925199433,
                    ID["EPSG",9122]]],
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["(N)",south,
            MERIDIAN[180,
                ANGLEUNIT["degree",0.0174532925199433,
                    ID["EPSG",9122]]],
            ORDER[2],
            LENGTHUNIT["metre",1]]]
Data axis to CRS axis mapping: 1,2
Origin = (0.000000000000000,-3650.000000000000000)
Pixel Size = (1.000000000000000,-1.000000000000000)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (       0.000,   -3650.000) (  0d 0' 0.01"E, 55d58'24.82"N)
Lower Left  (       0.000,   -4415.000) (  0d 0' 0.01"E, 49d21'43.40"N)
Upper Right (     700.000,   -3650.000) ( 10d51'23.09"E, 55d23'20.17"N)
Lower Right (     700.000,   -4415.000) (  9d 0'33.39"E, 48d53'43.07"N)
Center      (     350.000,   -4032.500) (  4d57'37.96"E, 52d30'20.09"N)
Band 1 Block=700x2 Type=Float32, ColorInterp=Gray
@mattijn
Copy link
Owner Author

mattijn commented Dec 1, 2020

gdalwarp using 2.4.4

gdalwarp -t_srs '+proj=sterea +lat_0=52.15616055555555 +lon_0=5.38763888888889 +k=0.9999079 +x_0=155000 +y_0=463000 +ellps=bessel +units=m +no_defs +<>' RAD_NL25_PCP_NA_201910281900_rad.tif RAD_NL25_PCP_NA_201910281900_rd.tif

Creating output file that is 768P x 826L.
Processing RAD_NL25_PCP_NA_201910281900_rad.tif [1/1] : 0...10...20...30...40...50...60...70...80...90...100 - done.

gdalwarp using 2.4.4

gdalwarp -t_srs '+proj=sterea +lat_0=52.15616055555555 +lon_0=5.38763888888889 +k=0.9999079 +x_0=155000 +y_0=463000 +ellps=bessel +units=m +no_defs +<>' RAD_NL25_PCP_NA_201910281900_rad.tif RAD_NL25_PCP_NA_201910281900_rd.tif

Creating output file that is 768P x 826L.
Processing RAD_NL25_PCP_NA_201910281900_rad.tif [1/1] : 0...10...20...30...40...50...60...70...80...90...100 - done.

gdalwarp using 3.x

gdalwarp -t_srs '+proj=sterea +lat_0=52.15616055555555 +lon_0=5.38763888888889 +k=0.9999079 +x_0=155000 +y_0=463000 +ellps=bessel +units=m +no_defs +<>' RAD_NL25_PCP_NA_201910281900_rad.tif RAD_NL25_PCP_NA_201910281900_rd.tif

Processing RAD_NL25_PCP_NA_201910281900_rad.tif [1/1] : 0ERROR 1: PROJ: proj_create_operations: Source and target ellipsoid do not belong to the same celestial body
ERROR 6: Cannot find coordinate operations from `PROJCRS["unnamed",BASEGEOGCRS["unnamed ellipse",DATUM["unknown",ELLIPSOID["unnamed",6378.137,298.252840776255,LENGTHUNIT["metre",1,ID["EPSG",9001]]]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433,ID["EPSG",9122]]]],CONVERSION["Polar Stereographic (variant B)",METHOD["Polar Stereographic (variant B)",ID["EPSG",9829]],PARAMETER["Latitude of standard parallel",60,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8832]],PARAMETER["Longitude of origin",0,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8833]],PARAMETER["False easting",0,LENGTHUNIT["metre",1],ID["EPSG",8806]],PARAMETER["False northing",0,LENGTHUNIT["metre",1],ID["EPSG",8807]]],CS[Cartesian,2],AXIS["(E)",south,MERIDIAN[90,ANGLEUNIT["degree",0.0174532925199433,ID["EPSG",9122]]],ORDER[1],LENGTHUNIT["metre",1]],AXIS["(N)",south,MERIDIAN[180,ANGLEUNIT["degree",0.0174532925199433,ID["EPSG",9122]]],ORDER[2],LENGTHUNIT["metre",1]]]' to `PROJCRS["unknown",BASEGEOGCRS["unknown",DATUM["Unknown based on Bessel 1841 ellipsoid",ELLIPSOID["Bessel 1841",6377397.155,299.1528128,LENGTHUNIT["metre",1,ID["EPSG",9001]]]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8901]]],CONVERSION["unknown",METHOD["Oblique Stereographic",ID["EPSG",9809]],PARAMETER["Latitude of natural origin",52.1561605555556,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8801]],PARAMETER["Longitude of natural origin",5.38763888888889,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8802]],PARAMETER["Scale factor at natural origin",0.9999079,SCALEUNIT["unity",1],ID["EPSG",8805]],PARAMETER["False easting",155000,LENGTHUNIT["metre",1],ID["EPSG",8806]],PARAMETER["False northing",463000,LENGTHUNIT["metre",1],ID["EPSG",8807]]],CS[Cartesian,2],AXIS["(E)",east,ORDER[1],LENGTHUNIT["metre",1,ID["EPSG",9001]]],AXIS["(N)",north,ORDER[2],LENGTHUNIT["metre",1,ID["EPSG",9001]]]]'

@mattijn
Copy link
Owner Author

mattijn commented Dec 1, 2020

gdalwarp using 3.x including -ct command. See rfc73

gdalwarp -t_srs '+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs' -ct '+proj=axisswap +order=2,1 +step' myfile.tif myfile_out.tif

Creating output file that is 765P x 700L.
Processing myfile.tif [1/1] : 0...10...20...30...40...50...60...70...80...90...100 - done.

@mattijn
Copy link
Owner Author

mattijn commented Dec 2, 2020

environments setup

conda create --name gdal2 python=3.8
conda install jupyterlab
conda install gdal=2
conda install rasterio
conda install h5py
conda install pygmt --> conflicts

conda create --name gdal3 python=3.8
conda install jupyterlab
conda install gdal=3.1 # 3.2 is available
conda install rasterio
conda install h5py
conda install pygmt --> OK

conda create --name gdal32 python=3.8

@mattijn
Copy link
Owner Author

mattijn commented Dec 2, 2020

window
gdal3

(gdal3) D:\KNMI-radar-PyGMT>echo "0.000 -3650.000" | gdaltransform -s_srs '+proj=stere +x_0=0 +y_0=0 +lat_0=90 +lon_0=0 +lat_ts=60 +a=6378.137 +b=6356.752' -t_srs EPSG:4326 --debug on
ERROR 1: PROJ: proj_create: unrecognized format / unknown name
ERROR 1: Translating source or target SRS failed:
'+proj=stere
GDAL: In GDALDestroy - unloading GDAL shared library.

gdal2

(gdal2) D:\KNMI-radar-PyGMT>echo "0.000 -3650.000" | gdaltransform -s_srs '+proj=stere +x_0=0 +y_0=0 +lat_0=90 +lon_0=0 +lat_ts=60 +a=6378.137 +b=6356.752' -t_srs EPSG:4326 --debug on
OGRCT: PROJ >= 4.8.0 features enabled
OGRCT: Using locale-safe proj version
OGR_PROJ4: Can't find +proj= in:
'+proj=stere
ERROR 1: Translating source or target SRS failed:
'+proj=stere
GDAL: In GDALDestroy - unloading GDAL shared library.

@mattijn
Copy link
Owner Author

mattijn commented Dec 2, 2020 •
edited
Loading

macos
gdal3

(gdal3) mattijnvanhoek@Mattijns-MacBook-Pro ~ % echo "0.000 -3650.000" | gdaltransform -s_srs '+proj=stere +x_0=0 +y_0=0 +lat_0=90 +lon_0=0 +lat_ts=60 +a=6378.137 +b=6356.752' -t_srs EPSG:4326 -ct '+proj=axisswap +order=2,1 +step' --debug on

0 -3650 0

gdal2

(gdal2) mattijnvanhoek@Mattijns-MacBook-Pro ~ % echo "0.000 -3650.000" | gdaltransform -s_srs '+proj=stere +x_0=0 +y_0=0 +lat_0=90 +lon_0=0 +lat_ts=60 +a=6378.137 +b=6356.752' -t_srs EPSG:4326 -ct '+proj=axisswap +order=2,1 +step' --debug on
Usage: gdaltransform [--help-general]
    [-i] [-s_srs srs_def] [-t_srs srs_def] [-to "NAME=VALUE"]
    [-order n] [-tps] [-rpc] [-geoloc] 
    [-gcp pixel line easting northing [elevation]]* [-output_xy]
    [srcfile [dstfile]]


FAILURE: Unknown option name '-ct'

(gdal2) mattijnvanhoek@Mattijns-MacBook-Pro ~ % echo "0.000 -3650.000" | gdaltransform -s_srs '+proj=stere +x_0=0 +y_0=0 +lat_0=90 +lon_0=0 +lat_ts=60 +a=6378.137 +b=6356.752' -t_srs EPSG:4326 --debug on  
OGRCT: Source: +proj=stere +lat_0=90 +lat_ts=60 +lon_0=0 +k=1 +x_0=0 +y_0=0 +a=6378.137 +b=6356.752 +units=m +no_defs
OGRCT: Target: +proj=longlat +datum=WGS84 +no_defs
OGRCT: Source: +proj=longlat +datum=WGS84 +no_defs
OGRCT: Target: +proj=stere +lat_0=90 +lat_ts=60 +lon_0=0 +k=1 +x_0=0 +y_0=0 +a=6378.137 +b=6356.752 +units=m +no_defs
0 55.9735620711571 0

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Assignees
No one assigned
Labels
None yet
Projects
None yet
Milestone
No milestone
Development

No branches or pull requests
1 participant
@mattijn