MassGIS Dataset: I'm using the MassGIS Data: USGS Color Ortho Imagery (2013/2014) aerial images of Boston. Basically, the dataset is divided into small patches of jpeg2000 images, where each of the images have a corresponding xml file. One important piece of data in the xml file is the bound of that image in latitude and longitude. They are using EPSG:26986 while reporting this. The images are ortho-rectified.
OSM Dataset: I donwloaded just the building polygons of Boston using the Overpass Turbo API. As far as I understand, they are denoted using WGS84 datum and the corresponding EPSG is 4326.
What I want to do: For some image from the MassGIS dataset, I want to draw the building polygons on top of it using python.
What I tried: For the given image, I retrieve the bounding latitudes and longitudes of the image from the corresponding xml file. EPSG:26986 and EPSG:4326 seems similar enough, so I don't convert one to another. Then for all of OSM building polygons that lie inside the MassGIS image bound, I find the corresponding pixel position of each node of the building polygons using a linear transformation-
x = (longitude - longitude_bound_min) / (longitude_bound_max - longitude_bound_min) * image_pixel_width y = (latitude - latitude_bound_min) / (latitude_bound_min - latitude_bound_max) * image_pixel_height
Doing the above, I am getting a mis-alignment between MassGIS images and OSM building polygons, as you can find here. My understanding was that since the images are ortho-rectified, this transformation should be enough. I would be really glad if anyone gave me any hint about what I am misunderstanding here.
Edit: For one of the MassGIS aerial images, I get the following denoted as a bound for the image in the xml file-
<bounding> <westbc>-71.240667</westbc> <eastbc>-71.222039</eastbc> <northbc>+42.196636</northbc> <southbc>+42.182783</southbc> </bounding>
There is also a corresponding .aux file for image, running `gdalinfo' on which gives me-
Driver: HFA/Erdas Imagine Images (.img) Files: 19TCG150725.aux Size is 5000, 5000 Coordinate System is: PROJCS["NAD_1983_UTM_Zone_19N", GEOGCS["GCS_North_American_1983", DATUM["North_American_Datum_1983", SPHEROID["GRS_1980",6378137,298.257222101]], PRIMEM["Greenwich",0], UNIT["Degree",0.017453292519943295]], PROJECTION["Transverse_Mercator"], PARAMETER["latitude_of_origin",0], PARAMETER["central_meridian",-69], PARAMETER["scale_factor",0.9996], PARAMETER["false_easting",500000], PARAMETER["false_northing",0], UNIT["Meter",1]] Origin = (315000.000000000000000,4674000.000000000000000) Pixel Size = (0.300000000000000,-0.300000000000000) Corner Coordinates: Upper Left ( 315000.000, 4674000.000) ( 71d14'26.40"W, 42d11'46.62"N) Lower Left ( 315000.000, 4672500.000) ( 71d14'24.68"W, 42d10'58.02"N) Upper Right ( 316500.000, 4674000.000) ( 71d13'21.04"W, 42d11'47.89"N) Lower Right ( 316500.000, 4672500.000) ( 71d13'19.34"W, 42d10'59.29"N) Center ( 315750.000, 4673250.000) ( 71d13'52.87"W, 42d11'22.96"N) Band 1 Block=64x64 Type=Byte, ColorInterp=Undefined Description = Layer_1 Metadata: LAYER_TYPE=athematic Band 2 Block=64x64 Type=Byte, ColorInterp=Undefined Description = Layer_2 Metadata: LAYER_TYPE=athematic Band 3 Block=64x64 Type=Byte, ColorInterp=Undefined Description = Layer_3 Metadata: LAYER_TYPE=athematic Band 4 Block=64x64 Type=Byte, ColorInterp=Undefined Description = Layer_4 Metadata: LAYER_TYPE=athematic