0

I downloaded the soilAtlas2M zip file, uploaded the .shp to arcgis, created a raster and clipped to my area of interest and made the map to show each soil type(map unit) to be a different colour, created a tiff and then asc file in ArcGIS Pro. When I read the asc into R the plot looks fine, but the extent is not in the format I need. (Data comes from: https://www.asris.csiro.au/themes/Atlas.html)

The asc extent/crs is in the wrong format - how can this be set in arcgis prior to exporting the asc? I have tried to fix the issue in r

I rastered the tif in r.

> s4<- raster("s4.tif")
> s4
class      : RasterLayer 
band       : 1  (of  4  bands)
dimensions : 299, 251, 75049  (nrow, ncol, ncell)
resolution : 5300, 5300  (x, y)
extent     : 14365436, 15695736, -4588206, -3003506  (xmin, xmax, ymin, ymax)
crs        : +proj=merc +a=6378137 +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null +wktext +no_defs 
source     : s4.tif 
names      : s4 
values     : 0, 252  (min, max)

The set the crs:

> s4crs<- projectRaster(s4, crs='+proj=longlat +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +no_defs')
> s4crs
class      : RasterLayer 
dimensions : 310, 261, 80910  (nrow, ncol, ncell)
resolution : 0.0476, 0.0403  (x, y)
extent     : 128.8089, 141.2325, -38.28692, -25.79392  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +ellps=GRS80 +no_defs 
source     : memory
names      : s4 
values     : 0, 255.617  (min, max)

Which affects the resolution. I can set the resolution manually to create another asc file, but lose much quality.

When I read the asc file into r:

s5<- read.asciigrid("~/s5.asc", as.image = FALSE, plot.image = TRUE)
It reads and displays:
Coordinate Reference System (CRS) arguments: NA

> s5rst<- raster(s5)
> s5rst
class      : RasterLayer 
dimensions : 299, 251, 75049  (nrow, ncol, ncell)
resolution : 5300, 5300  (x, y)
extent     : 14365436, 15695736, -4588206, -3003506  (xmin, xmax, ymin, ymax)
crs        : NA 
source     : memory
names      : s5.asc 
values     : 0, 252  (min, max)

The extent of the asc when I make it into a raster:

14365436, 15695736, -4588206, -3003506 (xmin, xmax, ymin, ymax)

The extent I need:

129.0018, 141.0028, -38.06215, -25.99639 (xmin, xmax, ymin, ymax)

Is there a way to define the extent/convert to longitude and latitude in ArcGIS Pro prior to exporting it, or a way to fix the issue in R? I believe ArcGIS has the coordinate system set to WGS_1984_Web_Mercator_auxillary_sphere. I ultimately need an asc file

2
  • 1
    The extent that it has (the big big numbers) and the extent you want to see (small numbers, lat-long degrees probably) are in different coordinate systems. You need to know what system your raster grid is in and then you can make sure your data has that set in its metadata and then you can transform if you need to. Read up on coordinate systems and projections for spatial data.
    – Spacedman
    Aug 2, 2022 at 20:45
  • Now I'm able to comment. Would you please provide the source of your file s4.tif or the file itself in order to make your example reproducible?
    – dimfalk
    Aug 3, 2022 at 19:25

3 Answers 3

1

"I ultimately need an asc file". This is your problem. ESRI "ASC" grid files do not contain any projection info. They are plain text files which look like this simple example of a 10x10 grid:

NCOLS 10
NROWS 10
XLLCENTER 1535380.000000
YLLCENTER 5083260.000000
DX 10
DY 10
NODATA_VALUE -9999
0 1 2 3 4 5 6 7 8 9
0 1 2 3 4 5 6 7 8 9
[etc]

Nowhere in that file does it say what coordinates the XLLCENTER and YLLCENTER are made from. If you try and read that file into another system you have to tell the system what coordinates they are (and in this case they are clearly not lat-long so you need to know the exact EPSG code or WKT or PROJ string).

However, if you can create a "sidecar" file with the same name as the ASC file but with a .prj extension then some software will look at this for projection info.

For example, without a .prj the above file info in gdalinfo looks like this:

Driver: AAIGrid/Arc/Info ASCII Grid
Files: tenbytenraster.asc
Size is 10, 10
Origin = (1535375.000000000000000,5083355.000000000000000)
Pixel Size = (10.000000000000000,-10.000000000000000)
Metadata:
  AREA_OR_POINT=Point
Corner Coordinates:
Upper Left  ( 1535375.000, 5083355.000) 
Lower Left  ( 1535375.000, 5083255.000) 
Upper Right ( 1535475.000, 5083355.000) 
Lower Right ( 1535475.000, 5083255.000) 
Center      ( 1535425.000, 5083305.000) 
Band 1 Block=10x1 Type=Float32, ColorInterp=Undefined
  NoData Value=-9999

But if I create a matching .prj file:

PROJCS["OSGB_1936_British_National_Grid",GEOGCS["GCS_OSGB 1936",
DATUM["D_OSGB_1936",SPHEROID["Airy_1830",6377563.396,299.3249646]],
PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]],
PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",49],
PARAMETER["central_meridian",-2],
PARAMETER["scale_factor",0.9996012717],PARAMETER["false_easting",400000],
PARAMETER["false_northing",-100000],UNIT["Meter",1]]

then now the info knows enough to compute its position in lat-long as well as its native coordinate system:

Driver: AAIGrid/Arc/Info ASCII Grid
Files: tenbytenraster.asc
       tenbytenraster.prj
Size is 10, 10
Coordinate System is:
PROJCRS["OSGB_1936_British_National_Grid",
    BASEGEOGCRS["GCS_OSGB 1936",
[etc]
                ID["EPSG",9001]]]]
Data axis to CRS axis mapping: 1,2
Origin = (1535375.000000000000000,5083355.000000000000000)
Pixel Size = (10.000000000000000,-10.000000000000000)
Metadata:
  AREA_OR_POINT=Point
Corner Coordinates:
Upper Left  ( 1535375.000, 5083355.000) (116d12'50.11"E, 78d29'59.07"N)
Lower Left  ( 1535375.000, 5083255.000) (116d12'36.14"E, 78d30' 0.60"N)
Upper Right ( 1535475.000, 5083355.000) (116d12'42.46"E, 78d29'56.29"N)
Lower Right ( 1535475.000, 5083255.000) (116d12'28.49"E, 78d29'57.81"N)
Center      ( 1535425.000, 5083305.000) (116d12'39.30"E, 78d29'58.44"N)
Band 1 Block=10x1 Type=Float32, ColorInterp=Undefined
  NoData Value=-9999

If you can avoid using ASCII grid files then do so, it will make your life easier if you can use GeoTIFFs everywhere since they contain the coordinate system info in one file with the data. If you cant and if the software creating them doesn't create matching projection files or the software reading them doesn't respect the projection files then you'll have to keep track of projections manually.

1
  • Thank you for the advice and examples. I do need to use .asc due to the modelling software I am using. However I believe I can use the s4.tif code I provided and manually set the resolution based on another tif file using the same extent and then write it as an asc. Thanks
    – pbt
    Aug 4, 2022 at 1:44
0

One of the suggestions could be: In ArcGIS, you can clip raster (present in Data management). Here is everything given in detail: https://pro.arcgis.com/en/pro-app/2.8/tool-reference/data-management/clip.htm

1
  • Your answer could be improved with additional supporting information. Please edit to add further details, such as citations or documentation, so that others can confirm that your answer is correct. You can find more information on how to write good answers in the help center.
    – Community Bot
    Aug 2, 2022 at 20:59
0
s4<- projectRaster(s4, crs='+proj=longlat +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +no_defs', res=0.04166667) #res is based on another tif file of the same extent

s4test<- mask(crop(s4, extent(mask_extent)), mask_extent) #mask_extent is the area of interest

testb<- mask(crop(saavg, extent(mask_extent)), mask_extent) #saavg is the other tif file with matching res

s<- resample(s4test, testb, method='bilinear')
#now I can create .asc of each with same geographic dimensions as is required by the modelling software

Your Answer

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

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