0

I have the following NetCDF file - I am trying to convert into raster but something is not right. The projection of the NetCDF file is not given but based on the software I received it from it should LatLong but might be cylindrical equal area. I tried both, but I keep getting this distortion which makes it impossible to query for the values at the right locations. I know the spacing of the grid is not even, not sure if that affects the end result (here visual from ArcGIS but in R it is the same problem unless plotted with levelplot function). enter image description here

library(raster)
library(ncdf4)
library(lattice)
library(RColorBrewer)

setwd("D:/Results")
climexncdf <- nc_open("ResultsSO_month.nc")

lon <- ncvar_get(climexncdf,"Longitude")
nlon <- dim(lon)
head(lon)

lat <- ncvar_get(climexncdf,"Latitude")
nlat <- dim(lat)
head(lat)

dname <- "Weekly Growth Index"

t <- ncvar_get(climexncdf,"Step")
tmp_array <- ncvar_get(climexncdf,dname)
tmp_stack <- vector("list",length(t))

for (i in 1:length(t)) {
        tmp_stack[[i]] <- tmp_array[,,i]
}

YearData <- vector("list",52)
for (i in 1:4) {
        YearData[[i]] <- tmp_array[,,i]
}   

Month1 <- YearData[c(1,2,3,4)]

# Calculate monthly averages
M1Avg <- Reduce("+",Month1)/length(Month1)

# Replace 0's with NA's
M1Avg[M1Avg==0] <- NA

# Piece of code that gives me what I need:
grid <- expand.grid(lon=lon, lat=lat)
cutpts <- seq(0,1,0.1)

# Convert to raster - work to include lat and long

M1Avg_reorder <- M1Avg[ ,order(lat) ]
M1Avg_reorder <- apply(t(M1Avg_reorder),2,rev)

M1AvgRaster <- raster(M1Avg_reorder,
                      xmn=min(lon),xmx=max(lon),
                      ymn=min(lat),ymx=max(lat),
                      crs=CRS("+proj=longlat +datum=WGS84"))
                        #crs=CRS("+proj=cea  +lat_0=0 +lon_0=0"))

r <- projectRaster(M1AvgRaster,crs=CRS("+proj=longlat +datum=WSG84"))

plot(M1AvgRaster)

# Location file not included but any locations can be entered
locations <- read.csv("Locations.csv", header=T)
coordinates(locations) <- c("y","x")

data <- extract(M1AvgRaster,locations)

writeRaster(M1AvgRaster, "M1AvgRaster_Globe_projWGSTest", format = "GTiff")
MIH
  • 1,083
  • 3
  • 14
  • 26
  • Are you sure the data file is OK? Latitude vales seem pretty random and as such it is definitely not a good file to work with. Latitude values start like: '-16.25, -20.75, -18.25, -14.25, -29.25, 51.75, 28.25, -44.25,'. After sorting they are very simple cylindrical with 0.5deg step, Longitude seems ok. – kakk11 May 10 '17 at 07:22
  • i know, that's how the software returns, it so i sorted it by M1Avg_reorder <- M1Avg[ ,order(lat) ] M1Avg_reorder <- apply(t(M1Avg_reorder),2,rev). Still however, I have issues with converting it to long lat, or any projection where i will be able to call locations by their coordinates correctly – MIH May 10 '17 at 08:41
  • Also, you may notice that step is not always 0.5 but it can be also 0.75 or 1 in some instances, in longitude and if i remember right in lat as well. – MIH May 10 '17 at 08:54
  • Ok, I see different size of grid step indeed. In python I can plot the data to map as long as coordinates are given, no matter what the projection is. Btw do you know why the data is being output in such weird coordinates? – kakk11 May 10 '17 at 09:18
  • can you query this location and see if you get some results? 25.7959° N, 80.2870° W. when i query, i get data in the ocean even though it should in theory get the location for Miami. – MIH May 10 '17 at 09:28
  • Would you mind sharing/posting the python code? – MIH May 10 '17 at 12:41

1 Answers1

1

the python version shows that after reordering at least the location of data seems correct. However, the data file seems strange, I saw data actually getting corrupted in the python netcdf library, which I've never seen before with quite a lot of different NetCDF files. Also, the chunking and compression settings are strange, better not to apply them at all.

But minimal python example to get the plot is here:

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap
from netCDF4 import Dataset

ff = Dataset('ResultsSO_month.nc')
test_var = np.copy(ff.variables['Maximum Temperature'][:])
## reorder latitudes
latindex = np.argsort(ff.variables['Latitude'][:])
## Set up map and compute map coordinates
m = Basemap(projection='cea', llcrnrlat=-90, urcrnrlat=90, 
llcrnrlon=-180, urcrnrlon=180, resolution='c')
grid_coords = np.meshgrid(ff.variables['Longitude'[:],ff.variables['Latitude'][latindex])
X,Y = m(grid_coords[0],grid_coords[1])
## Plot
m.pcolormesh(X,Y,test_var[0,latindex,:])
m.drawcoastlines()
plt.colorbar()
plt.show()

enter image description here

kakk11
  • 898
  • 8
  • 21
  • thanks so much! Is it possible to add an extra line or two to convert it to raster that will be read fine in ArcGIS and have a specific projection? – MIH May 11 '17 at 08:15
  • You wan't to end up with GeoTIFF in lonlat projection and WSG84, right? As much as I know, you would actually need consistent grid steps for that, GTiff does not save coordinates separately AFAIK, but only the trasformation. So you need to throw away some lines of data or interpolate something in between, I guess. – kakk11 May 11 '17 at 08:51
  • Yeah, tiff or img with a projection - any really would do, equal area cylindrical for example. But it doesnt seem like an option without adding extra columns right? i think thats what i will need to do. what does AFAIK stand for? – MIH May 11 '17 at 09:46
  • AFAIK -- "as far as I know". And yes, I believe you need to add extra columns, if correcting the initial data processing is not an option. If no actual science is being done in ArcGIS, then simple interpolation might work. – kakk11 May 11 '17 at 09:59