I'm currently working to import CANGRID climate data (provided as Surfer Grid ascii, ".grd" files) into ArcGIS. The grid is 95 rows by 125 columns in size.Metadata provides lat/lon of origin (lower left corner), cell size (50km) as well as notes projection as polar stereographic with central meridian (110 degrees W) and latitude of origin (60 degrees N).
After first attempting to convert the .grd to rasters as .ascii and .flt unsuccessfully, I've managed to use GDAL to set the extent and projection, however the dataset does not correctly align with the boundaries of the intended area. See below image.
Is there an accepted geotransformation for polar stereographic that could explain this lack of alignment?
For instance, is there a specific conversion factor or rotation that I should be using?
An example file from the dataset is here: "t201113.grd"
Here's the code I am currently using in GDAL
ds = gdal.Open("t201113.grd")
array = ds.ReadAsArray()
x_rotation = 0
y_rotation = 0
xres = 1
yres = -1
llx = -129.8530
lly = 40.0451
ulx = -175.144
uly = 71.385
input_osr = osr.SpatialReference()
input_osr.ImportFromWkt(ds.GetProjection())
wgs84_osr = osr.SpatialReference()
wgs84_osr.ImportFromEPSG(4326)
wgs_to_nps_trans = osr.CoordinateTransformation(wgs84_osr, input_osr)
x, y, z = wgs_to_nps_trans.TransformPoint(llx,lly)
geo_transform = [ x, xres, x_rotation, y, y_rotation, yres ]
ncol = ds.RasterXSize
nrow = ds.RasterYSize
out_driver = gdal.GetDriverByName("HFA")
out_ds = out_driver.Create(t201113.img", ncol, nrow, 1, gdal.GDT_Float32)
out_ds.SetGeoTransform(geo_transform)
out_prj = 'PROJCS["North_Pole_Stereographic",GEOGCS["GCS_WGS_1984",DATUM["WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Stereographic"],PARAMETER["False_Easting",0.0],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",-110.0],PARAMETER["Scale_Factor",1.0],PARAMETER["Latitude_Of_Origin",60.0],UNIT["50_Kilometers",50000.0]]'
out_ds.SetProjection(out_prj)
out_ds.GetRasterBand(1).WriteArray(array)
out_ds.GetRasterBand(1).SetNoDataValue(1.70141e+038)
out_ds.FlushCache()
out_ds = None
`
Also, here's the projection info, as defined by the input, i.e. from "GetProjection()":
'PROJCS["North_Pole_Stereographic",GEOGCS["GCS_WGS_1984",DATUM["WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Stereographic"],PARAMETER["False_Easting",0.0],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",0.0],PARAMETER["Scale_Factor",1.0],PARAMETER["Latitude_Of_Origin",90.0],UNIT["50_Kilometers",50000.0]]'
And the input GeoTransform:
(-0.5, 1.0, 0.0, 94.5, 0.0, -1.0)
Lat, longs of the grid coordinates are also provided, and when view in the projected coordinate system look like below. When the geotransform is defined by coordinates of the lower left (yellow) or upper right (pink) cordinate, I can effectively set the extent, but there remains alignment issues throughout the raster.
No comments:
Post a Comment