search

mattijn / KNMI-radar-PyGMT

3 stars 2 forks source link

gdal diffs #1

Open mattijn opened 3 years ago

mattijn commented 3 years ago

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 commented 3 years ago

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 commented 3 years ago

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 commented 3 years ago

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 commented 3 years ago

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 commented 3 years ago

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