texts: Creating raster with same size as DTM raster containing certain percentage of zeros and ones using QGIS

texts:
I've got a DTM and would like to add 1m to a percentage of the points, so I thought of creating a new raster containing either zeros or 1m and add the two raster together.
I didn't find any tools that could do it on SAGA or GRASS.
Do you know a way?


texts: Crop a raster by a grid shapefile with different extensions in R

texts:
I need crop a raster by a grid shapefile to obtain raster files cutted according to grid cells (which are named "PageNumber"). My function is:

  g=shapefile("50x50km_final.shp")

  f=raster("declivcerrado.tif")

  cell=nrow(g)

  i=1

  for (i in 1:cell){
    print(paste("Processando paisagem",i,"de",cell))
    tmp=subset(g,g$PageNumber==i)
    tmp1=crop(f,tmp)
    writeRaster(tmp1,paste("paisag",i,sep=""),format='HFA',overwrite=T)
  }

But I got this error:

  Error in .local(x, y, ...) : extents do not overlap

The files extent and projection are:

  extent(g)
  class       : Extent
  xmin        : -654873.4
  xmax        : 1381884
  ymin        : -1416856
  ymax        : 1056824

  extent(f)
  class       : Extent
  xmin        : -60.40702
  xmax        : -40.48814
  ymin        : -24.70977
  ymax        : -2.266753

  projection(g)
  [1] "+proj=aea +lat_1=-2 +lat_2=-22 +lat_0=-12 +lon_0=-54 +x_0=0 +y_0=0 +ellps=GRS80 +units=m +no_defs"

  projection(f)
  [1] "+proj=longlat +ellps=GRS80 +no_defs

How can I solve this problem?


texts: GeoKettle, adding Source Values into Value Mapper

texts:
I'm testing GeoKettle to see if it can help me transfer shapefiles from one schema to another.  I'm trying to use the value mapper but am finding that all Source Values and corresponding Target Values have to be entered manually - is there a way to make this process more efficient.

I know that in FME you can scan a feature to automatically read in
source values - this doesn't seem to be possible in GeoKettle.
I also understand that FME allows you to somehow attach an external
lookup table - I can't see a way to do this in GeoKettle.
I have tried copying values from a spreadsheet or text document but
they do not paste into the ValueMapper dialog in the correct way.


texts: What is the most efficient method to delete points?

texts:
I have a script that takes a point featureclass (FC) through the following workflow:

Select points that lie within individual polygons
From the selected points within each polygon, randomly subset X% of
points
Merge the selected points from each polygon feature sub-selection
Erase the original points with the merged subset and write to a new FC

I am using arcpy.Erase_analysis() to do the point deletion.

The subset of data I am testing this script on has ~35000 point features.  The point subset that I would like to use as the erase feature has ~3000 point features.  It takes approximately 3-5 minutes to complete processing.  Seeing that I will be scaling this up to a dataset containing ~30 million points, is there a more efficient method to delete points than by using the arcpy.Erase_analysis() method?  Perhaps, programmatically deleting points in an editing session?  Or, a method that does not involve writing a new point FC?


texts: Mapinfo Line and polygons from different layers to be snaped

texts:
I have a layer that I need to edit [ kerb and gutter (KNG)] which consists only of lines that represent the kerb and gutters on both sides of the road. I have separated the left and right side of the KNG and label them accordingly.

The KNG naturally should be sitting inside roads. We have a layer called road section that contains polygons that traces the road.

I want to make each line of the KNG layer to be completely within their respective road sections (think the lines are continuous and connected but is snapped at the end of each road section) so that when I run a SQL with completely within, it should all show up.

There seems to be an issue with the bounds.

When I save the two tables without bounds and then pack the table, (I realised that packing the table after saving without bounds changes the bounds slightly) I am able to achieve what I want albeit needing to do a little bit of manual editing.

This, however, would not really work as I want to be able to get the KNG table into the corporate layer so that everyone can use it.

When I change the bound of the KNG so that it matches the bound of the road section layer, no matter what I do, the nodes on the KNG and road section would not snap.

In the without bounds version, I did:

Overlay node onto road section, then overlay nodes onto KNG lines. This somehow combines the two nodes. I can then run some SQL or do some manual edit to get the result I want.

The unbounded tables seem to just be allocated to a default bound, which is unnecessarily big for us.

In the bounded version:
Did the same as above, I can see that the two nodes are off by a few mm not a bit deal for us in the grand scheme of things but some SQL would not work.

Is there any function in MapInfo to allow me to make the end of a line snap with a node on a different layer within a threshold say 50cm or so?


texts: Dynamically cutting shapefiles based on an eliptical or polygonal form

texts:
I have been trying to cut a shapefile in small parts in order to reduce the processing of the graph created based on that map.
I did find the ogr2ogr clipsrc option however, as far as i know, that one uses only a bounding box which is ok, but still includes unnecessary data that i want to get rid off.

Basically i need to cut  the shapefile in an elliptical form and that ellipse would have the major axis determined by the 2 points that i would insert. Is there any way i can accomplish that?


texts: How to create dissolved buffers in QGIS with PostgreSQL/PostGIS?

texts:
Creating a dynamic buffer in QGIS I use the following query:

CREATE VIEW buffer40units AS
SELECT gid, st_buffer(geom,40)::geometry(polygon, 31492) AS geom
FROM point;

Now I want to dissolve overlapping buffers but all queries I've tried aren't working.

CREATE VIEW buffer40units AS
SELECT gid, st_buffer(st_collect(geom),40)::geometry(polygon, 31492) AS geom
FROM point;

CREATE VIEW buffer40units AS
SELECT st_buffer(st_collect(geom),40)::geometry(polygon, 31492) AS geom
FROM point;

It seems that QGIS needs a (new) gid for every row in the view. Can anybody please help me?


texts: Sieve function in QGIS not working properly?

texts:
I am trying to eliminate single pixels in a (binary) raster with the gdal_sieve function. However, regardless of the settings I choose (I tried threshold=1, 2 and 70 and pixel connection 8), the result is identical to the original raster.
Any insight on this issue?


texts: Web Mercator to 3d globe projection problems

texts:
I've been trying to use map tile images from services like OpenStreetMap to create a globe, and I've mostly succeeded. The problem is that when I compare my globe to one using the same image source made with Cesium.js, mine looks weird and stretched.

Here's mine:

Here's the one that Cesium.js makes

As you can see, the one that I've made fits all the images together properly, but they seem slightly too large.

What could be causing this?

My methodology for texturing the globe is as follows

Get the latitude and longitude of the vertex
Project lat/lng to a x/y value in tile-space using leaflet.js
Use that x,y coordinate to sample the image and put it on my sphere

Those three steps seem pretty straightforward to me, but I'm clearly doing something wrong and I'm not sure what it is.


texts: How to union polygons with common ID in QGIS?

texts:
I have a shapefile with lots of polygons and I need to union all polygons which have the same identifier-id into a new (multi-)polygon, ideally with the same identifier-id on another (new) layer.

As there are lots of polygons I don't want to do this manually. Is there a simple way I can do this automatically in QGIS?  I suppose This can be done with Spatialite, but that's only a guess.

Please help a rookie. Any hints appreciated. Thank you!


texts: How do I combine geometries in a shapefile based on a grouping variable?

texts:
Here is a simple example. I have a shapefile of three US states:

This is the code that creates the shapefile.

library(USAboundaries)
library(sf)
library(dplyr)

states <- us_states(map_date = "2000-01-01", resolution = "high", states = c("CA", "OR", "WA")) %>%
  mutate(group_var = if_else(state_abbr == "CA", 1, 2))
plot(st_geometry(states))

Now, I'm trying to combine the geometries based on the value of group_var, so California should stand alone while Oregon and Washington get lumped together into a single geometry. Unfortunately, st_combine doesn't take a grouping variable, and although the aggregate function in sf looks promising, code like this throws an error that group_var cannot be found.

x <- aggregate(states, group_var, mean)

Furthermore, aggregate requires a function as a third argument, presumably because aggregate in the stats package applies the function to the data, but in this case, there isn't any data to apply a function to. I'm just trying to combine/aggregate the shapefiles.


texts: Where to look for oceanographical/environmental data for the Arctic (Baffin Bay/Davis Strait area)

texts:
I am looking for oceanographical/environmental variables that are available for Baffin Bay/Davis Strait area over the time period of 1995-2011, however, I am not familiar with oceanographical/environmental datasets.

Perhaps some of you may be able to point me in the right direction.
Which oceanographical/environmental variables/layers are available from this part of the world and for which time period? Where can I access these?

Thank you very much for your help.


texts: Seeking library that reads/writes shapefiles to use with .net core 2.2?

texts: Previously I used to use Easy Gis .net. Now we have a .netcore 2.2 project and I need to read/write Shapefiles. I have googled for 2 hours and I didn't find anything. I also tried to use dotspatial and easy gis on .netcore 2.2 and it failed.


texts: arcpy CSV file to raster attribute table?

texts:
Using TableToTable_conversion I can convert a CSV file to a table in a geodatabase, but how can make this table be the default attribute table for a raster in the geodatabase?

I can't seem to find anything on this...


texts: Adding circular string feature using PyQGIS?

texts:
I am attempting to add a circular string feature (in this case an arc of a circle).  The following code runs without error but doesn't create a feature.

I am using version 3.4.2.

vl = QgsVectorLayer("Polygon", "temp", "memory")
from qgis.PyQt.QtCore import QVariant
pr = vl.dataProvider()
pr.addAttributes([QgsField("name", QVariant.String)])
vl.updateFields()
cString=QgsCircularString()
cString.setPoints([QgsPoint(306580,167317),QgsPoint(306680,167317),QgsPoint(306680,167817)])
print (cString.hasCurvedSegments())
geom_cString=QgsGeometry(cString)
f = QgsFeature()
f.setGeometry(geom_cString)
f.setAttributes(["One"])
pr.addFeature(f)
vl.updateExtents()
QgsProject.instance().addMapLayer(vl)


texts: QGIS Map Tips Table, style every other row

texts:
I simply want to style every other row in a table I've generated using map tips with a different coloured background. I've got the correct generic CSS (I think) but nothing is changing.

tr:nth-child (odd) {
background:grey!important;
}

Below is the CSS styling for the table, I feel the contents of the table is irrelevant to show as code but I have shown an image below.

<style>

body {
background-color: white!important;
font-family: sans-serif;
width: 200px;
}

table {
border-collapse: collapse;
}

th {
font-size: 20px;
text-transform: uppercase;
padding: 10px;
background: grey;
color: white;
}

tr {
border: 2px solid white;
}

tr:nth-child (odd) {
background:grey!important;
}

td {
font-size: 12px;
white-space: nowrap;
padding: 5px;
}

td.h2{
font-size: 15px;
float: left!important;
font-weight: bold;
background:grey;
}

td.bold {
font-size: 12px;
font-weight: bold;
}

td.gap {
background-color:white;
padding:1px;
}

</style>


texts: Combining two TileWMS with different axisOrientation using OpenLayers3?

texts:
I have two TileWMS, one that has map data, with terrain, roads and stuff. And one that only is images that should be placed upon the map. Both uses the same projection, ESPG:3006.

The one hosting the map is a WMS of version 1.1.1 and the other 1.3.

Is it possible to set these to to work?

I have managed to setup a map where both have the same axisOrientation, specifically "neu" but I have not managed to get it to run with different axisOrientations.


texts: Choosing Java or .net for GIS application

texts: I currently use python for GIS scripting in arcgis and image processing using GDAL library. Now, I want to improve my programming skills for   GIS application development both for desktop as well as internet applications using either Java or .net. I am a bit confused at this moment which will be better for learning as well as better scope for industry jobs.


texts: QGIS function layer_name in fieldcalculator batch-mode

texts:
Working with QGIS 3.10, I want to use the @layer_name function in batch mode for all layers in a project. In individual cases, I do this via the field calculator.
I also found the script for the field calculator in batch mode in the toolbox
but somehow I don't get that configured.
Can someone explain how to do this?
What I mean is that fieldcalculation in batch-mode
for all layers:


texts: Create an Identify Route Locations tool that converts M values

texts: I have a data set with many maritime routes that I would like to see the m values (time) for. ArcGIS's  Identify Route Locations tool works exactly as I want it to, however, the m value displayed is in Unix time which is not easily read, if at all. I need to convert the Unix times to something readable. I found that using datetime.datetime.fromtimestamp().isoformat() works great but I don't want the user to have to convert the time by hand every time they click along the line. Is there a way to create a tool that works like Identify Route Locations but also converts the m values before displaying them? I can't figure out how to access them. I've looked around and the only solutions I've found were to break down the routs back into points and recalculate the m values, but I'd rather avoid that.


texts:
I'm working on site accessibility. Having various locations, I can identify the LCP going from one to any other. So, by default, LCP calculation is giving me the first most least-cost path only. I can use various programs to get a few different options, but they generally point to the same general direction.

I wonder if anyone knows a way to get other paths which, despite being more costly, could be used as an alternative. As I see it, it could imply using a simple XY chart (X being each alternative and Y the cost) and splitting it in natural breaks or quartile, or whatever method. I suppose it would need to calculate a lot more.

Maybe it already exists in some road GPS applications (for an example, pointing out various roads to avoid jams or something) ? If so, how could it be used with elevation rasters?


texts: Reclassifying NLCD data, specifically open water

texts:
I've been working with NLCD data at a county level. My goal is to reclassify the pixels that fall under the open water classification. I've got a hydrography shapefile for the same geographic location so I can compare and areas and see what the open water actually classifies as, whether it's river, lake etc. Is there an easy tool to reclassify these based on the hydrography. For some reason I'm drawing a blank on the most efficient way to do this? Pretty much I want to take the open water cells from my raster and change them from open water to river lake or whatever based on what the hydrography attribute table says.

NLCD : National Land Cover Database. It's published by the USGS every five years in raster format 30x30m pixels. Each pixel is classified by a number which corresponds to a legend representing ladncover type. For instance, any pixel that is classified as 11 means open water. So basically I want to classify all the 11 pixels more in depth to indicate whether they are a lake, river or actually open water since I am looking at an area around the puget sound


texts: How to create a raster from NoData values using GDAL?

texts:
I have 1,000 very large single-band rasters and many of them have nodata values. I want to take those nodata values and create a binary map to visualize all of the areas where I do not have data coverage.
Even after researching for quite awhile I don't really understand alpha bands/masks, where they exist, how to examine them. I thought I could just add an alpha band in GDAL but when I export the band I just get a raster full of 255s (checking with np.unique(tif))
# Adds alpha band
gdalbuildvrt out_vrt src_dir/*.tif -b 1 -vrtnodata -9999 -addalpha

# Deletes the first band, leaving just the alpha
gdal_translate in_vrt out_tif -b 2

My next idea was to use gdal_calc to create a binary tif where nodata=1 and data=0
# Both of these attempts left me with just a raster of zeroes
# Attempt #1
gdal_calc.py -A src_tif.tif --A_band=1 --calc="(A<=-9999)" --outfile out_tif.tif --co compress=Deflate --type=Byte

# Attempt #2
gdal_calc.py -A src_tif.tif --A_band=1 --calc="(A*0)" --NoDataValue=1 --outfile out_tif.tif --co compress=Deflate --type=Byte

Any tips? This seems like it should be so simple, it's driving me crazy


texts:
I have a series of JPEG images taken from a UAV that have geospatial info associated with them. I'm trying to pull out that info in python but I'm having some trouble. In the past when I've done this the images are either GeoTIFF files that have all the necessary info in the header, or they had an associated world file (.tfw or .jpw) that contained the location of the top left pixel, along with a rotation or heading so you can what direction the image is facing. The image I'm using right now doesn't have a world file. I'm able to use GDAL to pull out the lat/long of the upper left pixel, but north isn't up in this image and I can't find that information in the image.

So far I've mostly used gdal:

src = gdal.Open(img_path)
src.GetMetadata_Dict()

With gdal I can easily get the corner location:

Which I can easily parse to degrees instead of minutes/seconds. But I'm not seeing any rotation or heading information. Similarly with PIL I can pull out the corner lat/long, but I'm not able to find what direction is up in the image. But the information must be in there because when I throw it into a GIS the image appears correctly geolocated and rotated.

Does anybody know how I can pull that information using Python?


texts: Get Transformation method and Parameters from coordinate system

texts:
I have data (vectors, and raster) that are in UTM coordinate system, and I need to reproject them to a local coordinate system. Unfortunately, I do not know the transformation method that I may use either the transformation parameters (there or seven parameters) :(
However, I can get get some points surrounding the AOI in both UTM and equivalent local coordinate system.
My question is, how can I get the transformation parameters using these points and the UTM and local projection parameters to get the transformation method and parameters?


texts: Simple Folder to folder batch python script

texts:
I am writing script to take shape files from a folder to Buffer it to given distance to Output folder with same name as input files. Here are block of codes.

# Name: Buffer_multiple file.py
# Usage : inFolder to outfolder Buffer analysis

# import system modules
import arcpy
from arcpy import env

# Set environment settings

arcpy.env.workspace= arcpy.GetParameterAsText(0)
arcpy.env.overwriteOutput = True

newFC = arcpy.ListFeatureClasses()

# get buffer size
buffer = arcpy.GetParameterAsText(1)

# get output feature class
outputFeatureClass = arcpy.GetParameterAsText(2)

for fc in newFC:
  arcpy.Buffer_analysis(fc, outputFeatureClass+fc, buffer)

It shows succeeded with No Output shapefile.

Executing: Script C:\Users\user1\Desktop\Input 20 C:\Users\user1\Desktop\buffered
Start Time: Sun May 17 12:39:12 2015
Running script Script...
Completed script Script...
Succeeded at Sun May 17 12:39:14 2015 (Elapsed Time: 2.00 seconds)


texts: Identifying rasters with specific attribute value that intersect polygon using ArcGIS for Desktop?

texts:
I have numerous raster files that each contain an attribute field "VALUE" that can take on the values 0 or 1.

How can I identify rasters that have any of their "VALUE=1" cells intersecting a focal polygon (this same polygon will be used for all rasters)?

I included a picture in case it helps. My bounding, spatial extent polygon is shown in tan and a single raster is displayed alongside it, whereby black ("0") cells do not contain the focal fish species of interest and blue ("1") cells do contain the focal species. I only care about raster-polygon intersection for "VALUE=1" cells in each raster (i.e., the fish species occurs within my focal study extent). I'm willing to convert all rasters to shapefiles if that leads to a more efficient tool.


texts: ArcMap python window refusing to execfile after first attempt, requires restart

texts:
I am attempting to execute a script mulitple times from the ArcMap python window (using different parameters that are hardcoded into the file).

The script runs the first time successfully, however it seems to "break" the python console, and other attempts at execfile are unsuccessful afterwards, until I restart arcgis.

Here's the a basic scenario:

Start ArcGIS, open python console, type:

execfile(r"C:\Users\xxx\Desktop\test.py")

Script executes successfully
Modify test.py script
Type execfile(r"C:\Users\xxx\Desktop\test.py")
Python console returns to prompt

Any ideas on what might be happening? Workarounds?


texts: Editing Specific Raster Cells Via CSV Values using GDAL?

texts:
I have a raster grid that I would like to modify. I am trying to manually edit specific raster cells' values via vales in a CSV file.

My CSV contains three columns:

VAL,LAT,LON
91.0,31.04250,-86.31167
91.2,30.21889,-96.37417
75.9,41.00928,-74.73628
86.1,36.61056,-83.73722
89.0,34.28028,-87.60028
64.9,43.77714,-71.75598
68.0,44.42010,-72.01930
37.0,45.60444,-103.54639
68.0,39.16052,-89.67419

Task Example

In the above example, I would try to change the middle cell value from 65 to 67.

I would like to use gdal to accomplish this task.


texts: Adding custom transformation parameters in Reproject raster chain function in Mosaic Dataset?

texts:
I have scanned maps georeferenced in custom projection with local datum. I have transformation parameters from the local datum to WGS84.

I created a mosaic dataset for the georeferenced maps with custom projection. Now I created another mosaic dataset referenced to first one to have new dataset with WGS84 datum.

I added Reproject function in raster chain functions and set the output projection as WGS84. There is no place to set transformation parameter from local datum to WGS84 and as a result the new raster dataset is not reprojected properly.

How can I overcome this issue?


texts: Disabling Extent Based Rendering using ESRI Javascript API

texts:
I am using ESRI's javascript API to generate a heatmap renderer for my feature layer. I am noticing that the heatmap calculation/visualization is dependent upon what is in the current zoom extent. For instance, when I observe the map at the World zoom extent, the map indicates there are clusters (see first attached image). However, when zooming into one of said clusters, the cluster disappears and is instead replaced with very faint highlighting (see second attached image).

I believe this is due to the fact that the renderer calculates the visualization of the heatmap based on which datapoints are contained within the current zoom extent. I am looking for a way to disable this functionality such that the visualization remains consistent regardless of what zoom level the user is on.


texts: Using GRASS v.delaunay with QGIS results with gaps in triangulation

texts:
When I use the GRASS v.delaunay algortithm in  QGIS 3.18 on my dataset the output has gaps in it. But when I use the QGIS delaunay algorithm there are no problems and all is triangulated.

The input that I use is a pointlayer with xyz geometry and the algorithms are called on from the processing toolbox. The settings for v.delaunay are shown below.

And these are the warnings.

Does anybody know how to solve this?


texts:
Is there a way to create a raster heatmap (aka 2d kernel density) from point data in R which can then be exported as a geoTIFF file?

I'm able to create a heatmap in R easily using ggplot2's geom_density_2d, but as far as I am aware this does not in any way allow me to export the resulting heatmap as a geoTIFF raster.


texts: Error when publishing Geoprocessing tool to ArcGIS Server

texts:
I've designed a model which compares some shapefiles to other shapefiles (I want to see if the files I get are the same as the files of my template).
When I try to publish the toolbox (where the model is found) on the server I get this message "Configuration can not be started. Server Object instance creation failed on machine....Geoprocessing service initialization failed. Invalid parameter: Tool..."

SOC has full access to the folder where my toolbox is so this can't be it. Is it because I'm working with shapefiles? Or why would the parameters be invalid?

I am using ArcGIS Server 10.0; ArcDesktop Editor 10.0.

Thank you


texts: Atmospheric correction for planet images using ArcGIS Desktop

texts:
I am very much new to the field of image processing. I have several images from planetscope. 20190722_155754_1049_3B_AnalyticMS is one of them. I don't know how to process these images in ArcGIS Desktop.

Can someone help me?


texts: QGIS 3: How to rename attribute/field names for Geopackage Layer

texts:
Issue: I need to rename several attributes/fields in a geopackage layer. How to do this?

Linked Issue: there are numerous pages on how to rename the attributes/fields of a shapefile. However there do not seem to be any pages specifically addressing the modification of attributes to geopackage layers, except a few about using Python scripts. Perhaps what is assumed at the moment (Feb 2019) is that you must be able to do Python programming in order to modify fields in Geopackage layers?

Before tagging this as a duplicate, please make sure you are not presuming that I am trying to modify a shapefile. I am trying to modify a geopackage layer.


texts: Error 000358 using SelectLayerByAttribut_management using escape characters

texts:
I am trying to use arcpy.SelectLayerByAttribut_management using the following code:

def filterData(self):
    self.filterBuilder()
    arcpy.MakeFeatureLayer_management (r"%s\tempclipprefilter"%self.workspace, "lyr")
    arcpy.SelectLayerByAttribute_management ("lyr", "NEW_SELECTION", '" ATTRIBUTE" = \'self.DataQuery\' ')
    arcpy.CopyFeatures_management("lyr", self.outFC)

I have also tried just

self.DataQuery

and

ATTRIBUTE - self.DataQuery

in place of

'"  ATTRIBUTE" = \'self.DataQuery\' '

I have looked at Runtime error : ERROR 000358: Invalid expression using SelectLayerByAttribute in ArcPy? but this has not solved my issue.

The code producing self.DataQuery is below:

    def filterBuilder(self):
        self.queryList = []
        if self.aField == 'true':
            self.queryList.append(" myDataColumn <> 'N' ")
        if self.bField == 'true':
            self.queryList.append(" class LIKE 'R%'")
        if self.cField == 'true':
            self.queryList.append(" class NOT LIKE 'L%'")
        if self.dField == 'true':
            self.queryList.append(" class LIKE 'CI%'")
        if len(self.queryList)>0:
            self.DataQuery = self.queryList[0]
            i=1
            while i < len(self.queryList):
                self.DataQuery= self.DataQuery + " AND %s" %self.queryList[i]
                i+=1
        else:
            self.DataQuery = ""

What is the error in my code or my approach?


texts: PHOTOMETRIC tag not producing expected result

texts:
Code takes a netCDF4 file and converts the data to GeoTiff:

import numpy as np, sys
from osgeo import gdal
from osgeo import osr
import netCDF4

f=netCDF4.Dataset("mydata.nc")

time = f.variables["time"]
lon = f.variables["lon"]
lat = f.variables["lat"]
data = f.variables["data"]
dims = f.dimensions
xsize = len(dims['lon'])
ysize = len(dims['lat'])

lonarray = np.array(lon)
latarray = np.array(lat)
dataarray = np.array(data)

# get bounds
xmin,ymin,xmax,ymax = [lonarray.min(),latarray.min(),lonarray.max(),latarray.max()]

# no. rows and columns
nrows,ncols = np.shape(dataarray)

# resolution of x/y
xres = (xmax-xmin)/float(ncols)
yres = (ymax-ymin)/float(nrows)

geot=(xmin,xres,0,ymax,0, -yres)

output_raster = gdal.GetDriverByName('GTiff').Create('myraster.tif',ncols,
nrows, 1 ,gdal.GDT_CInt16, options = ['PHOTOMETRIC=MINISWHITE'])  # Open the file

output_raster.SetGeoTransform(geot)            # Specify its coordinates

srs = osr.SpatialReference()                           # Establish its coordinate encoding

srs.ImportFromEPSG(4326)                               # This one specifies WGS84 lat long.

output_raster.SetProjection(srs.ExportToWkt())         # Exports the coordinate system
                                                       # to the file
output_raster.GetRasterBand(1).WriteArray(dataarray)   # Writes my array to the raster

output_raster.FlushCache()

del output_raster

After I load the Geotif into Arc Globe, the min value is black, ie. the  ['PHOTOMETRIC=MINISWHITE'] tag has not done what I expected.


texts: Sorting/unsorting fields displayed in Identify window of ArcMap?

texts:
I'm using ArcMap 10.1, and upon bringing in a feature class, the identify tool window shows the attributes in alphabetical order of the field headings.

I didn't try to do this and haven't touched any settings.

Is there a way to make them display once again in the order that they were created (as they still are in the att table)?


texts: Looking for specific polygon boundaries of DMA regions

texts:
I am trying to create a map in CartoDB and I need 2 specific polygons for it. The regions are:

DMA_id, DMA_Name
743 Anchorage, AK
744 Honolulu, HI

I searched every source I have and those two are missing in every file. Anyone may help here?

Edit: I managed to find a shapefile of all the areas in Honolulu. I join them and create a new, one polygon that contains all the Honolulu. So, only Anchorage remains.


texts: Averaging variables of smaller polygons to larger polygons in R using the SF package

texts:
Still getting used to the language of GIS with R so please bear with me.

I have this public dataset.

It has two levels of political boundaries: ADM1 and ADM2.

I am currently doing everything (i.e. plotting, calculations, etc.) for ADM2 which is each of the 2,440+ municipalities in Mexico because that is how the geometry is determined. However, I would like to calculate averages per states (ADM1) and start working with those by sort of "dissolving", if you will, the geometry of municipalities (ADM2) into geometries of (ADM1).

So far I have imported it:

library (sf)
library (dplyr)

mexi_sf<-read_sf( dsn = getwd()
       , layer = "MexicoCaseStudy"
       , stringsAsFactors = FALSE )

Fixed some labels and changed some classes:

 mexi_sf$ADM1NAME<-str_replace_all(mexi_sf$ADM1NAME,
                c("M.xico" = "México", "Nuevo Le.n" = "Nuevo León","San Luis Potos."="San Luis Potosí","Quer.taro de Arteaga"="Querétaro de Arteaga"))
mexi_sf$ADM1NAME<-as.factor(mexi_sf$ADM1NAME)
mexi_sf$ADM2NAME<-as.factor(mexi_sf$ADM2NAME)
mexi_sf$ADM1CODE<-as.numeric(mexi_sf$ADM1CODE)
mexi_sf$ADM2CODE<-as.numeric(mexi_sf$ADM2CODE)

Changed some outliers:

  for (row in 1:nrow(mexi_sf)){
    if (mexi_sf$FOODEXP[row]==-9999){
     mexi_sf$FOODEXP[row]=-1
    }}

Here is my attempt to average the variable FOODEXP by ADM1NAME.

  average_food_exp<-mexi_sf%>%
  group_by(ADM1NAME)%>%
  mutate(meanFOODEXP=mean(FOODEXP))

And obvs something is wrong because I can't plot it:

    big_food_map_sf_trial<-tm_shape(average_food_exp,style="quantiles",n=8)+
  tm_fill(col="meanFOODEXP",palette="Reds",title="mean FOODEXP (MXN 2004)")+
  tm_text("ADM1NAME",size=0.5)+
  tm_layout(legend.outside =FALSE)+
   tm_borders(col="black")+
  tm_credits("Source: CIMMYT, June 2004, Mexican municipalities below food poverty line, year 2000.",position=c("left","bottom"),size=0.75)

save_tmap(tm=big_food_map_sf_trial,filename="big_food_map_sf_trial.png",width = 17,height=7)

Any ideas on how to proceed?


texts:
I am bringing in leaflet routing control into my map and I would like to have the route icons show up in the center of the current view. using the code below I can get the icons to show up at the center of the inital view but I cannot find a way to get the center to update

Javascript:

    var route = L.Routing.control({
  waypoints: [
  L.latLng(map.getCenter()),
  L.latLng(map.getCenter())
  ]
});

JQuery:

$("#route-btn").click(function(event) {
  event.preventDefault();
  $('.leaflet-routing-container').is(':visible') ? route.removeFrom(map) : route.addTo(map);l
});

like most of the mistakes I make it is probably something simple, any help would be welcomed.


texts: How to set compression option correctly in GDAL raster calculator in PyQgs?

texts:
I perform some calculations on a set of approx. 700 rasters. The calculation works in my loop, but the new raster files are now 500 MB instead of 2 MB before the calculations. By default, there is no copression set in the GDAL raster calculator.
I have tried to add the compression information to the 'OPTIONS' input, but it does not work:
    calc_expression = '(255>A>={0})*A)'.format(some_value)
    parameters = {'INPUT_A' : raster_layer,
        'BAND_A' : 1,
        'FORMULA' : calc_expression,
        'RTYPE': 2,
        'NO_DATA': 0,
        'OPTIONS': '--co="COMPRESS=LZW"',
        'OUTPUT' : output_raster}

    processing.runAndLoadResults('gdal:rastercalculator', parameters)

I have also tried by calling gdal_calc.py directly as follows:
    pythonpath = 'C:/Program Files/QGIS 3.10/apps/Python37/Scripts/gdal_calc.py'
    gdal_calc = 'python {0} ' \
        '-A {1} ' \
        '--outfile={2} ' \
        '--calc={3} ' \
        '--overwrite'.format(pythonpath, raster_file, output_raster, calc_expression)

    os.system(gdal_calc)

But this won't create any raster at all (giving "1" as output).
Can anybody help me to set the compression option in the first code snippet correctly (priority 1), or alternatively to get the code in the 2nd snippet running (priority 2)?


texts: Copying set of points perpendicular along polyline using PyQGIS

texts:
I have two layers, 'river (line)' and  'points' for a set of points laying on a line which is perpendicular to the river (for example red points in the first image).

Given an interval (5.0 m for example), how can I copy those points along the river (blue sets of points in the next image)? But every set of points must be perpendicular to the related segment of river.

But I have just the following result using the next script.

from numpy import arange

L = lambda name: QgsProject.instance().mapLayersByName(name)[0]

point_layer = L("points")
river_layer = L("river")

points = list(point_layer.getFeatures())
river_geom = list(river_layer.getFeatures())[0].geometry()

# Find the nearest point on the river to the points
# The points are laying on the same line perpendicular to the river,
# it's sufficient to use any point.
nearest_point = river_geom.nearestPoint(points[0].geometry())

interval = 5.0
start = river_geom.lineLocatePoint(nearest_point) % interval
end = int(river_geom.length()/interval)*interval + start
for i in arange(start, end + interval, interval):
    p = river_geom.interpolate(i)
    dx = p.asPoint().x() - nearest_point.asPoint().x()
    dy = p.asPoint().y() - nearest_point.asPoint().y()

    feats = []
    for point in points:
        feat = QgsFeature()
        geom = point.geometry()
        geom.translate(dx, dy)
        feat.setGeometry(geom)
        feats.append(feat)

    point_layer.dataProvider().addFeatures(feats)

point_layer.updateExtents()
iface.mapCanvas().refresh()


texts: How to relate flood events (points) with flood conditioning rasters for Machine Learning?

texts:
I have point shapefile of flood events which has either flood occur or not.
On the other hand, I have 13 flood conditioning rasters i.e. slope, elevation, rainfall etc. these are all at 30 meter resolution.
Now I want to relate each point with these rasters spatially. I just want to confirm that should I take corresponding pixel values or there is a better approach.
Thanks.


texts: How to query for a specific point, within a specific polygon?

texts: I have a shapefile points layer, with each point representing either Marginal, suboptimal, or optimal conditions. I want to find only the marginal points within a specific polygon from my catchments polygon layer. So far I have been able to select all marginal points, or the wanted polygon, but not the marginal points WITHIN the specified polygon. Is this possible? Thank you.


texts: Calculating volume of building in ArcGIS Pro?

texts:
I want to calculate the volume of buildings in arcgis.
The data consists of 3D multipatches (shapefile).

One building consists of different multipolygons. I can calculate the volume of the single multipolygons with the built-in functions of the 3D-Analyst (minimum bounding volume, convex hull).

The problem is that one building is composed of several multipolygons that aren't related to the building via an ID.

Therefore I need to do a spatial match between the different polygons of a building.

Any suggestions of how to do that using ArcGIS Pro?

I was not able to execute 3Dunion of the 3D analyst).


texts: Unsupported Operand Type for Field Query

texts:
Following on from Correct SQL expression for SelectbyAttribute Tool, I have my code:

import arcpy

plantFile = r"U:\Users\K\Plants.shp"
field = "UID"

arcpy.MakeFeatureLayer_management (plantFile, "plant_temp")
expression = round("UID" / 2, 0) <> "UID" / 2
arcpy.SelectLayerByAttribute_management("plant_temp", "NEW_SELECTION", expression)

I encounter -- TypeError: unsupported operand type(s) for /: 'str' and 'int'

"UID" is a long integer field, so I'm not understanding why this is happening. I've looked up methods for changing types for individual values, but not for an entire field. I will also be running this on over 20,000 records, so processing time is a concern.

EDIT
I applied the suggested expression expression = '''round("UID" / 2, 0) <> "UID" / 2'''

which returned ExecuteError: ERROR 000358: Invalid expression


texts: Get information from geotiff with lat/lon

texts:
I have a geotiff file containing data like temperature, now I need to extract the temperature if a real-world lat/lon pair is given. I am using python and GDAL, following is a code I found and it worked.

However, when I used QGIS (I used value tool plugin) to validate the result, for a given lon/lat, the output value from QGIS is not consistent with the following code. The cell of my code is actually some 1- or 2-distance neighbors of what QGIS looks into. I don't see any problem of my code, is QGIS value tool has some precision issues that I should not trust?

from osgeo import gdal,ogr
from osgeo.gdalconst import *
import struct
import sys

lat = float(sys.argv[2])
lon = float(sys.argv[3])

def pt2fmt(pt):
    fmttypes = {
        GDT_Byte: 'B',
        GDT_Int16: 'h',
        GDT_UInt16: 'H',
        GDT_Int32: 'i',
        GDT_UInt32: 'I',
        GDT_Float32: 'f',
        GDT_Float64: 'f'
        }
    return fmttypes.get(pt, 'x')

ds = gdal.Open(sys.argv[1], GA_ReadOnly)
if ds is None:
    print 'Failed open file'
    sys.exit(1)

transf = ds.GetGeoTransform()
cols = ds.RasterXSize
rows = ds.RasterYSize
bands = ds.RasterCount #1
band = ds.GetRasterBand(1)
bandtype = gdal.GetDataTypeName(band.DataType) #Int16
driver = ds.GetDriver().LongName #'GeoTIFF'
transfInv = gdal.InvGeoTransform(transf)
px = (lon-transf[0])/transf[1]
py = (lat-transf[3])/transf[5]
#px, py = gdal.ApplyGeoTransform(transfInv, lon, lat)
structval = band.ReadRaster(int(px), int(py), 1, 1,  buf_type = band.DataType )
fmt = pt2fmt(band.DataType)
intval = struct.unpack(fmt , structval)
print "value: ", intval[0]


texts: Cross-schema function call with PostGIS dependencies

texts: I have an issue where I am unable to make a cross-schema function call because of that function's use of PostGIS. I have both PostGIS and my function in the public schema, but I'm trying to call the function from my ORM which, by necessity, is pointed at another schema. I am able to reference the function by using SELECT * FROM public.functionName, but the function execution fails when it reaches PostGIS-relevant functions like ST_MakeLine. What is the best solution in this case, move PostGIS and my function out of public? Edit the function itself to call public.PostGISFunction?  Something else entirely?


texts: SPOT and LANDSAT images on the same grid

texts:
Let us introduce ourselves, we are a GIS team of the Cartographic division of INDEC (the national census bureau of Argentina). We are working on the next national Agricultural Census 2018, which will be carried out on Tablets with Android O.S., running Leaflet to show the maps.

We have strong storage restrictions (SD 16 GB memory) and the whole program (images included) has to run local on the tablet (no internet signal in the most areas to be covered).

We get a province-wide 30mts, LANDSAT mosaic, resampled at 15mts with the monochromatic band (each pair: census taker/tablet) are restricted to only one province.

But for some zones (intensive agricultural use) we need SPOT 2,5 mts images.

We get to generate the tiled pyramid of the 15mts landsat, using GDAL gdal2tile

gdal2tiles.py -t "BS AS" -c INDEC -z 0-8 -e -d -a 0 -p raster -v -s EPSG:22185 -r cubic /BsAsGK5.tif /tiles_bs_as

After that we convert the png files to jpg, using:

find tiles_bsas -iname '.png' | while read i; do mogrify -format jpg "$i" && rm "$i"; echo "Converted $i to ${i%.}.jpg"; done

Here is our question: How should we generate the tiled pyramids of the 2,5 mts SPOT images in a way that both images use the same grid without 'resampling' the whole landsat to 2,5mts, of course. :)

If you are willing to give us some advice, we will send you the links to some sample images.

Link to LANDSAT 15mts. example (prov of BsAs) 6.5GB

Link to SPOT 2,5mts. example (Department of San Pedro) 2.3GB


texts:
I created slope map with DEM model in ArcGIS Desktop. But it is only categorize in to two sections. I think it is problem with the horizontal and vertical scale. I used WGS84 coordinate system.

How do I overcome this problem?


texts: Get lat/lon point of intersection a route and a polygon object on google maps

texts:
I have a shape (polygon object) and a route, such as on the attached google map.

How I get waypoint (lat/lon, C point on map) where intersect the route and the polygon object?


texts: Adding raster layers to an existing GeoPackage in R

texts:
I have a GPKG file with a few vector layers. I have a few processing steps in my script and every now and then a new layer is added to the file. This works fine for vector layers using something like

sf::st_write(vectorlayer, "output.gpkg", "layername", update = TRUE)

but I have no idea how to add raster layers. I can use raster::writeRaster to create a new GPKG, but I'm not able to extend an existing one. Any ideas what I'm missing here?


texts: Symbol missing in legend when adding to layout map

texts:
I am new to QGIS and I am trying to make a map of India displaying some values. I want to have a legend, and I customised the legend in layers panel, the Layers panel clearly displays my legend. But the issue is when I try to make a map  with that, legend displays no symbols, only values are displayed with the legend. I am adding the screenshot of my Layers panel showing the symbols and layout map without the symbols. I am using QGIS version 3.18.1.
Layers panel showing symbols:

Layout map legend without symbols:


texts: How to re-project only one band from a hdf?

texts:
I have this file which contains:

gdalinfo product.hdf
Coordinate System is `'
Metadata:
calendar=standard
center=gsfc
Subdatasets:
Corner Coordinates:
Upper Left  (    0.0,    0.0)
Lower Left  (    0.0,  512.0)
Upper Right (  512.0,    0.0)
Lower Right (  512.0,  512.0)
Center      (  256.0,  256.0)

as you can see, there are 24 bands What I need is just to extract the first band  of the sub-dataset.


texts: What does 'building topology' mean for pgRouting?

texts:
From a couple of answers at gis.stackexchange.com and at the pgRouting mailing list I have found out that pgRouting 'builds topology'. What does that mean?

If for example we import an .osm file into PostGIS and then apply all the pgRouting scripts, what does happen?

I thought that pgRouting just creates a table of road links (edges) from the OSM nodes, but after following the FOSS4G workshop I see the following tables in my database:

List of relations
 Schema |       Name       | Type  |  Owner
--------+------------------+-------+----------
 public | classes          | table | postgres
 public | geometry_columns | table | postgres
 public | nodes            | table | postgres
 public | spatial_ref_sys  | table | postgres
 public | types            | table | postgres
 public | vertices_tmp     | table | postgres
 public | ways             | table | postgres

I know what geometry_columns and spatial_ref_sys serve for, but what about the other tables and how do they relate to each other?

What is topology in pgRouting's parlance?


texts: ArcGIS Pro 2.5 Python Package Manager not displaying

texts:
I am trying to install some new packages to work with and my Python Package Manager in ArcGIS Pro 2.5 is completely blank; will not display. It was displaying last week when I was working from home, but now it will not show on my work network. I even tried to pull it up at home last night but it would not show. I haven't changed anything that would effect the installation files to my knowledge. Has anyone else run into this, or at least the instance where the manager would not display?

Here is all I see:


texts: Shapefile will not load in QGIS but does in ArcGIS

texts: Using QGIS version 3.14.0-Pi on Windows 10.  I'm trying to download a shapefile (Layer, Create Layer, New Shapefile layer).  In the resulting dialog I use: File encoding UTF-8, Geometry type "Polygon" Add. dimensions EPSG:4269 - NAD83.  The original file is in NAD 83 UTM 12N.  The map I"m using is using EPSG:4326 WGS 84. When I say OK I get a message "The layer already exists.  Are you sure you want to overwrite the existing file?"  This happened the first time I tried to make the layer and occurs each time I make a new attempt.  When I say "yes" I get a new layer but nothing on the map and the attribute file is empty.  The same shapefile opens without any issues in ArcGIS.


texts: Open directories in interaction with a Shapefile

texts:
I am trying to write a python code whose purpose is to open windows directories corresponding to the communes I select in a "common" Shape layer (with Qgis software).
I specify that the common layer contains a column X with the paths to the directories of the communes (example of the path: V: / 02_doc / commune1 /)

Here is my current code. It does not work as intended (it opens ... the windows directory "'my documents"

import os
X = os.path.dirname( unicode( qgis.utils.iface.activeLayer().dataProvider().dataSourceUri() ) )
os.system("explorer " + "[%"X"%]")
layer = qgis.utils.iface.activeLayer()
qgis.utils.iface.setActiveLayer(layer)

this code is used to indicate a single path and point on when I click on the layer.

import subprocess
layer = iface.activeLayer()
x = layer.source()
subprocess.Popen(r'explorer /select,' + x)

that's not what I want.
fcka

`with this code:

import os
X = os.path.dirname( unicode( qgis.utils.iface.activeLayer().dataProvider().dataSourceUri() ) )
os.system("explorer "+posplu)
layer = qgis.utils.iface.activeLayer()
qgis.utils.iface.setActiveLayer(layer)

there is a change, when I click on my layer in Qgis, I point the directory where my layer is.

that's not what I want.

My layer contains an x column.
The latter includes paths pointing repertory directories to cities. I would like to click on a city A and open the directory A, on a city B and open the directory B.

fcka

0
down vote
accept

i found

import os layer = qgis.utils.iface.activeLayer() import subprocess subprocess.Popen(r'explorer /select r,"[% X %]"')

x is the column containing directory paths fcka


texts: Calculating ratio of polygon area inside another polygon area (ArcInfo 10.5)

texts:
I am trying to calculate cover ratio for each plot in a neighbourhood and thus there are 2 shapefiles, one with parcel and other with building outline.

So I want to:

Assign each building a plot according to location. I would need some tolerance option since some buildings slightly overextend plots.
Add up the area of building in a plot.
Divide by plot area.

I am unable to think how to do step one.

I am using ArcInfo, with Spatial Analyst License.


texts: Select points North of line across polygon using ArcGIS for Desktop?

texts:
I need to select all points north of a line. I'm have trouble conceptualizing what tools would allow me to do this. Any advice? I'm in ArcGIS 10.2.2 if that matters.

In the poorly drawn picture I would like to select all red points north of the green line. Any advice?


texts: Write stacks of GeoTIFF image to NetCDF by chunk size

texts:
I converted GeoTIFF image stacks to NetCDF file with writing data by Index order (layer by layer) (this was a main solution).
The problem is, it writes slowly (for 25 images - 1Gb takes 90 seconds). One of the solution is to write the data chunk by chunk not the whole image at once.
Does somebody know how to do it or any other solution !?

Code (I used memory to keep images as much as possible then write them all at once)

def tiff_to_netcdf(imagesPath, output, chunksize=None):

    images = []
    for paths, subdirs, files in os.walk(imagesPath):
        for name in files:
            images.append(os.path.join(paths, name))

    ds = gdal.Open(images[0])

    a = ds.ReadAsArray()
    nlat, nlon = np.shape(a)

    b = ds.GetGeoTransform()  # bbox, interval
    lon = np.arange(nlon)*b[1]+b[0]
    lat = np.arange(nlat)*b[5]+b[3]

    basedate = dt.datetime(2006, 01, 11, 0, 0, 0)

    # Create NetCDF file
    # Clobber- Overwrite any existing file with the same name
    nco = netCDF4.Dataset(output+'.nc', 'w', clobber=True)

    # Create dimensions, variables and attributes:
    nco.createDimension('lon', nlon)
    nco.createDimension('lat', nlat)
    nco.createDimension('time', None)

    timeo = nco.createVariable('time', 'f4', ('time',))
    timeo.units = 'days since 2006-01-11 00:00:00'
    timeo.standard_name = 'time'

    lono = nco.createVariable('lon', 'f4', ('lon',))
    lono.standard_name = 'longitude'

    lato = nco.createVariable('lat', 'f4', ('lat',))
    lato.standard_name = 'latitude'

    # Create container variable for CRS: lon/lat WGS84 datum
    crso = nco.createVariable('crs', 'i4')
    crso.long_name = 'Lon/Lat Coords in WGS84'
    crso.grid_mapping_name = 'latitude_longitude'
    crso.longitude_of_prime_meridian = 0.0
    crso.semi_major_axis = 6378137.0
    crso.inverse_flattening = 298.257223563

    # Create short integer variable with chunking
    tmno = nco.createVariable('tmn', 'i2',  ('time', 'lat', 'lon'),
                              zlib=True, chunksizes=chunksize, fill_value=-9999)
    tmno.scale_factor = 0.01
    tmno.add_offset = 0.00
    tmno.grid_mapping = 'crs'
    tmno.set_auto_maskandscale(False)

    nco.Conventions = 'CF-1.6'

    # Write lon,lat
    lono[:] = lon
    lato[:] = lat

    pat = re.compile(r'^(\w{1})(\d{4})(\d{2})(\d{2})_\d{6}___SIG0.*.tif$')
    itime = 0
    written = 0
    data = np.empty((3, 8000, 8000), dtype=np.float)

    # Step through data, writing time and data to NetCDF
    for root, dirs, files in os.walk(imagesPath):
        dirs.sort()
        files.sort()
        for f in files:

            match = re.match(pat, f)
            if match:

                if itime  % 3 == 0 and itime != 0:
                    tmno[written:itime, :, :] = data
                    written = itime
                # read the time values by parsing the filename
                year = int(match.group(2))
                mon = int(match.group(3))
                day = int(match.group(4))
                date = dt.datetime(year, mon, day, 0, 0, 0)
                print(date)
                dtime = (date-basedate).total_seconds()/86400.
                timeo[itime] = dtime
                tmn_path = os.path.join(root, f)
                print(tmn_path)
                tmn = gdal.Open(tmn_path)
                a = tmn.ReadAsArray()  # data
                data[itime%3, :, :] = a
                itime = itime+1

    nco.close()


texts: Tile rendering time measurement in MapServer?

texts:
I try to speed up rendering time of tiles in mapserver app. How can I specify rendering engine in mapserver? The only thing about I found is: outputformat

Tiles are generating very slow - about 4-5 sec (with postgres query - about 1,5 sec).
How can I measure rendering time for tiles and what to do to speed up this process?


texts: Python Script Zonal Stats as Table Loop Question

texts:
I'm trying to write a small script loop script which retrieves rasters from each folder, performs zonal statistics using watershed polygons and saves it as a table. I've written similar loop scripts, but this seems to be confusing me. Here's what I have so far:

import arcpy
from arcpy import env
from arcpy.sa import *
env.workspace = "F:/lu_crop_rast"
arcpy.env.overwriteoutput = 1

watershedFeat = "F:/watersheds.shp"
outDir = "F:/lu_crop_rast/tables/"

for raster in arcpy.ListRasters():
    outTable = outDir + raster + "_TBL.dbf"
    arcpy.sa.ZonalStatisticsAsTable(watershedFeat,"Name",raster,outTable,"NODATA","MEAN")

This is the error I'm getting:

Runtime error  Traceback (most recent call last):
File "<string>", line 15, in <module>   File "c:\program files (x86)\arcgis\desktop10.3\arcpy\arcpy\sa\Functions.py", line 6232, in ZonalStatisticsAsTable     statistics_type)
File "c:\program files (x86)\arcgis\desktop10.3\arcpy\arcpy\sa\Utils.py", line 53, in swapper     result = wrapper(*args, **kwargs)   File "c:\program files (x86)\arcgis\desktop10.3\arcpy\arcpy\sa\Functions.py", line 6224, in Wrapper     statistics_type)
File "c:\program files (x86)\arcgis\desktop10.3\arcpy\arcpy\geoprocessing\_base.py", line 504, in <lambda>     return lambda *args: val(*gp_fixargs(args, True)) ExecuteError: ERROR 999999: Error executing function. Create output table failed Failed to execute (ZonalStatisticsAsTable).

I've checked and my rasters and shapefile aren't corrupt. I've looked at this thread as well and made the following changes:

arcpy.gp.ZonalStatisticsAsTable_sa(watershedFeat,"Name",raster,outTable,"NODATA","MEAN")

But so far nothing seems to work.


texts: Trying to use the Least Cost Path plug in in QGIS

texts: I am trying to calculate the least cost path among one central point (start point layer) and another layer made of several points, representing other sites (end-point(s) layers). I am using the Least Cost Path plug in on QGIS, and have produced a slope layer for my study area, representing my Cost raster layer. When I run it, this comes up:  'The end-point layer contains no legal point'. I have no idea what that means, can someone explain it to me? I have researched online a little bit, and it says that on QGIS that it is not possible to find the least cost path when there's multiple points. If that is true, why does this plug in say point(s)?


texts: Assign factor values and make calculations by using Corine raster dataset as basis in ArcMap

texts:
I am using ArcMap.
I have the Corine raster map which is consolidated by the the class codes (example 121) and some additional information.
The purpose is to calculate the Crop/Management factor for importing it in the RUSLE model.
What I would like to do to achieve this calculation is to:

Assign manually the values of the factor (according to literature) when the code belongs to the arable classes (which can be done easily by adding a field in the attribute table and then edit it according to the code) &
(The difficult part) Calculate the factor for the non-arable classes (code from 221-334) by using NDVI.

In particular the formula is C = minC + range(C)* (1-NDVI)
However, each class/code has a different min, range and NDVI value so I cannot just use the raster calculator.
I need to select somehow the value of each code and according to the respecting min and range values to calculate the formula by using the NDVI that corresponds to each pixel.


texts: Color correcting/balancing for mosaic in ArcGIS for Desktop?

texts:
I have about a hundred of individually georefered RGB aerial images which I'd like to mosaic to have a nice, seamless final raster. I need a method to make the mosaic lines as invisible as possible. l Here is a pic:

I tried simply mosaic it, but lines would be clearly visible like above. The best result so far is with mosaic color correction, with dodging, color grid option:

But it's still not the best.

I tried ERDAS Imagine MosaicPro with dodging and or color correcting, but didn't get an acceptable result.

Is there a way to do a good looking mosaic in ArcGIS for Desktop?

I have an Advanced level license available.


texts:
I have a big raster image with medium details, but for some regions I do have far more details as a seperate raster image.

Now I'd like to display both of them in my map that I render with Tilemill.
So now I would need to cut a hole out of the bigger raster (exactly the shape of the small one) file so that both layers will be properly displayed on my map. (Note: I don't display it with an opacity of 100%)

I found that using the gdalwarp -cutline shape.shp input.tif output.tif command I can cut a raster. But I need the inverse of that.
Is there a parameter for gdal which I am missing?

Thanks a lot


texts: Supervised Classification Plugin QGIS: change input image; use another raster layer

texts:
I'm doing supervised classification on a sugar cane field for identifying vegetation, soil, and soil and litter.

The plugin is only taking in the GeoTiff taken from the multispectral camera we used (RedEdge, MicaSense) for creation of ROIs and the training shapefile. However, this is not the best input I can use because parts of the field have plenty of shadows from adjacent trees. Thus, the shade throws off my ROIs and my classification.

I have a drastically better view of vegetative biomass if I use the raster calculator to display NDVI ((Band 5 - Band 3)/ (Band 5 + 3)). And would greatly improve my classification process if I could use that raster layer instead of the original GeoTiff raster layer. Alas, the SC-Plugin will only let me use the GeoTiff for the classification process.

I imagine the plugin only takes in <.tif> files. My NDVI is <.tif.aux.xml> so I tried saving it as a Tif file. However, QGIS returned an error message saying that my raster layer (a .tif file) was not a supported data source.

How can I get the Supervised Classfication Plugin to accept my NDVI raster layer?


texts: Deleting specific vertex in polygon in ArcGIS Desktop?

texts: I use ArcGIS 10.3 and I have many  Building consist of window and wall but my problems that windows are not on straighter the wall so I search tool help me to delete corner vertex of window because i have more buildings


texts: Why EPSG codes in the QGIS database are not the same as in the EPSG registry?

texts:
QGIS (current 3.4.4 version) have almost all its CRS definitions as EPSG Authority codes.

But many of them are not registered in the https://www.epsg-registry.org/ database, so many of them really have not an EPSG Authority.

By the side, at this moment that lack of support is causing QGIS not recognizing the CRS of the layers exported to non-EPSG authority codes, as can be seen in the following posts:
 - Several projection issues QGIS 3.4.3
 - CRS information not recognized when geopackage is loaded

I'm trying to understand the QGIS CRS codes (sr_id and auth_name fields of the tbl_srs table, in srs.db) database, and how QGIS writes and reads the CRS of the layers and try to match them with its database.

And what I can not find is the reason: Why almost all the codes in the QGIS database carry the EPSG Authority when really many of them are not in the EPSG registry?


texts:
I have a list of 18 properties, which are all flats. They are split between 3 blocks, 6 flats each.
Therefore there are 6 points with the same co ordinates, but each of the flats has their own attributes.

In number “Block 2”, two of the six flats have different attributes to the other four.
(e.g. different number of bedrooms per flat)
When I use the symbology tool within the Layer Properties to show the differences within,  it is only going to show me whichever is on top.

I would like to find a way to represent this?

I know I could move these manually, but with the amount of properties I am going to have to do this to it is not the way forward.

Using: Arcmap 10.1


texts: Caching layers on the fly and seeding combination

texts:
I have created a map in geoserver Im am using geoserver 2.5.1 with the intergrated geowebcache.

The map seems to be working fine and the seeding also seems to be working fine provided I don't use a bounding box for some strange reason.
But also want it to cache on the fly as well so I can seed the main sections of the map that would be used and have on the fly caching in the areas that arn't being used as much.

Am I mistaken in thinking that it saves those certain tiles like the seeding does or does it do it virtually or save in a temp tile folder and then goes away instantly tht person either closes the map program.

I have set everything that I could set but for some reason on the fly caching does not seem to work.

Am i missing something in the settings?
Or am misunder standing something with regards to it as in the geoserver documention states that it can combine both on the fly caching as well as seeding the layers.


texts: ArcGIS API - Creating Popups for Geoprocessing Result

texts:
With the ArcGIS API for Javascript I created a web application that accesses my before published geoprocessing service. The result of the service are a set of points and are displayed on the map. Here is the way how I display the points (symbol is a simple Marker Symbol):

var features = results[0].value.features;
for (var i=0; i<features.length; i++)
{
  var feature = features[i];
  feature.setSymbol(symbol);
  map.graphics.add(feature);
}

Now I would like to enable popups for these points, which include some further attribute information that can be accessed with e.g. feature.attribute.name.

How can I, for these points, add a popup?


texts: Is it possible to customize the identify form in QGIS (I am using version 3) to create Balloon like windows and add a picture in addition to the features extracted from the table?


texts: How to prevent RGB composite flatten to 1 band after raster calculation

texts:
I have a series of landsat wanting to process in raster calculator, but after doing the tool it happened these output images only contain one single band. They used to be RGB 3-band images.

It is on arcmap 10.3. I was using raster calculator to remove a certain value in my 3-band landsat image. After I type in the conditional statement to remove all 0 value from the inpur RGB image, it only returns an image of 1band, black and white. How can I preserve the RGB value in output.


texts: Is there a raster data format that can store dates instead of numeric values?

texts:
When storing spatiotemporal data, there is a need to store dates or time data. I know raster files can store numeric data (integer or float), but what about storing of dates data?
I'm currently looking at the following options:

Raster — Data are numeric values of the passed seconds from the Unix Epoch, 1970-01-01 00:00:00 being the zeroth second
ENVI .hdr labelled raster — Labels are in YYYY-MM-DD date format
Raster + CSV file — CSV file contains YYYY-MM-DD date format along with numerical value lookup for the raster
Raster — Data are numeric day-of-the-year, negative values and values larger than 365 are offsets and will be treated a different year

Is there a raster data format that can store dates instead of numeric values or utilizing labels and/or adding extra files?
If there aren't, at least, is there a de facto convention of how to store them on a raster file?


texts: Displaying point data to reflect time/age

texts:
I have a layer of 3,000 points for a state, with each indicating the age of a specific type of building. What is a recommended visual display that effectively and equally highlights the age of these structures.

I thought of using a heat map, but that could lead to some bias, as in a small vicinity you can have 3 new buildings and 1 old one or vice versa taking away from the truthfulness. However, I'd like a visual aid that when looked at briefly or a snapshot can tell a story of the age of these structures.


texts: How to refine a gpx track and set additional points on existing road network? (aka map-matching)

texts:
I am trying to refine/populate gpx tracks recorded by a GPS device. I would like to add 1 or 2 points between each gpx point recorded and snap these additional points to roads.
I am using essentially QGIS. Any suggestion how I could do this?
This issue is also discussed in this topic how-to-create-an-exact-tracking-route-from-gps-data.
As mentionned in one answer, the problem is called map matching.
For reference: some ideas are discussed on this post map-matching-links-and-ideas


texts:
I have a field "istat" that contains multiple values of the istat code (municipal code):

and a table containing the centroids of the municipalities

I am looking for a way to join the "istat"-fields of the two tables in order to create a table with multipoint geometry with the centroids of the common.
I thought of creating as many "istat" columns as the maximum number of "istat"-codes present in a cell, doing a multiple join.
Finally combine all the geom fields creating a multipoint field (if it is possible to do such a thing).


texts:
I'm trying to load multiple shapefiles as GeoJSON and display them on OpenLayers.

This is the code:

        var map = new ol.Map({
        layers: [new ol.layer.Tile({
            source: new ol.source.OSM()
        })],
        target: 'map',
        controls: ol.control.defaults({
            attributionOptions: {
                collapsible: false
            }
        }),
        overlays: [overlay],
        view: new ol.View({
            projection: 'EPSG:4326',
            center: center,
            zoom: zoom,
            rotation: rotation
        })
    });

    loadshp({
            url: './10tnvillage.zip',
            encoding: 'big5',
            EPSG: 3826
        },
        function (data) {

            var vectorLayer = new ol.layer.Vector({
                source: new ol.source.Vector({
                    features: (new ol.format.GeoJSON()).readFeatures(data,{
                        featureProjection: 'EPSG:4326'
                    })
                }),
                style: setStyle
            });

            map.addLayer(vectorLayer);
        });

    loadshp({
            url: './abadan.zip',
            encoding: 'big5',
            EPSG: 3826
        },
        function (data1) {

            var vectorLayer1 = new ol.layer.Vector({
                source: new ol.source.Vector({
                    features: (new ol.format.GeoJSON()).readFeatures(data1,{
                        featureProjection: 'EPSG:4326'
                    })
                }),
                style: setStyle
            });

            map.addLayer(vectorLayer1);
        });

It's show the Tile layer and the latest shapefile, but not all of them.

I did used a Array List and pass it to Map object, but didn't work.

        var layers = [];
    layers.push(new ol.layer.Tile({
        source: new ol.source.OSM()
    }));
    ...
    layers.push(vectorLayer);
    ...
    layers.push(vectorLayer1);
    ....
    var map = new ol.Map({
        layers: layers,
    ...

So, how can I display all my layers in OpenLayers?


texts: What Are the Minimum KML File Requirement to Draw Multiple Polygons?

texts:
My application is trying to export a KML file. The file should contain a series of polygons that represent sales territories. The territories cover the whole of the US.

I've followed the instruction for creating polygons on the KML tutorial. The file output from my application can be found here. I've tried to import this into QGIS and MyMaps. Both application parse the file but only show some of the territories.

My question is what am I missing? What are the minimum fields required to show multiple shapes with multiple polygons (and holes).


texts: Set Extents of MXD to new layer data source?

texts:
Alright here I go again this forum has been good to me. I add a layer and replace its data with a shapefile I created using the following bit of python:

layer = arcpy.mapping.Layer("C:/GIS/images.lyr")
arcpy.mapping.AddLayer(df, layer)
lyr = arcpy.mapping.ListLayers(mxd,"images",df)[0]
lyr.replaceDataSource(projectlocation+projectname+"/Aerial/"+"/"+projectdate,"SHAPEFILE_WORKSPACE","images")
lyr.name = "Images"

I cannot figure out how to set the extents of the mxd to the new source data I can just get it to set it to the original extents of the shp the .lyr from saved from. Any help is greatly appreciated.


texts: Only half of intersection layer saved due to OGR error

texts:
I am trying to intersect two shape files, one is a vector of .5 x.5 degree rectangles, the other is a vector of 20 x 10 degree rectangles.

I want the ID values of the 20 x 10 degree polygons to be carried over to the .5 x .5 degree polygons which are contained within the larger shapes.

The intersection produces what looks like the product I want for the whole globe, however when saving this I get the following error:

Feature creation error (OGR error: Attempt to write non-polygon (LINESTRING) geometry to POLYGON type shapefile.)

AND half of the output is successfully saved, but not the remainder (see image).


texts: $area gives NULL for some polygons in QGIS's Field Calculator

texts:
I have the map of Cambodia(map) and the economic land concessions(elc). In trying to find the proportion of district covered by land concesions, I need the area variable. But when I create the variable as $area, there are some entries that are NULL - it is the case for both layers. I have tried 'fix geometries.' It mostly fixed the problem for elc layer but still is NULL for one layer as following. Smaller piece and the bigger piece are seperate entries but area is calculated only for the small area.

The area calculation issue is not fixed at all for the elc layer. What is the next thing I could try to fix this problem?
I am following the steps from this link:
Finding overlay percentage between two polygon shapefiles
Data using:
GEO2_KH (Cambodia shapefile)
https://international.ipums.org/international/gis_harmonized_2nd.shtml
ELCs: All economic land concessions shapefile
https://data.opendevelopmentcambodia.net/dataset/economiclandconcessions


texts: How to remove the leaflet L.control.layers

texts:
I'd like to remove the layers control from my map, but I just can't find the solution...I tried almost everything I found on the internet

This what i tried so far:

var layerscontrol = L.control.layers(baseLayers, overlays).addTo(map);
...
layerscontrol.removeFrom(map);

or

layerscontrol.remove(map);

or

map.removeControl(layerscontrol);

what have I done wrong?


texts:
It is possible to compute shortest paths with GeoTools with its graph package.
Is there a package (or extension) for computing isochrone lines or areas, that is the area which can be reached from a point in less than X minutes by travelling on a transport network ?


texts: Raster mosaic error: Format driver GMT does not support creation and piecewise writing

texts:
I am novice in using QGIS and other GIS software and am trying to create mosaics of NOAA raster nautical charts.  Using pulldown menus, I entered raster>misc>merge and chose two NOAA charts and then specified a output file. The resulting command is shown below.

When I submit the command I then get the following error. Doing a web search, I found someone describing a similar problem with the solution being to specify the output file as GTiff.  When I specified the name testmosaic.GTiff, the command actually appeared as testmosaic.GTiff.tif.

Two questions:

If saving the file as GTiff is the solution to my original error,
how do I specify that format? I don't see how I can specify that in
the menu screen.
Is there some better way to get the results I want and avoid the
error

Command:

gdal_merge.bat -co COMPRESS=PACKBITS -of GMT -o "C:/Users/hamp/Documents/BC & AK/NOAA RNC ALaska/large area charts/testmosaic.tif" "C:\\Users\\hamp\\Documents\\BC & AK\\NOAA RNC ALaska\\large area charts\\Grouped Files\\17400_1.KAP" "C:\\Users\\hamp\\Documents\\BC & AK\\NOAA RNC ALaska\\large area charts\\Grouped Files\\17420_1.KAP"

Error:

 Format driver GMT does not support creation and piecewise writing.
 Please select a format that does, such as GTiff (the default) or HFA (Erdas Imagine).


texts: Convert ETRS89_Portugal_TM06 coordinate to Google Map Coordinate

texts:
I am trying convert an coordinate (point = {-26841,7714704367;42953,3965652831}) from shapefile with follow prj file:

PROJCS["ETRS89_Portugal_TM06",GEOGCS["GCS_ETRS_1989",DATUM["D_ETRS_1989",SPHEROID["GRS_1980",6378137,298.257222101]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",39.66825833333333],PARAMETER["central_meridian",-8.133108333333334],PARAMETER["scale_factor",1],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["Meter",1]]

To a Google Maps coordenate. I am reading about it and I found that I need to transform this coordenate into EPSG 3857. I am using the follow code to do that:

const string latlongwkt = "PROJCS[\"Popular Visualisation CRS / Mercator\", GEOGCS[\"Popular Visualisation CRS\", DATUM[\"WGS84\", SPHEROID[\"WGS84\", 6378137.0, 298.257223563, AUTHORITY[\"EPSG\",\"7059\"]], AUTHORITY[\"EPSG\",\"6055\"]], PRIMEM[\"Greenwich\", 0, AUTHORITY[\"EPSG\", \"8901\"]], UNIT[\"degree\", 0.0174532925199433, AUTHORITY[\"EPSG\", \"9102\"]], AXIS[\"E\", EAST], AXIS[\"N\", NORTH], AUTHORITY[\"EPSG\",\"4055\"]], PROJECTION[\"Mercator\"], PARAMETER[\"semi_minor\",6378137], PARAMETER[\"False_Easting\", 0], PARAMETER[\"False_Northing\", 0], PARAMETER[\"Central_Meridian\", 0], PARAMETER[\"Latitude_of_origin\", 0], UNIT[\"metre\", 1, AUTHORITY[\"EPSG\", \"9001\"]], AXIS[\"East\", EAST], AXIS[\"North\", NORTH], AUTHORITY[\"EPSG\",\"3785\"]]";
IProjectedCoordinateSystem googleProjectedCoordinates = CoordinateSystemWktReader.Parse(latlongwkt) as IProjectedCoordinateSystem;

var wktstring = "PROJCS[\"ETRS89_Portugal_TM06\",GEOGCS[\"GCS_ETRS_1989\",DATUM[\"D_ETRS_1989\",SPHEROID[\"GRS_1980\",6378137,298.257222101]],PRIMEM[\"Greenwich\",0],UNIT[\"Degree\",0.017453292519943295]],PROJECTION[\"Transverse_Mercator\"],PARAMETER[\"latitude_of_origin\",39.66825833333333],PARAMETER[\"central_meridian\",-8.133108333333334],PARAMETER[\"scale_factor\",1],PARAMETER[\"false_easting\",0],PARAMETER[\"false_northing\",0],UNIT[\"Meter\",1]]";
IProjectedCoordinateSystem sourceCS = CoordinateSystemWktReader.Parse(wktstring) as IProjectedCoordinateSystem;

CoordinateTransformationFactory ctFactory = new CoordinateTransformationFactory();
ICoordinateTransformation transformer1 = ctFactory.CreateFromCoordinateSystems(sourceCS, googleProjectedCoordinates);

double[] toPoint = transformer1.MathTransform.Transform(fromPoint);

But it still give me a wrong coordinate.

I also tried to replace IProjectedCoordinateSystem googleProjectedCoordinates to:ICoordinateSystem geographicCoordinates = ProjNet.CoordinateSystems.GeographicCoordinateSystem.WGS84

And I got and lat/log coordinate but with wrong location in Google Map.

I already spent a lot of time and I'm stuck on it. Could anyone give me some help? I really can't figure out what is wrong and why it didn't work.

I am using the Proj.Net (http://projnet.codeplex.com/) to apply the projection.

I am using the wkt (that I think is compatible with Google Maps coordinates, correct me if I am wrong) provided in this article (https://alastaira.wordpress.com/2011/01/23/the-google-maps-bing-maps-spherical-mercator-projection/)

I made some changes in code, I changed proj.net to DotSpatial that it is based in Proj.4 and it was recently updated (Thu Apr 28, 2016 at 8:00 AM).

So the code is the follow:

Shapefile nycBoroughs = Shapefile.OpenFile(filename);
nycBoroughs.Reproject(KnownCoordinateSystems.Geographic.Europe.ETRS1989);

and now I got the coordinates as you said:
coordenate = {(-8,44768822459892, 40,0546873641873)}

I also apply the formulas and I got the coordinates also as you said:
x = -940392.35141178139, y = 4873892.4828391289

But now I can't use this coordinates directly in Google Maps, I need to convert again to ETRS1989 to get the coordinates in degree and use it with google maps?

Sorry, probably it is a simple question, but I am a programmer and i don't have enough know-how in system coordinates to understand it easily.


texts: Importing KML file into Google Earth Pro says Junk after document element?

texts:
I'm very new to using GIS platforms.

I've sourced some spatial data for a archaeological project from the 'Portable Antiquities Scheme' website. There was a KML file available to download.

I'm not really familiar with KML file formats at all. From my understanding they don't import easily into QGIS (which is the GIS platform I'm using for my project), so you have to open them in Google Earth Pro and then save as a shapefile which you can then import into QGIS?

I've encountered a problem though in that I can't bring the KML file into Google Earth. I get a error a message which states:

  ... Junk after document element.

From a general google search I understand there might be something wrong with the formatting of the KML file possibly?


texts: Enrich a featurecollection with more data

texts:
I have generated a vector feature collection which is great but I would like to add to this data before saving it out as a CSV. So the data is areas of forest loss in the county of Gwent and I want to add extra fields like average, min and max elevation to each polygon.

How would I do that before I save out the CSV?

var fromList = ee.FeatureCollection("FAO/GAUL/2015/level2");

var filter = ee.Filter.inList('ADM2_NAME', ['Gwent']);
var filteredArea = fromList.filter(filter);

Map.setCenter(1.1, 51.5, 8);
Map.addLayer(filteredArea, {}, 'Regions');

var gfc2018= ee.Image('UMD/hansen/global_forest_change_2018_v1_6').clip(filteredArea);

var lossImage = gfc2018.select(['lossyear']).gte(7);
var lossBand = lossImage.and(gfc2018.select(['lossyear']).lte(20));
var zones = lossBand.updateMask(lossBand.neq(0));

// Add the loss layer in red.
Map.addLayer(lossImage.updateMask(lossBand),
            {palette: ['FF0000']}, 'Loss 2007 onwards');

var area_threshold = 300; / about 100 hectares?

// Define whole county as a region
var geom = filteredArea.geometry();

// Convert the zones of the forest loss areas to vectors.
var vectors = zones.addBands(zones).reduceToVectors({
  geometry: geom,
  crs: lossBand.projection(),
  scale: 30,
  geometryType: 'polygon',
  eightConnected: true,
  labelProperty: 'zone',
  reducer: ee.Reducer.count()
}).filterMetadata("count","greater_than", area_threshold);

// Make a display image for the vectors, add it to the map.
var display = ee.Image(0).updateMask(0).paint(vectors, '000000', 1);
Map.addLayer(display, {palette: '000000'}, 'vectors');

// Export the FeatureCollection to a CSV file.
Export.table.toDrive({
  collection: vectors,
  description: 'Forest',
  fileFormat: 'CSV'
});


texts: PostGIS layers as WMS service with QGIS Server

texts:
I have QGIS, PostGIS and QGIS Server running under Windows 7.

When I make a QGIS project with different Postgis layers and save the project under my QGIS Server path to publish it as an WMS service, i can see no layers at all.
The qgis2web plugin for example also shows me only the map legend but no map at all.
My WMS connection has also only the project name in it, no layers to select.
My KBS is set correct, thats not the problem i think.

Is it not possible to make a QGIS Project with PostGIS data in it, save it, and then publish it as a WMS service?


texts: How to rotate tiles in GeoServer ImageMosaic

texts:
I am starting a project relatively "fresh" to GeoServer. We have a large number of tiled (seabed) images to load up to GeoServer using the ImageMosaic plugin in order to make a seabed mosaic. We know the geolocation of all the images as well as the scaling. We can calculate the required rotation offsets.

When I add a rotation factor to the world file or GeoTIFF header I expect the images to be rendered rotated or "warped" but do not seem to be able to get anything other than unrotated images.

How should I prepare my image mosaic folder so that all of the images are rendered rotated? Do I need to rotate them myself before uploading to GeoServer or can I tell GeoServer/ImageMosaic to do this for me?


texts: How do I find the points of intersection between a small circle and a great circle?

texts:
I only actually need to intersect a great circle arc i.e. the shortest path between two points, with a given small circle.

I found this tutorial for intersecting two great circle arcs, but unfortunately it doesn't cover small circles.


texts: osm2pgrouting --> option '--conf' cannot be specified more than once

texts:
I'm attempting to import some OSM files using osm2pgrouting via command line on a windows 10.

When I specify the mapconfig.xml with the -c switch, I am getting the following error:

option '--conf' cannot be specified more than once error

If I attempt the same without the -c switch, I'm getting:

Opening configuration file: lean
Parsing configuration

Error opening lean: No such file or directory
Failed to open / parse config file lean

These are the steps I've taken so far:

1) Confirmed the mapconfig.xml file was installed under the PostgreSQL/9.5/share directory.

2) Uninstalled and reinstalled postgis 2.2 via the EnterpriseDB Stack Builder.

3) Placed a copy of mapconfig.xml in the same directory as the osm files.

4) Attempted to use just the file and the file path to PostgreSQL/9.5/share.

5) Confirmed the correct pgrouting files were placed under the bin and lib directories.

6) Confirmed that the Postgis and pgrouting extension were created in the DB.

The one item I've not been able to find is a config file for osm2pgrouting to change the default config directory.

If I'm correct, if you specify a file with the -c switch, the default should be ignored, which in this case, it appears that both are being called.


texts: Adding F_ZLEV and T_ZLEV to road network - do I make all ground roads 0?

texts:
I am trying to create a network dataset from TIGER for a very small city. Of course there is no elevation/Z data. I am watching a tutorial that shows how to add Z values (1 for bridge, 0 for ground, -1 for tunnel). They go in detail on how to treat highway ramps, but they do not say if I should do a calculate field for all "ground" roads to set them to 0. Currently they are blank/null. Using ArcMap 10.4 Network Analyst.
TLDR; should I "calculate field" for all of my ground roads to be equal to 0?


texts: Esri file geodatabase data migrated to PostGIS through FME becomes non-editable in QGIS

texts:
I loaded a layer (point) from Esri file geodatabase to PostGIS (PostgreSQL) using FME and when I load it in QGIS, I cannot edit the layer.

Am I missing any step in FME data migration process?

I used the 'drop&create' table handling option in FME.


texts: Is it possible to combine two input types into one ArcGIS Geoprocessing Tool for web?

texts:
I'm working on a project where I basically need to combine the different types of input parameters into one geoprocessing widget, rather than have multiple widgets for the different input types.  Is this possible?

For example, one widget has an input where the user simply clicks a point on a map and outputs a single circle buffer.  another has an input where the user selects the input from a dropdown box and outputs a different buffer.


texts: What is the name of the popular buffer algorithm?

texts:
There is a well-known algorithm that I believe is used for the "buffer" function in some GIS libraries and applications. I have a vivid memory of reading the paper that describes the algorithm and illustrates how it works, but now I am quite embarrassed: I seem to have misplaced the PDF, and I can't remember the name of the algorithm either. The paper was very interesting and I was hoping to read it again.

Can anyone help with a suggestion of what the algorithm is called that I may be (poorly) remembering? I checked the source code for a couple of libraries but didn't see a reference to the algorithm name.

Of course there is probably more than one such algorithm, but I will recognize the one I was thinking of when I see the name and look at the paper again.


texts: Is it possible to show Leaflet Marker Cluster with JSON API?

texts:
I have a problem with displaying Leaflet Marker Cluster working with Leaflet-layerJSON.
There are markers, but I cannot show the clusters. This works with GeoJSON, but now I have different data structure and I'm not sure is it possible to show clusters with JSON (not GeoJSON) data structure:

[{
    "latitude":"51.5615",
    "longitude":"-0.0731"
 },

 {
    "latitude":"51.545",
    "longitude":"-0.070"
 },
 {
    "latitude":"51.4638",
    "longitude":"-0.1677"
 }]

Please take a look at the other (not working) example without angular leaflet directive.


texts: Tools for (semi) automated woodland tracing from satelite images

texts:
I'm trying to draw woodlands in OpenStreetMap, using Yahoo sattelite images.

JOSM editor has some plugins that try to automate the process - you have to click inside the area and the plugin finds the boundries. But the quality is pretty bad.

I'm looking for some libraries/algorithms to get good quality boundaries.

The images i'm working with look like this: http://maps.yahoo.com/#mvt=s&lat=56.907056&lon=24.597595&zoom=14


texts: Removing Trailing Zeros from Numeric Field in ArcMap

texts:
Working with a feature class in an enterprise GDB (SQL Server) which has a numeric field with type Double - precision 10, scale 3 - but is displaying trailing zeros on a lot of the values e.g. 265.360000.

What can be done to remove the trailing zeros from the field?

Edit: The issue appears to only be in the display of the values within ArcMap and not how the values are being stored in the database.


texts: Multi choice parameter data type in ArcMap Toolbox

texts:
I'm a beginner ArcMap user.

I created a toolbox in ArcMap/mytoolboxes which will run an specific Python script from my hard drive.

In my python script I need user to choose between 2 options. I want to add a 2 choice radio buttons in running toolbox page. If user choose one of those radio buttons, the other will be deselected and vice-versa. The problem is that I do not know which data type for this parameter is right for me.


texts: Coordinate Reference System with units in meters for central Germany?

texts:
I found several questions in regard of reprojecting data into a CRS with unit of meters and what CRS to use. Usually the counter question is, "for what region, on what scale" and there the discussions stops. I have tried finding a good CRS myself, but as I am no expert on these things, this whole CRS business is extremely confusing.

So, I am looking for a good CRS with units in meters for:
Central Germany, a town called Paderborn. The Region of interest is quite small, ca. 15km around Paderborn.

I create noise simulations for wind turbine parks. So, I have point data (currently in EPSG:3857) for wind turbines. Then, in a python script, I create a equidistant grid around those points, calculate the noise levels on those grid points, use contour-methods to create "isophones" (lines with identical noise level, e.g. 45dBA) and save this data as line or polygon data.

Generally, this all works very nice and gives beautiful results.

But the really dirty part is this:
In my script I take the EPSG:3857 data and transform it into EPSG:4326. Then I use a self coded haversine transform to use the already transformed data into "my grid" with an origin I choose.
I run the simulation, do the analysis, use again a self coded inverse haversine transform to EPSG:4326 and the transform again back to EPSG:3857.

That probably is not a good idea... Any recommendations?


texts: Display ArcGIS Online Hosted Feature Layer in Google Earth Pro?

texts:
Disclaimer: I am much less familiar with network/server/sharing nuances (much more experienced with desktop work).

I work for a nonprofit and am trying to map out where work is done. I have a large FGDB/HFL (hosted by ESRI in ArcGIS online) with a bunch of polygons of river basins (where the work is done). A lot of people we work with don't have ESRI accounts (and this will soon be a public layer).

My goal: have the polygons display in Google Earth Pro. While a conversion to a KML would be simple, the work territory is constantly changing, so I really need a link to share so that when I make updates, everyone can have the most updated version. I know Google Earth Pro has gotten more flexible in what it can handle and have tried adding the service URL for the HFL into the Add>Network Link option. Also tried creating a WFS. Neither working.

Any easy way to create a KML that will regularly update off a FGDB/HFL?


texts:
I'm familiar with using gdal_calc.py to sum individual raster bands. However, I am curious whether there is an easy way to programmatically sum all bands of a raster into a single value. I have a list of multi-band .vrts, each with a variable number of bands (FYI, each .vrt represents harvested area of a number of different crops in different groups such as cereals, fruits, etc. So each .vrt may have differing numbers of bands depending on how many crops are in that group.) For each individual .vrt, I want to sum up all the bands into a single value and write it to another .vrt. What is the recommended syntax for this?
This question is sort of a duplicate of this. But the linked question is >6 years old so I was hoping to give some new attention to the topic, in case there is a new and better answer.


texts: Implementating Analytical Hierarchy Process with ArcGIS 10.2 for Desktop?

texts:
I want to use the AHP in arcmap 10.2. I have tried using the arcscript available from

http://arcscripts.esri.com/details.asp?dbid=13764

Which was designed for Arcmap 9.2.

I was able to get as far as selecting the input rasters, but when I attempted to proceed with the next step I was met with the following error message:

  Run-time error '339':

  Component 'COMDLG32.OCX' or one of its dependencies is no correctly
  registered: a file is missing or invalid

Ideally, I would like to know how to get this script to work in arcmap 10.2.


texts:
I am using google map in leaflet like this

            var osmUrl = 'http://{s}.google.com/vt/lyrs=s&x={x}&y={y}&z={z}',
        osmAttrib = '&copy; OpenStreetMap contributors',
        osm = L.tileLayer(osmUrl, {maxZoom: 18, subdomains:['mt0','mt1','mt2','mt3'],attribution: osmAttrib}),
        leafletMap = new L.Map('mapDiv', {layers: [osm], center: new L.LatLng(20.6266667,78.749305), zoom: 18 });

i want to know how i can download that tile .so that i can genrate offline map
is there any way for that


texts: Does anyone knows if there is a flat file database of points of interest in a json format somewhere available for download? If not how would one set off to create such a database? Are there any tools available for such a task?


texts: Which preprocessing comes first for pansharpened imagery- atmospheric correction or orthorectification?

texts: For pansharpened imagery, what preprocessing method to do first- is it atmospheric correction or orthorectification? From my readings, it seems the former but the rational polynomial coefficients or RPCs needed for orthorectification are then removed in the corrected imagery which makes it impossible to do the latter. If the reverse is done with orthorectification first, wavelengths needed for atmospheric correction are removed as well making the latter impossible to accomplish. What gives?


texts: Surface Temperature + Air Temperature in Landsat image comparisons

texts:
Quick question about relating Landsat Surface Temperature to locally recorded Air Temperature on a given day:

Are differences in local air temperature captured by the Landsat ARD ST imaging process?

Would I need to adjust for those differences in order to compare Surface Temperature data over time?

In the images below I have visualized surface temp. data from before and after a development with a high requirement for green roofs/ infrastructure (2006 and 2018). The locally measured air temperatures have a difference of 5 degrees celsius on those days, so I boosted the pixel values for the lower range by the difference in temperature. Although this skews the absolute data for the image, my thinking is that it should theoretically allow for a better visualization of spatially relative ST differences. Was this step necessary? Or is it the nature of ST imaging that differences in local air temperature are not necessarily captured in the process?


texts: QGIS2web layer search for more than 1 layer

texts:
I would like to expand my layer search on another 2 or 3 layers.
In the qgis2web.js file, I found the following section:
    var searchLayer = new SearchLayer({
     layer: lyr_Area_5B_Chamber_14,
     colName: 'Chamber No',
     zoom: 10,
     collapsed: true,
     map: map
    });

which corresponds to 1 layer only.
I would like to add another layer.
The thread here is for Leaflet, but states about the Layergroup option:
GeoJSON search box for multiple layers
but if I apply to my example it doesn't work
  var searchLayers = L.layerGroup([lyr_Area_5A_Pole_13, lyr_Area_5B_Chamber_14])

  var searchLayer = new SearchLayer({
  layer: searchLayers,
  colName: 'Chamber No', 'Pole No.',
  zoom: 10,
  collapsed: true,
  map: map
  });

as well as in the following case:
  var searchLayer = new SearchLayer({
  layer: lyr_Area_5A_Pole_13, lyr_Area_5B_Chamber_14,
  colName: 'Chamber No', 'Pole No.',
  zoom: 10,
  collapsed: true,
  map: map
  });

Is there any way to expand the layer search in QGIS2web plugin?
UPDATE:
I tried sth like this
  var searchLayer = new SearchLayer({
  layer: lyr_Area_5B_Chamber_14 || lyr_Area_5A_Pole_13,
  colName: 'Chamber No' || 'Pole No.',
  zoom: 10,
  collapsed: true,
  map: map
  });

with the OR operator. Unfortunately only 1 layer comes


texts:
Is there any to order the Foot Print fields created in Raster Mosaic Dataset?

I know, we can edit, or reorder fields of any FeatureClass, or Shapefile using X-Ray ArcCatalog, ArcGIS Diagrammer.,

But I am looking for reordering fields of FootPrint creating in Mosaic Dataset. Like, I particularly want to add field of default field created during MosaicDataset..

Please refer the below figure. The Field "Positional Accuracy" have been added manually, which used be placed in the end of all fields. Raster(Field Order 1), and Name (Field Order2) Fields created automatically, while we create Raster Mosaic Dataset. I need to reorder the Positional Accuracy field (Field Order2), and Name field (Field Order 3)


texts: Looking for Python module to do lat lon calculations

texts:
I'm looking for a Python module that can provide some / all of the following services using latitude/longitude to define points:

Distance between two points
Point in polygon
Distance along a multi-point path
Identifier of closest point from a dict of points with identifiers
Distance from a point to each point in a dict containing points with identifiers
etc

You get the idea...  Any suggestions?


texts:
I'm looking for a complete, freely available file showing UTM Grid zones. Preferably not a simplified version, but one that shows the full grid zones and all 'exceptions' (e.g. in Northern Europe) - as per http://whatutmzoneamiin.blogspot.com/

Does anyone know where I can download such a file?


texts:
Is it possible to hide icons in the QGIS print composer?

I have used the Customization page in QGIS to create an ini file to simplify the QGIS interface for my users, but I can't see how to do the same for Print Composer.

I know I can turn off entire toolbars, but I need to keep a few icons to allow printing, zooming and changing text.


texts: Adding features with attribution to Feature Service using Jupyter Notebooks and ArcGIS Python API?

texts:
I want to have a user add points by clicking on a map in a Jupyter Notebook and those points get pushed to a feature service.

I'm able to accomplish that task, but I would like the user to be able to customize the attribution for the points they just added.

I'm envisioning a pop up with the various fields that appears when the user adds a point on the map and they just type in a string. I've tried using input(), as well as ipywidgets.

Is there another way to accomplish this?

Most of the examples I see online are for the ArcGIS JavaScript API.

def create_feature(_map, location):

    try:
        object_id = 1
        point = geometry.Point(location)
        feature = features.Feature(
            geometry=point,
            attributes={
                'OBJECTID': object_id,
                'PARK_NAME': 'My Park',
                'TRL_NAME': 'Foobar Trail',
                'ELEV_FT': '5000'
            }
        )

        trailheads_layer.edit_features(adds=[feature])
        _map.draw(point)

    except Exception as e:
        print("Couldn't create the feature. {}".format(str(e)))

m.on_click(create_feature)