[GIS] Reprojecting from Lambert Conformal Conic into EPSG:3857

cs2csgdalgdal-translategdalwarpweb-mercator

I'm trying to reproject a plain png image file into EPSG:3857 for use in an OpenLayers application. I'm pulling the image from the NWS site at: http://rapidrefresh.noaa.gov/HRRR/for_web/hrrr_ncep_jet/2015032517/full/cref_sfc_f00.png

According to that site the map parameters are:

    :CEN_LAT = 38.5f ;
    :CEN_LON = -97.5f ;
    :TRUELAT1 = 38.5f ;
    :TRUELAT2 = 38.5f ;
    :MOAD_CEN_LAT = 38.5f ;
    :STAND_LON = -97.5f ;
    :POLE_LAT = 90.f ;
    :POLE_LON = 0.f ;

Below are corner points for WRF-HRRR mass points, SW-NW-NE-SE

    :corner_lats = 21.13812f, 47.84364f, 47.84364f, 21.13812f;
    :corner_lons = -122.7195f, -134.0986f, -60.90137f, -72.28046f;

I've referenced Making geotiff and reprojecting to EPSG:4326 using GDAL? . I've cropped the frame down to just the mapped area and I've tried running cs2cs to get the map units I need to run gdal_translate and gdalwarp but I can't seem to get the correct (WebMercator) image out of gdalwarp. I think my problem is that I don't have the projections quite right in cs2cs. I'm running:

 cs2cs +init=epsg:3857 +to +proj=lcc +lat_1=38.5 +lat_2=38.5 +lat_0=38.5 +lon_0=-97.5 +k_0=1 +a=6371200 +b=6371200 +units=m +no_defs
 -134.0986 47.84364
 10996563.15    1837306.24 0.00
 -72.28046 21.13812
 10996629.05    1837356.50 0.00

Then 

 gdal_translate -a_srs "+proj=lcc +lat_1=38.5 +lat_2=38.5 +lat_0=38.5 +lon_0=-97.5 +x_0=0 +y_0=0 +k=1.0 +a=6378137 +b=6378137 +units=m +no_defs" -a_ullr 10996563.15    1837306.24 10996629.05  1837356.50 cropped2.tif cropped2_translated6.tif

Then finally: 

 gdalwarp -t_srs "+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=-90.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs" cropped2_translated6.tif cropped2_warped.tif

My final warped image comes out with N at lower left S at Upper Right E at Upper left and W and lower right so the US looks like it is falling over backwards.

A gdalinfo of this warped file shows:

 gdalinfo cropped2_warped.tif
 Driver: GTiff/GeoTIFF
 Files: cropped2_warped.tif
 Size is 850, 917
 Coordinate System is:
 PROJCS["unnamed",
  GEOGCS["unnamed ellipse",
    DATUM["unknown",
        SPHEROID["unnamed",6378137,0]],
    PRIMEM["Greenwich",0],
    UNIT["degree",0.0174532925199433]],
 PROJECTION["Mercator_1SP"],
 PARAMETER["central_meridian",-90],
 PARAMETER["scale_factor",1],
 PARAMETER["false_easting",0],
 PARAMETER["false_northing",0],
 UNIT["metre",1,
    AUTHORITY["EPSG","9001"]]]
 Origin = (10012443.303476206958294,7747.939850974368710)
 Pixel Size = (0.072684385553582,-0.072684385553582)
 Metadata:
   AREA_OR_POINT=Area
 Image Structure Metadata:
   INTERLEAVE=PIXEL
 Corner Coordinates:
 Upper Left  (10012443.303,    7747.940) (  0d 3'24.09"W,  0d 4'10.56"N)
 Lower Left  (10012443.303,    7681.288) (  0d 3'24.09"W,  0d 4' 8.41"N)
 Upper Right (10012505.085,    7747.940) (  0d 3'22.09"W,  0d 4'10.56"N)
 Lower Right (10012505.085,    7681.288) (  0d 3'22.09"W,  0d 4' 8.41"N)
 Center      (10012474.194,    7714.614) (  0d 3'23.09"W,  0d 4' 9.49"N)
 Band 1 Block=850x3 Type=Byte, ColorInterp=Red
 Band 2 Block=850x3 Type=Byte, ColorInterp=Green
 Band 3 Block=850x3 Type=Byte, ColorInterp=Blue

What am I doing wrong?

Best Answer

Your input values of your first cs2cs statement are latitude/longitude (decimal degrees), not EPSG:3857. And they're probably relative to the sphere that you're referencing, not WGS84 (EPSG:4326).

You first need to convert the latitude/longitude values to the Lambert conformal conic (LCC) coordinate system. Label them as lat/lon in the cs2cs statement. Right now you're saying they're in 3857 which is a projected coordinate system. It should start with "+proj=latlong +a=6371200 +b=6371200...". with the same sphere radius as the LCC system.

Related Question