Adding geo-reference data back into ArcGIS REST API exported image?

Community Help
1

Hi Community,

I am trying to figure out how to export TIFF from data served on an ArcGIS REST MapSever:
http://geoapps.icimod.org/icimodarcgis/rest/services/Nepal/Agriculture_Mask/MapServer using JavaScrpit.

I can make a request to the server that returns the TIFF using default parameters as follows:

fetch(cors_proxy + 'http://geoapps.icimod.org/icimodarcgis/rest/services/Nepal/Agriculture_Mask/MapServer/export?bbox=77.94771290110823,26.081245580749236,91.38253135336103,30.0000338114638&imageSR=4326&format=tif&transparent=true&size=&dpi=&f=pjson')

The query provides a URL link the TIFF file:

http://geoapps.icimod.org/icimodarcgis/rest/directories/arcgisoutput/Nepal/Agriculture_Mask_MapServer/_ags_mapd80853cfa46742cf8fd1adae08e613a4.tif

Despite the native projection of the data and the specified projection of this TIFF being 4326, when loading the exported TIFF in QGIS to view it, the image doesnā€™t appear to be spatially referenced.

projection-not-aligned.png
projection-not-aligned.png853Ɨ451 71.9 KB

Any insight as to what I am doing wrong or how to ensure the the export TIFF is georeferenced?

I have a sneaking suspicion that the exported TIFF is not a GeoTIFF, and that there are some additional steps necessary to append the projection data, but I havenā€™t been successful in finding a way to do this.

Here is an Observable notebook where I am playing around with different ways to load and and render GeoTIFFā€¦ with varying degrees of success. I welcome all types of feedback and guidance.

My goal is to obtain the highest resolution version of the raster data so that I can transform it into vector data for d3 maps.

Otherwise, ESRI Leaflet will render the layer using just the defaults:

image-in-leaflet.png
image-in-leaflet.png711Ɨ445 522 KB

But this is ā€˜not enoughā€™ :slight_smile:

Thank you for your time and help!

2

If you know the size of each pixel and the origin coordinates you could construct a world file:

to georeference. It is simpler to achieve as it appears from the article.

If the geoblaze or other libraries cannot accept world sidecar files, you would need to use QGIS to export an actual geotiff after loading the raster.

1 Like
3

Thanks @andreasplesch ! Iā€™ll have a go at this. Itā€™s weird, because I can plug the MapServer URL into QGISā€™ ArcGIS MapServer connector and the resulting (png32) image is spatially projected. The response from the ArcGIS server wraps a reference to the image in a JSON object that carries its geo-referencing information, but since this information isnā€™t written into the exported TIFF, Iā€™m still trying to figure out how to re-work it so that I can use it properly. Iā€™ve also tried playing around with the SVG exports and using the projection information to add it into a D3-generated map. Just canā€™t seem to get there!

4

It sounds like the ArcGIS map server just does not generate GeoTiff. Esri-Leaflet probably uses additional metadata from the mapserver to georeference the image.

looks like it is using the bounds from the request to georeference the returned image.

geoblaze parse uses georaster which lets you specify georeference information for a png:

Depending on what you want to do, that may be enough and you would not need qgis. geotiffjs appears to have a writearraybuffer function which essentially lets you generate a geotiff from a georaster:

1 Like
5

Thanks for helping me work through this. Thereā€™s definitely a solution in here! Iā€™ve learned from all of this that I donā€™t necessarily need GeoTIFF to achieve a spatially-referenced image, and that the ArcGIS REST server isnā€™t providing me a GeoTIFFā€“just a regular TIFF. Since we already know all the spatial reference information, and since one of the export options allowed by the server it SVG, it likely makes the most sense to try to just map the SVG image into our desired d3 map. I found some pointers on similar processes on Stack Overflow that maybe will help (just dropping the link here not to lose it). While I can kinda conceptually imagine how this can work, I still need to work through the materials @andreasplesch shared and also to play around more. Hopefully I figure something out over the next couple of weeks :slight_smile:

ā€¦ More soon :slight_smile: Dropping this in as inspiration to pick up the problem again later and to say that I keep plugging away at this.

6

I concur, these are the important takeaways. Corner coordinates in ā€œboundsā€ and fully defined projections suffice to georeference an image. Geotiffs just embed these data points directly in the image (header tags).

After looking more at georaster, it may be that PNG/JPG support is still pending. In that case, the qgis with world file (.pgw,jgw) approach may be easiest.

The returned SVG would still embed a raster image. Not sure if D3 can properly reproject SVG raster images (perhaps filled rectangles), as most examples deal with vector graphics. It is not quite enough to reproject just the corners and distort the embedded image accordingly. While a good approximation for small areas (up to perhaps one degree), each pixel should be reprojected individually (warping in QGIS).

1 Like
7

Thanks Andreas! I read the world file documentation. I have not yet tried manually running these calculations, as I was hoping that I could figure out how to ask QGIS to run the export. This step of manual intervention isnā€™t ideal, but Iā€™d happily enjoy a projected raster of any kind that to help me get started :slight_smile: All my attempts to export from QGIS seem to fail:

image
image1641Ɨ395 43.3 KB

Since that first notebook was a bit off conceptually, I tried narrowing down the problem here (though from your reply above I might be similarly misguided on my new SVG-focused attempt):

Iā€™ll keep trying! Thank you for sharing these insights and avenues for learning! :slight_smile:

1 Like
8

Constructing the world file boils down to these six lines in a text file:

pixel size along x, in degrees
0
0
negative pixel size along y, in degrees
x-coordinate of the center of the original imageā€™s upper left pixel, in degrees
y-coordinate of the center of the original imageā€™s upper left pixel, in degrees

for Spatial Reference: 4326 (lat/long), since the skew is usually 0.

I found that

struggles with a similar issue.

Alternatively, to export the mapserver layer from QGIS, the easiest option would be the Project - Import/Export - Export Map to Image function, essentially a georeferenced screenshot. You can define the extent and resolution which you could base on the original source.

1 Like
9

Playing with the REST API:

http://geoapps.icimod.org/icimodarcgis/rest/services/Nepal/Agriculture_Mask/MapServer/export?bbox=79.86836655581855%2C23.342939621840408%2C89.46187769865067%2C32.738305956214724&bboxSR=&layers=&layerDefs=&size=4096%2C+4096&imageSR=&historicMoment=&format=png&transparent=false&dpi=&time=&layerTimeOptions=&dynamicLayers=&gdbVersion=&mapScale=&rotation=&datumTransformations=&layerParameterValues=&mapRangeValues=&layerRangeValues=&clipping=&spatialFilter=&f=image

returns the max. 4096x4096 size image for the default extent, which is a bit clipped from the total extent.

Here is the full extent:

http://geoapps.icimod.org/icimodarcgis/rest/services/Nepal/Agriculture_Mask/MapServer/export?bbox=79.86836655581855%2C23.342939621840408%2C89.46187769865067%2C32.738305956214724&bboxSR=4326&layers=&layerDefs=&size=&imageSR=&historicMoment=&format=png&transparent=false&dpi=&time=&layerTimeOptions=&dynamicLayers=&gdbVersion=&mapScale=&rotation=&datumTransformations=&layerParameterValues=&mapRangeValues=&layerRangeValues=&clipping=&spatialFilter=&f=html

If json is requested:

http://geoapps.icimod.org/icimodarcgis/rest/services/Nepal/Agriculture_Mask/MapServer/export?bbox=77.94771290110823%2C+26.081245580749236%2C+91.38253135336103%2C+30.0000338114638&bboxSR=&layers=&layerDefs=&size=4096%2C+4096&imageSR=&historicMoment=&format=png&transparent=true&dpi=&time=&layerTimeOptions=&dynamicLayers=&gdbVersion=&mapScale=&rotation=&datumTransformations=&layerParameterValues=&mapRangeValues=&layerRangeValues=&clipping=&spatialFilter=&f=pjson

the response is:

{
 "href": "http://geoapps.icimod.org/icimodarcgis/rest/directories/arcgisoutput/Nepal/Agriculture_Mask_MapServer/_ags_map4be8f03f3ce94088ada4e9384e35359b.png",
 "width": 4096,
 "height": 4096,
 "extent": {
  "xmin": 77.947712901108247,
  "ymin": 21.323230469980118,
  "xmax": 91.382531353361045,
  "ymax": 34.758048922232916,
  "spatialReference": {
   "wkid": 4326,
   "latestWkid": 4326
  }
 },
 "scale": 1378457.1235492998
}

I think the href link may be temporary, to be consumed soon.

This is more useful since it gives you the actual extent of the returned image:

  "xmin": 77.947712901108247,
  "ymin": 21.323230469980118,
  "xmax": 91.382531353361045,
  "ymax": 34.758048922232916,

So the upper left pixel is at

xmin, ymax

and the pixel size is

(xmax - xmin)/(4096 - 1) in x
-(ymax - ymin)/(4096 - 1) in y

(I think subtracting one is correct if you consider a 2 pixel wide image pixel center to center).

So you could give this a try in a world file.

1 Like
10

Wow @andreas! Thank you! This really helps.

I am nearly there. I just misinterpreted something in your breakdown of how the World File is constructed :

image
image897Ɨ563 174 KB

I interpreted what you wrote as:

${(xmax - xmin)/(4096 - 1)} // pixel size along x, in degrees
0
0
${(ymax - ymin)/(4096 - 1)}// pixel size along y, in degrees
${xmin} // x-coordinate of the center of the original imageā€™s upper left pixel, in degrees
${ymax} // y-coordinate of the center of the original imageā€™s upper left pixel, in degrees

which in my case works out to:

0.0032807859468260803
0
0
0.0009569690429095396
77.94771290110823
30.0000338114638

Iā€™ll keep reading. Youā€™re very generous and helpful! Thank you!!

11

Clearly, something is off. x seems alright but y is flipped and the pixel size too small. Make sure not to use the input bbox but the returned extent from the json.

Ah, I noticed from the worldfile doc. and I remembered that the y pixel size should be negative. That explains the flipping at least.

I also noticed that using the full extent for bbox does not quite return the full resolution of the source data (which I do not think is available). It would be necessary to request smaller maps at 4k resolution and then stitch those together. I think this size is getting to full resolution:

http://geoapps.icimod.org/icimodarcgis/rest/services/Nepal/Agriculture_Mask/MapServer/export?bbox=77.94771290110823%2C+26.081245580749236%2C+81.38253135336103%2C+30.0000338114638&bboxSR=&layers=&layerDefs=&size=4096%2C+4096&imageSR=&historicMoment=&format=png&transparent=true&dpi=&time=&layerTimeOptions=&dynamicLayers=&gdbVersion=&mapScale=&rotation=&datumTransformations=&layerParameterValues=&mapRangeValues=&layerRangeValues=&clipping=&spatialFilter=&f=html

1 Like
12

After loading the raster into QGIS, you will need to assign manually a projection since qgis cannot know what it should be: Set CRS - Set Layer CRS (use 4326).

1 Like
13

Looking at my example with the full extent, it appears that requesting a square (4kx4k) image results in an extent which is also square in degrees (13.43 x 13.43). So the y pixel size should be the same or very similar to the x pixel size. In your case, the y pixel size would be about 3x which would work to stretch the squeezed map to about the right size, judging from the screenshot.

1 Like
14

Thank you @andreasplesch - I think you just solved it! :slight_smile:

complete_image = fetch('https://geoapps.icimod.org/icimodarcgis/rest/services/Nepal/Agriculture_Mask/MapServer/export?bbox=77.94771290110823%2C+26.081245580749236%2C+91.38253135336103%2C+30.0000338114638&bboxSR=&layers=&layerDefs=&size=4096%2C+4096&imageSR=&historicMoment=&format=png&transparent=true&dpi=&time=&layerTimeOptions=&dynamicLayers=&gdbVersion=&mapScale=&rotation=&datumTransformations=&layerParameterValues=&mapRangeValues=&layerRangeValues=&clipping=&spatialFilter=&f=pjson').then(response => response.json())

world_file=`
${(complete_image.extent.xmax - complete_image.extent.xmin)/(complete_image.width - 1)}
0
0
-${(complete_image.extent.ymax - complete_image.extent.ymin)/(complete_image.width - 1)}
${complete_image.extent.xmin}
${complete_image.extent.ymax}
`

image
image899Ɨ569 67.2 KB

I am dancing :dancer: :dancing_men: :dancing_women: :dancer: :dancing_men: :dancing_women: :dancer:

Thank you so much!!!

1 Like
15

Super. Nice use of observable to do the world file :wink:

1 Like
16

For completeness, searching for ā€˜reprojectionā€™ brings up quite a few notebooks with raster reprojections, for example Raster reprojection / Jake Low / Observable .
Perhaps one could use d3 to associate a 4326 projection with the image, probably with custom offsets and scale, now that these are known. This would effectively work to georegister, eg. retrieve geo-coordinates from image/canvas coordinates. Then one could use image processing such as contouring to retrieve georeferenced features, just like one would in gis.

1 Like