[GIS] Geotransformation for polar stereographic

arcgis-desktopcoordinate systemgdalpython

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.

enter image description here

Best Answer

Too long for a comment, this is to accompany @Matej's answer.

  1. Add the “.grd” data into ArcGIS.

  2. Use the “Raster to Other Format” function and convert your .grd file into an ESRII GRID format. This is important because most of the raster functions in ArcGIS are accessible for this format only, either that or it usually gets too slow when you use it on other formats.

  3. As it already has the projection file associated to the file. Rather than projecting the new converted data, define its projection. ArcToolbox > Data Management Tools > Projections and Transformations > Define Projection. You can navigate to the pre-defined polar stereo graphic ESRII projection and see if its parameters match the one given in the metadata (it doesn’t), so you can modify it as per @Matej. Only here - rather than modifying, create a new one based on the NPS projection with the central meridian and Latitude of origin changed and save it on disk, then navigate to the new projection and use that when defining your projection. This is because your on the fly modification won't be available later on when you want to use it to set the coordinate system for your data frame, which is necessary in order to align your other boundary files to the same modified NPS system.

Related Question