2

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?

1

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.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.