[gdal-dev] Projecting points on a mercator projection

[Brice Lambi]

Also, this is the gdalinfo of the image I am trying to match:

Driver: GTiff/GeoTIFF
Files: /tmp/1725.tif
Size is 1073, 689
Coordinate System is:
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["Mercator_1SP"],
    PARAMETER["central_meridian",0],
    PARAMETER["scale_factor",1],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]]]
Origin = (-14483681.357462916523218,6915896.954714655876160)
Pixel Size = (7192.681255611132656,-6738.891395318105424)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (-14483681.357, 6915896.955) (130d 6'32.84"W, 52d49'15.77"N)
Lower Left  (-14483681.357, 2272800.783) (130d 6'32.84"W, 20d 7'19.33"N)
Upper Right (-6765934.370, 6915896.955) ( 60d46'45.92"W, 52d49'15.77"N)
Lower Right (-6765934.370, 2272800.783) ( 60d46'45.92"W, 20d 7'19.33"N)
Center      (-10624807.864, 4594348.869) ( 95d26'39.38"W, 38d17'33.03"N)
Band 1 Block=1073x1 Type=Float64, ColorInterp=Gray
  NoData Value=-9999


lambi at ucar.edu wrote:
> Hello,
> 
> I'm trying to figure out how to take a list of points and put them on a
> mercator projected image using the gdal python bindings.  I'm trying to
> match it with a NDFD grid that I've been able to create and put on google
> maps from a GRIB2 file.  This is what I have:
> 
> geo = [-14483681.357462916523218,
>           7192.681255611132656,
>           0,
>           6915896.954714655876160,
>           0,
>           -6738.891395318105424]
> nx = 1073
> ny = 689
> driver = gdal.GetDriverByName("GTiff")
> dst_ds = driver.Create(fname, nx, ny, 4, gdal.GDT_Byte)
> dst_ds.SetGeoTransform(geo)
> oproj = osr.SpatialReference()
> oproj.ImportFromWkt('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["Mercator_1SP"],PARAMETER["central_meridian",0],PARAMETER["scale_factor",1],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]]]')
> 
> iproj = osr.SpatialReference()
> iproj.ImportFromWkt('GEOGCS["GCS_North_American_1983",DATUM["D_North_American_1983",SPHEROID["GRS_1980",6378137.0,298.257222101]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]]')
> 
> for i in range(1,5):
>     dst_ds.GetRasterBand(i).WriteArray(numpy.zeros((ny, nx)))
> 
> transform = osr.CoordinateTransformation(iproj,oproj)
> dst_ds.SetProjection(self.oproj.ExportToWkt())
> Rband = numpy.zeros([dst_ds.RasterYSize,dst_ds.RasterXSize])
> Aband = numpy.zeros([dst_ds.RasterYSize,dst_ds.RasterXSize])
> geom = gdal.InvGeoTransform(geo)[1]
> 
> for point in points:
>     lat = point[0]
>     lon = point[1]
>     (gx, gy, gz) = transform.TransformPoint(lon,lat)
>     (gx, gy) = gdal.ApplyGeoTransform(geom, gx, gy)
>     Rband[gy][gx] = 255
>     Aband[gy][gx] = 255
> dst_ds.GetRasterBand(1).WriteArray(Rband)
> dst_ds.GetRasterBand(4).WriteArray(Aband)
> 
> This results in a semi-correct image but it also has lots of points that
> are negative that wrap around the image and don't look right.  I also know
> this isn't taking into account that the image size is 1073x689.  What I
> expect from ApplyGeoTransform with the inverse geotransform is (x,y) where
> (x,y) are all positive and within the 1073x689 domain or can be scaled to
> fit in that domain but I don't understand the negative points.  Any
> suggestions?
> 
> Thanks,
> Brice
> 
> 
> 
> _______________________________________________
> gdal-dev mailing list
> gdal-dev at lists.osgeo.org
> http://lists.osgeo.org/mailman/listinfo/gdal-dev