0

I imported an image from a Arc GIS REST Map service and the dimensions of the image as projected don't match the transform, but they do line up correctly in Google Earth with the expected building. The measured width in Google Earth is about 82.5 meters, but the expected width is 100 meters. I am trying to calculate the area for each color, so this discrepancy is messing up my area calculation.

The well known ID of the service is EPSG 3857.

I am new to working with raster data, so I know there is a lot that could be improved here.

import rasterio
import urllib
from requests import get

Here is the call to the API to get the image.

url_base=r'https://public.gis.lacounty.gov/arcgis/rest/services/Solar/Solar_Point_Cache/MapServer/export?'
bbox = '-13165211.089884153,4037019.7962294333,-13165111.089884153,4037119.7962294333'
params={'bbox':bbox, 'bboxSR':'3857','imageSR':'3857','format':'png24',
       'transparent':'true','f':'json'}
url = url_base+urllib.parse.urlencode(params)
response=get(url)

Save the response to a file. However, the data has no spatial reference in the image itself, and the image isn't returned by the API response. The API response includes a link to retrieve the image. The spatial reference data is all in the json response.

data = json.loads(response.text)
df_temp = pd.json_normalize(data)
img=Image.open(BytesIO(get(df_temp.href.loc[0]).content))
img.save(facility+'.png')

This is the value of the json response from the requests.get call

{'href': 'https://public.gis.lacounty.gov/arcgis/rest/directories/arcgisoutput/Solar/Solar_Point_Cache_MapServer/_ags_mapccd98bd8c5694ef2aa4e0875b5ff1865.png',
 'width': 400,
 'height': 400,
 'extent': {'xmin': -13165211.089884153,
  'ymin': 4037019.7962294333,
  'xmax': -13165111.089884153,
  'ymax': 4037119.7962294333,
  'spatialReference': {'wkid': 102100, 'latestWkid': 3857}},
 'scale': 944.8799999999998}

Define the transform and the crs

transform = rasterio.transform.from_bounds(df_temp['extent.xmin'].loc[0], df_temp['extent.ymin'].loc[0],\
                                           df_temp['extent.xmax'].loc[0],
                                           df_temp['extent.ymax'].loc[0], dataimg.shape[1], dataimg.shape[2])
crs = {'init': 'epsg:'+str(df_temp['extent.spatialReference.latestWkid'].loc[0])}

So here is the transform and the CRS. Based on this I think that each pixel should have a width of 0.25 meters in the x and y directions. However, the image as rendered in Google Earth has a width of about 0.206 meters in the x and y. Why are they different?

Affine(0.25, 0.0, -13165211.089884153,
       0.0, -0.25, 4037119.7962294333)
{'init': 'epsg:3857'}
2
  • The image from ArcGIS server is in EPSG:3857 and a meter in EPSG:3857 corresponds to one meter on ground only on top of the equator. See the indicatrix at openlayers.org/en/latest/examples/tissot.html. In Google Earth the scale is about correct everywhere.
    – user30184
    May 31 at 18:50
  • 1
    Thanks. I changed my projection to EPSG 2874, which is specifically for the US, and it worked like a charm.
    – Foggy
    May 31 at 20:46

1 Answer 1

2

This answer is based on @user30184's feedback.

The problem was that the projection I was using, 3857, which was the one used by the REST API, has a lot of distortion as I moved away from the equator. This is visible here, https://openlayers.org/en/latest/examples/tissot.html, which was suggested by @user30184.

In projection 3857 my image was 100 meters wide, but when it was projected by Google Earth at my location it was only 82 meters wide or so.

I changed my bbox projection to one specifically for the US, 2874, and changed the appropriate parameters as well. Now when I plot the image in Google Earth and measure it, the measurement is 399 feet. Case closed.

New parameters and bounding box

bbox='6481312.353165579,1844685.4502506657,6481712.353165579,1845085.4502506657'
params={'bbox':bbox, 'bboxSR':'2874','imageSR':'2874','format':'png24',
       'transparent':'true','f':'json'}

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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