0

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
Improve this question
3
  • 2
    I think you had better to use ogr2ogr to change coordinates of the buildings from EPSG:4326 into EPSG:26986 or use gdal_translate for the rasters.
    – Zoltan
    Commented Feb 28, 2018 at 6:49
  • They're not off by much (relatively), so the misalignment looks like a NAD83 <> WGS84 problem.
    – mkennedy
    Commented Feb 28, 2018 at 17:49
  • My bad, turns out the bounds were probably not given in NAD83 datum, but instead WGS83 datum. Do you have any idea why this misalignment might occur then?
    – crysoberil
    Commented Mar 1, 2018 at 0:33

1 Answer 1

Reset to default
1

I think your assumption that EPSG:26986 and EPSG:4326 are "similar enough" is wrong. One is in metres and the other is in degrees.

You should use either gdalwarp to transform the images into EPSG:4326 or ogr2ogr to transform the OSM vectors into EPSG:26986.

Improve this answer
3
  • Interesting point. The bounds are represented in this way- <westbc>-71.240667</westbc> <eastbc>-71.222039</eastbc> <northbc>+42.196636</northbc> <southbc>+42.182783</southbc> Does that mean they are already represented in the WGS84 degrees?
    – crysoberil
    Commented Feb 28, 2018 at 23:50
  • depends on how you got the bounds - please add the gdalinfo result for one of them to your question
    – Ian Turton
    Commented Mar 1, 2018 at 10:06
  • I did add that to the question. I suppose, in the xml they already put in the WGS84 coordinates as bounds(they don't mention it anywhere in their dataset documentation though). But in that case, I still don't understnad why I'm getting the misalignment.
    – crysoberil
    Commented Mar 1, 2018 at 12:19

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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