TransWikia.com

Going from degrees to meters in projection conversion using GDAL?

Geographic Information Systems Asked on July 14, 2021

I have raster data in a Lambert-Conformal conic projection created from a NetCDF file by running sub-processes in a Python app that runs the following (Linux) command line argument:

gdal_translate -ot "Float32" -of "GTiff" -a_srs ${MY_WKT} nc_file.nc tif_file.tif

where:

]$ echo ${M_WKT}
]$ PROJCS["unnamed",
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"]],
PROJECTION["Lambert_Conformal_Conic_2SP"],
PARAMETER["standard_parallel_1",13.34249973297119],
PARAMETER["standard_parallel_2",47.31750106811523],
PARAMETER["latitude_of_origin",30.33000183105469],
PARAMETER["central_meridian",-86.16000366210938],
PARAMETER["false_easting",0],
PARAMETER["false_northing",0],
UNIT["metre",1,
    AUTHORITY["EPSG","9001"]]]

and a snippet of the "gdalinfo" dump is:

Attributes:
grid_mapping_name:              lambert_conformal_conic
earth_radius:                   6370000.0
standard_parallel:              (13.3425, 47.3175)
longitude_of_central_meridian:  -86.16
latitude_of_projection_origin:  30.330002


pressure_interp#coordinates=XLAT Time level XTIME XLONG
  pressure_interp#FieldType=104
  pressure_interp#grid_mapping=LambertConformal
  pressure_interp#missing_value=9.969209968386869e+36
  pressure_interp#units=hPa
  pressure_interp#vert_units=m
  pressure_interp#_FillValue=9.96921e+36
  south_north#axis=Y
  south_north#units=m
  south_north#_FillValue=nan
  west_east#axis=X
  west_east#units=m
  west_east#_FillValue=nan
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  ( -937500.000,  937500.000) ( 97d20'18.12"W, 38d41'55.14"N)
Lower Left  ( -937500.000, -937500.000) ( 95d28'43.54"W, 21d 9'29.33"N)
Upper Right (  937500.000,  937500.000) ( 74d58'53.91"W, 38d41'55.14"N)
Lower Right (  937500.000, -937500.000) ( 76d50'28.49"W, 21d 9'29.33"N)
Center      (   0.0000000,   0.0000000) ( 86d 9'36.01"W, 30d19'48.01"N)
Band 1 Block=150x13 Type=Float32, ColorInterp=Gray
  NoData Value=9.96920996838686905e+36
  Unit Type: hPa

I want to both be able to apply a color scheme for the pixels based on the floating-point data and convert the projection to web mercator (EPSG:3857), so the goal is to use "gdalwarp" and then "gdaldem color-relief" with my color file. I seem to be missing a step as I try the following:

gdalwarp -ot "UInt16" -t_srs "EPSG:3857" tif_file.tif new_tif.tif

and receive the following error:

ERROR 1: The transformation is already "north up" or a transformation between pixel/line
and georeferenced coordinates cannot be computed for /path/to/tif_file.tif

I was suggested to add "-to SRC_METHOD=NO_GEOTRANSFORM" but this is not what I am trying to accomplish (and doing so just to see yielded the same error). When I use "gdal_translate" instead the same error comes up because it creates a virtual warped file. Adding "-a_ullr" seemed to work and the color shading code is fine, but it wants meters and not degrees. Thus putting the degree extent of the data as the corners as arguments for "a_ullr" puts the image basically at the Equator/PM intersection since -97 longitude is considered 97 meters west of the PM. Is there an intermediate step I need where this projection can be converted into meters and then using "gdal_translate" or "gdalwarp"?

Add your own answers!

Ask a Question

Get help from others!

© 2024 TransWikia.com. All rights reserved. Sites we Love: PCI Database, UKBizDB, Menu Kuliner, Sharing RPP