texts: Earth Engine Python API, Extracting NDVI values from an Image

texts:
I am using the Earth Engine Python API and attempting to extract date and NDVI values from an image collection. I have added an NDVI band to each image in the collection but I am having trouble figuring out how to extract the actual NDVI value from each band in order to produce a dataframe of an NDVI time series.

I am creating an NDVI band with the code below

    nir = img.select('B5')
    red = img.select('B4')
    ndvi = nir.subtract(red).divide(nir.add(red)).rename('NDVI')
    print('NDVI value:', ndvi.select('NDVI').getInfo())

and I am receiving info from the NDVI band in the form of

NDVI value: {'type': 'Image', 'bands': [{'id': 'NDVI', 'data_type': {'type': 'PixelType', 'precision': 'float'}, 'dimensions': [492, 492], 'origin': [2368, 3104], 'crs': 'EPSG:32612', 'crs_transform': [30, 0, 175185, 0, -30, 5063115]}]}

The documentation on accessing computed NDVI values is limited and I have been scanning tutorials and other exchange posts for ee api related solutions but have found nothing. Any help would be greatly appreciated!


texts: Derived field dropped from layer after restarting session in QGIS

texts:
This is a question about QGIS: I am working locally (no connection to servers or databases) with some maps data which I have imported as a non-visible CSV layer. These were merged against some census data and a derived variable (proportion) was created using the field calculator. I used this to shade the census tracts in a choropleth. I saved my session and called it a day.

When I came back to reload my session, all my previous layers still existed, and the geocoded points remained. However, the choropleth shadings were absent. I opened up the field calculator to look at those data and did not see my derived field. When I tried to regenerate this value, the field calculator throws an error that the field name is not unique. But still I cannot find the field.

Is there a way for me to recover this field?


texts: How to get the projection of ui.Map in Google Earth Engine?

texts:
I know it's 3857, but I would really like to not hard-code it, only so that the script fails later when the GEE team decides to switch.
What I tried is to get the projection from a displayed image of a layer:

Map.layers().get(0).getEeObject().projection()

But it just returns 4326, which is obviously not the map's, nor the displayed layer's projection. It's actually not even the image's projection, but that's another story.

Any other suggestions?


texts: Modelling two-way single-track/one-lane width road in ArcGIS Network Analyst?

texts:
I don't know how make one-lane roads(only room for one car at a time) for one car waiting for another car to pass by in ArcGIS.

Yesterday, I asked about the wait time when the two cars meet in the road, one car should wait for the other car to pass that as Calculating wait time on one-lane road with two-way travel using ArcGIS Network Analyst?

But I realized that I can't even make a one-lane road...

How can I make this road in ArcGIS Network analyst (especially, VRP (Vehicle Routing Problem))?


texts: Creating Polyline Buffer in OpenLayers 2?

texts: I was wondering if there was a way to create a buffer around a polyline feature in OpenLayers 2. I have created a tool that allows the user to draw a circle buffer zone of whatever size and style they want around a single point (using this handy walkthough ), but was just wondering if I could do the same with Polylines. The user would define the size of their buffer zone and then draw their polyline on to the map. Once finished drawing the buffer would be added to the map, like the example I made in MapInfo below.


texts: Why does my "Export to CAD" script work in Python 3.7, but not in 2.7?

texts:
I've written a script that iterates through all feature classes in a dataset in a geodatabase and writes them to separate dwg-files. In this case 230 feature classes. If I run the script in Python 2.7.18 it only converts some of the files (and the amount seems to be random depending on during which day/session I run it). The other files get the following error:

Failed to execute. Parameters are not valid.                           ERROR 000732:
Input Features: Dataset xxxxx does not exist or is not
supported

I tried everything I could think of, and find on the internet, but nothing helped. It only started working correctly when I changed the script to Python 3.7.9. I've run the 2.7.18 script from IDLE and 3.7.9 from IDLE (ArcGIS Pro), if that matters.
I based my script on the ArcGIS desktop "Export to CAD" example code. This is the code I ran (Python 2.7.18):
# Import system modules
import arcpy
from arcpy import env

# Set workspace
env.workspace = r'C:\\test\\test.gdb\\dataset'

# Set local variables
output_type = "DWG_R2018"
output_filepath = r'C:\\test\\CAD_\\'

for feature in arcpy.ListFeatureClasses():
    print feature
    try:
        # Process: Export to CAD
        in_features = feature
        output_file = output_filepath + feature +".dwg"
        arcpy.ExportCAD_conversion(in_features, output_type, output_file, "USE_FILENAMES_IN_TABLES", "OVERWRITE_EXISTING_FILES", "")
        print arcpy.GetMessages()
    except:
       # If an error occurs while running a tool print the message
        print arcpy.GetMessages()

The only thing I changed when running it in 3.7.9 was the print syntax, i.e. use parenthesis: "print (feature)".
This drives me nuts. Why does this code work in Python 3.7.9, but only partially in 2.7.18?


texts: OpenLayers 4.6.5 addLayer circular structure & layer.getLayerArray is not a function

texts:
Currently working with OL 4.6.5 and Angular 6 and for the life of me have read the documentation and API references on the OL site http://openlayers.org/en/v4.6.5/doc/

And despite following there examples of how to add layers such as

  var osmLayer = new ol.layer.Tile({source: osmSource});
  map.addLayer(osmLayer);

  @Input( 'coreMap' ) coreMap;
  vector : VectorLayer;
  source;
  ngOnInit(){
    console.log(this.coreMap);
    console.log(this.coreMap.layerManager);```

      var vector = new ol.layer.Vector({
        source: new ol.layer.Vector(),
        style: new ol.style.Style({
          fill: new ol.style.Fill({
            color: 'rgba(255, 255, 255, 0.2)'
          }),
          stroke: new ol.style.Stroke({
            color: '#ffcc33',
            width: 2
          }),
          image: new ol.style.Circle({
            radius: 7,
            fill: new ol.style.Fill({
              color: '#ffcc33'
            })
          })
        })
      });
      this.vector = vector;
    }
  addLayer(){
    this.coreMap.map.addLayer(this.vector)
  }

The coreMap being the base level of the map (It was created by someone else on the project so I'm picking up from how they have set it up)

I can see elsewhere they have written it like so

this.coreMap.layerManager.addLayer( this.drawLayer )

Which I've also tried and yielding the same result of circular structure error & layer.getLayerArray is not a function.


texts: Why does inline variable output overwrite occurs with Reorder fields tool?

texts:
I am trying to use this Reorder fields modelbuilder tool but for some reason the inline variable (%name%) for the output gets overwritten with each file pass until the last one. Otherwise it works fine, however, I need to process hundreds of shapefiles and rearrange the field order.

Any help or corrections would be greatly appreciated.

I've attached an image from the modelbuilder I am currently using. Here's the script:

Galindo

import os,string,sys

try:
    import arcgisscripting
    gp = arcgisscripting.create()
except:
    import win32com.client
    gp = win32com.client.Dispatch("esriGeoprocessing.GpDispatch.1")

FC_In = sys.argv[1]
NewOrder = sys.argv[2]
FC_Out_WS = sys.argv[3]
FC_Out_Name = sys.argv[4]

strFields = gp.listfields(FC_In)
strNewOrder = NewOrder
Layers = [""]

# Get Field Collection into List
Fields = gp.listfields(FC_In)

NewOrder = string.split(strNewOrder,";")
print Fields
print NewOrder

x=0 #NewOrder enumerator
y=0 #Layers enumerator (first)

# Cycles through the input feature class creating field collections
while x<len(NewOrder):
    Fields.Reset()
    Field = Fields.Next()
    while Field:
        FieldName = str(Field.name)
        if x<len(NewOrder) and NewOrder[x].lower() == FieldName.lower():
            Layers[y] = Layers[y] + FieldName + " " + NewOrder[x] + " VISIBLE"
            x=x+1
        else:
            Layers[y] = Layers[y] + FieldName + " " + FieldName + " HIDDEN"

        Field = Fields.Next()

        if Field:
            Layers[y] = Layers[y] + "; "
            FieldName = str(Field.name)
        else:
            print Layers[y]
            Layers.append("")
            y=y+1
            Layers[y]=""

# Use Field collections to create Layers with proper order until all fields have been used
z=0 #Layers enumerator (second)
LayerList = ""
for Layer in Layers:
    if Layer <> "":
        if LayerList <> "":
            LayerList=LayerList + ";"
        LayerList=LayerList + "Layer" + str(z)
        print LayerList
        gp.MakeFeatureLayer(FC_In, "Layer" + str(z),"#","#",Layer)
        gp.addmessage(gp.getmessages()+ "\n")
        z=z+1

#Collect Info about Input Feature Class
dsc = gp.describe(FC_In)
fc_in_stype = dsc.shapetype
fc_in_sr = dsc.spatialreference

# Create new empty Feature class with desired fields
try:
    gp.CreateFeatureClass(FC_Out_WS,FC_Out_Name,fc_in_stype,LayerList,"#","#",fc_in_sr)
    gp.addmessage(gp.getmessages()+ "\n")
except:
    gp.addmessage(gp.getmessages()+ "\n")

# Append Records to new feature class
print FC_Out_WS + "\\" + FC_Out_Name
print FC_In
gp.append_management(FC_In, FC_Out_WS + "\\" + FC_Out_Name, "NO_TEST")
gp.addmessage(gp.getmessages()+ "\n")

del gp


texts:
I'm trying to run a console script, but it gives an error:
:

  AttributeError:'NoneType' object has no attribute 'identify'

Yet it works well in QGIS 3.7

Here is the script:

import sys
import time
import os
import PyQt5
import numpy

# Qgis modules
from osgeo import gdal
from qgis.core import QgsRasterLayer

#from qgis.core import *
import qgis.utils
from PyQt5.QtCore import QFileInfo, QSettings
from qgis.gui import *
from qgis.utils import *
from PyQt5.QtCore import *
from PyQt5.QtGui import *
from qgis import *
from PyQt5 import *
from qgis.core import QgsRasterDataProvider

chemin = 'D:/Projet_Python_Soil/tests/'

for element in os.listdir(chemin):
    if element.endswith('.tif'):
        fileInfo = QFileInfo(element)
        baseName = fileInfo.baseName()
        rlayer = QgsRasterLayer(chemin + element, baseName)
        print (baseName)
        ident = rlayer.dataProvider().identify(QgsPoint(368371, -406601), QgsRaster.IdentifyFormatValue)
        #os.system("pause")
        print (ident.results())
        if ident.isValid():
            iden = {}
            iden = ident.results()
            print (iden[1])


texts: OpenStreetMap pbf files- extracts all street/roads intersections in a bounding box

texts:
I wrote a query which extracts all street/roads intersections in a bounding box from OpenStreetMaps, as can be seen here:
    <osm-script>
<query type="way" into="hw">
  <has-kv k="highway"/>
  <has-kv k="highway" modv="not" regv="footway|cycleway|path|service|track"/>
   <bbox-query e="51.0" n="50.748" s="50.746" w="7.154"/>
</query>

<foreach from="hw" into="w">
  <recurse from="w" type="way-node" into="ns"/>
  <recurse from="ns" type="node-way" into="w2"/>
  <query type="way" into="w2">
    <item set="w2"/>
    <has-kv k="highway"/>
    <has-kv k="highway" modv="not" regv="footway|cycleway|path|service|track"/>
  </query>
  <difference into="wd">
    <item set="w2"/>
    <item set="w"/>
  </difference>
  <recurse from="wd" type="way-node" into="n2"/>
  <recurse from="w"  type="way-node" into="n3"/>
  <query type="node">
    <item set="n2"/>
    <item set="n3"/>
  </query>
  <print/>
</foreach>
  </osm-script>

However, I want to work locally on PBF files so I need to figure out how to do it. is it possible to extract all street/roads intersections in a bounding box using OSMFILTER/OSMCONVERT/Osmosis/OTHER?


texts:
I have windows 7 in Arabic. Now I dowloaded QGIS 2.14.3-2 in French.
But when I installed it shows Arabic menu/tools...
How can I change the menu language to French?


texts: How do I use a layer of Bing / Google a "WebApplication" in "ArcGIS Server Manager"?

texts:
I'm using "ArcGIS Server."
I want to publish a "WebApplication".
I created a map of my city.
I want to have as the background world map with satellite images of "BingMaps" or "GoogleMaps" or "ArcGIS World Map free".
See the image below. I tried it. But it bring to me this error: "Unable to add layer. Request failed with HTTP status 417: Expectation Failed."

https://sites.google.com/site/christiancrca/erro%20BING%20maps%20Arcgis%20Christian%20Rogerio.jpg

Any idea about how to use as background a service of image satellite maps with my map?


texts:
Using the Qgis atlas functions im having trouble with the atlas filtering out the features that do not show up in the map extent.

For example i want the full color ramp to show up, even
though there are some features that are not within the atlas map extent.

Below is what the atlas is currently doing:

This is what i need it to look like

Everything is set up and working in the atlas coverage layer to dictate the layers that should be on and off. The legend needs to change with the atlas feature, but not limit the features displayed

I feel like there is some options in the Filter legend by expression options, but im not sure what expression can get me to the result i need.


texts:
I successfully imported a 7 gb raster into my PostGIS database using the following script
raster2pgsql.exe -I -c -d -t 100x100  -s  <SRID> <file_path> public.<table_name> | <psql_connexion>
The problem I encountered is when I try to load the raster in QGIS it takes a lifetime and sometimes it freezes.
Do you have any idea how can I load my raster quickly?
Do I have to change something in the loading parameters specially the -t?


texts: Improving visualization of streets intersections in ArcGIS

texts:
I need to improve the street network intersections.
The image describes what I want.
These streets are all on the same level (elevation).


texts: Choosing Topographic Wetness Index or Flow Accumulation raster?

texts:
What is the better raster to use qualitatively - Topographic Wetness Index or Flow Accumulation - to clearly see potential areas of surface flow/wetness across a landscape?

As an example, what's "better" for the farmer who'd like to see potential flow of water in his fields, under wet climatic conditions.

The rasters were created from the same base LiDAR dataset, but TWI was interpolated to a cell size about 3.2 ft2 and the Flow Acc raster cell size is 1 m2. TWI and Flow Acc Rasters were not generated from the same DEM.

A couple pictures in different landscape types. Yellow equals TWI raster. The orange lines represent flow accumulation thresholded at > 1 acre and and blue represents thresholded raster of 0.25-1 acre:

Urban Landscape:

Terraced and Non-Terraced Agricultural Field:

My current line of thinking is TWI has more practical use for the lay person, as opposed to a Flow Accumulation layer that might be more useful in a hydrologic model. But perhaps I'm not approaching this correctly, as TWI and Flow Acc are not exactly apples to apples raster products.

Ultimately, I want to make sure I'm understanding their differences and end uses clearly before making a suggestion.

EDIT:
To investigate Jeffrey's comments about contour bias, I thought it would be helpful to share a picture that compares the TWI raster that results from using the DEM that was generated from LiDAR data processed by our internal IDW/smoothing process versus using the DEM that was generated from the same base LiDAR data but generated as part of the final production of elevation data on PA's spatial data portal.  NOTE: The PA Spatial Database DEM only covered a small watershed of interest (for faster processing times), while our internal process generated TWI from a DEM covering an entire county.

Yellow = TWI from PA spatial database DEM
Red = TWI from our internally generated DEM.


texts: ArcGIS Desktop cannot display TIFF image manipulated with GIMP (graphic software)?

texts:
I used the graphic software program GIMP to manipulate a geoTIFF image (GIMP cannot handle geoTIFFs, so after manipulating the image in GIMP the geospatial information is lost).

When I now try to import the manipulated TIFF into ArcGIS 10.3.1, for some strange reason Arc displays the original TIFF image without any of the alterations I made in GIMP (outside ArcGIS, I see the manipulated version though)


texts: positionAlongLine() does not give equal distances on curve using ArcPy?

texts:
I'm trying to use positionAlongLine() to create equally spaced points down a polyline as part of a loop. I also need to calculate the distance at the same time.

However, on curves there is a mismatch between what the distance SHOULD BE and what the distance IS when the points are measured using the measure tool in arcmap.

See example image - 18 points on a 660.979m curved line should all be 36.72m apart but instead vary from 35 to 40m each (but the calculated distance in the script is correct as 36.72m).

I'm not sure why this is not working for curves - does the positionAlongLine() take the curve into account to get the correct position? If not, is there a way of adapting it?

def fillInPoints(prevID):

prevStreetLen = streetLengths[prevID]  # How long was the street before (before we add our point)?  # 18/04/18 - REMOVED THE (*1000)
streetVertexCount = len(xyDict[prevID]) - 1  # How many vertices on street according to the geometry?
distBetween = prevStreetLen /streetVertexCount  # Average distance between vertices

if distBetween > 100:  # 100 is our max distance, fill in between if greater
    streetVertexCount += floor(distBetween / 100)  # add in the min number to bring distbetween down below 100 (e.g. 300m = 3)

thisDirection = streetSeq_dict[streetSeq-1]['direction']  # Direction (FT/TF) of PREVIOUS street

# Here we are filling in the gaps where the distance between the start and end of the segment are too long
if streetVertexCount > 1:  # More than just the start and end (NB we've already added vertices if streets were >100m)

    # Loop through each vertex
    for n in range(1,int(streetVertexCount)):
        position = n * (1 / float(streetVertexCount))  # Get the position (e.g. 0.5 if only 1 point added)

        streetShape = streetShapeDict[prevID]
        thisDistance = 0

        # Subtracting from the current distance as we're looking back to the street before
        if thisDirection == "FT":
            thisPoint = streetShape.positionAlongLine(position,True).firstPoint
            thisDistance = distanceFromStart - ((1 - position) * prevStreetLen)

        else:
            thisPoint = streetShape.positionAlongLine((1 - position),True).firstPoint
            thisDistance = distanceFromStart - ((1 - position) * prevStreetLen)

        outputList.append([thisPoint.X,thisPoint.Y,thisDistance])


texts: How do you get NoData to show up transparent in a .TIF file?

texts:
I am creating a handful of raster images using ArcGIS 10. I've been exporting the rasters as .TIF files. From what I've read, even though you can set the background to No Color, this actually doesn't do anything since .TIF files do not support transparency. My background ends up gray in color for NoData areas when I select this option.

I have tried setting the .TIF file's properties to display NoData as 255,255,255 and that didn't work either. I also tried reclassifying the data and setting the NoData color to transparent and that failed as well. I even tried using QGIS and selecting the pixel for NoData and changed that to 100% transparent and that didn't work either.

My goal is to get the .TIF loaded into TileMill, then export the .TIF to MapBox Tiles and display the image using MapBox's Online Editor. I have done this successfully in the past exporting the image from data view in ArcGIS, but for whatever reason can't seem to get it to work today. Any thoughts or suggestions?


texts: How do I change label color depending on feature intersection in QGIS?

texts:
I have two polygon layers. Layer 1 is the coverage layer for an atlas. I would like the color of labels for layer 2 to change based on whether or not the feature overlaps with the current atlas controlling feature in layer 1. Here is what I've tried as the expression controlling label color in layer 2:

case when intersects_bbox($geometry, atlas_geometry)
then color_rgb(0,0,0) else color_rgb(255,255,255) end

But all the labels are black no matter what I try. What am I doing wrong? (and is there a better way to do this?


texts: Linking to PDFs via attribute field of shapefile for ArcGIS Online?

texts:
I am working with some data that cannot be presented in an attribute table.  Therefore I would like to attach a specific PDF to each point within a shapefile.

Is there a way this can be done in ArcGIS for Desktop and then uploaded to ArcGIS Online so that an interactive map can be published with this information?


texts: Can you show a WMTS layer above largest scale in QGIS?

texts:
When zooming in on a WMTS layer, the layer disappears at a set large scale. Does QGIS have the ability to show the last scale available when zoomed in further?
I can see it expanded at the new scale before it disappears.
For example, the layer is visible at 1:10000. At 1:9000 it is gone. I'd like to see the layer even if it is pixelated. Am I missing an option somewhere?


texts: Not enough memory in the file group "Primary": SDE Database error with MSSQL

texts:
I got this error message whilst trying to load data into a raster in ArcMap 10.0.

Underlying DBMS error [Microsoft SQL Server Native Client 10.0: Eine neue Seite für die sde-Datenbank konnte nicht belegt werden, weil in der Dateigruppe 'PRIMARY' nicht genügend Speicherplatz verfügbar ist. Speicherplatz kann durch Löschen von Objekten in der Dateigruppe, Hinzufügen von Dateien zur Dateigruppe oder Festlegen der automatischen Vergrößerung für vorhandene Dateien in der Dateigruppe gewonnen werden.]
Es ist kein Raumbezug vorhanden.
Fehler beim Ausführen von (CreateRasterDataset).

A rough translation:

" A new page for the sde database couldn´t be allocated, because in the filegroup "PRIMARY" there is not enough disk space. Memory can be gained by deleting objects, adding files to the filegroup oder determining the automatic enlargement for existing Files."

The Database is an sde 9.3 File-Geodatabase which is on my harddrive and not in the MSSQL database.


texts: Which year change does NLCD Land Cover Change Index show?

texts:
After downloading the National Land Cover Dataset (NLCD) Land Cover Change Index from the Multi-Resolution Land Characteristics Consortium (MRLCC) website, it only contains one raster file. However, the website states that the data is from a series of years (2016,2013,2011,2008,2006,2004,2001).

Does that mean that this file indicates the total change from the year 2001 to 2016 or is there a way to view the change between each set of years?


texts: Getting ISpatialReference from ArcObject mapcontrol in VBA

texts:
I am working with ArcObjects and ESRI for the first time, and I am trying to do what is probably a simple task.

I have a map control in VBA that is associated with a mxd.

I have a table with XY data, and XY data is non-spatial, but the XY is in the same coordinate system as the mxd.

I can get the XY data into a ITable, but I cannot seem to get the ISpatialReference from the map control.

I would like to create a symbolized layer from the non-spatial coordinates, and add it to the map control.

Any samples in VBA of how to do this would be very appreciated.

Thanks again for you time and patience.


texts: Unable to publish BuildingSceneLayer from revit file in ArcGIS Pro?

texts:
Unable to publish a Building Scene Layer in ArcGIS Pro version-2.3 application running on Windows-64 bit operating system, from a Revit file.

When I try to do so I get:

  Your portal does not support publishing this dataset type

How to check and modify the dataset types are supported in the portal if possible?

  The revit file was downloaded from example from the Autodesk Revit website example given below:

[http://www.autodesk.com/revit-rac-basic-sample-project-2018-enu?_ga=2.195217900.1145496523.1548998945-1799372949.1548830014][1]

The version of ArcGIS Portal being used is 10.6


texts: Creating buffer of points with pre-defined area?

texts:
I need to create a buffer with a predefined area, my inputs are points, at first I tried it by calculating the area that each buffer may have and then giving the radius as an input, but then I have to dissolve the buffers polygons so the area change.

Is there some tool that can do this?


texts: What does Image width is not in range Valid range is from 1 to 2048 mean?

texts:
We are facing an issue on ESRI Silverlight (v 3.1) and ArcGIS Server 10 SP4 based web application. Few of the users complaining that map is not being displayed in their machine. When I investigated fiddler log, map service request throws 404 error . Below is the detailed

  Error Specified format is not supported  Code: 400 Image width is not in range. Valid range is from 1 to 2048.

Users are using dual monitors and different browser zoom levels, so few people are OK and few others facing trouble in viewing map layers, how to fix this issue for different screen resolutions and browser zoom levels?

I am looking for comprehensive solution for this issue.


texts: How to create a unique identifier in the attribute table based on spatial location within another layer?

texts: I am currently working in ArcGIS 10.2.2 for Desktop and would like to know if there is a workflow for creating a unique ID for polygons within a shapefile based on their location. I am using data driven pages that are clipped to individual Township/Range/Sections of a city. There is polygon information that spans the entire city, but would like to create some type of unique identifier that allows users (of the exported attribute table) to quickly see where a polygon is located based on the identifier. What I currently have in mind is to list the Township/Range/Section and add a numeric suffix. For example, in the image below the first polygon in Township26-Range42-Section27 (highlighted in red) would be 26-42-27-1. The Township/Range/Section information is saved in another shapefile, so I am not sure the best process to append that info to the individual polygons within this other shapefile based on location.


texts: Passing tab file paths into QuickImport?

texts:
I'm just trying to add tab file paths into Quickimport to convert to shapefiles. I generate the tab file paths below but just not sure how to load them into Quickimport using for loop?? No python errors reported.

import arcpy
from arcpy import env
import os

path = "C:/temp1"

     for root, dirs, files in os.walk(path,topdown=False):
      for file in files:
         if file.endswith(".TAB"):
             print(os.path.join(root, file))

#Execute Quick Import
  # arcpy.QuickImport_interop("**How do I insert the tab file paths from the above list into here using for loop?**")

When its run, there's no errors reported. It prints the list of tab file paths like this..

C:/temp\Boundary.TAB
C:/temp\circle.TAB
C:/temp\square.TAB
etc,


texts: How to create a copy of a feature class that is vertically and horizontally offset by a given distance?

texts: I want to duplicate a polygon feature class and offset all of the polygons by about 10 feet in both the x and y directions. Is this possible to do using arcpy, or are there any geoprocessing tools that could help me do this?


texts: Finding close and parallel lines using ArcPy?

texts:
I need help with script for QA of water network. I need to find all lines that are almost parallel and close to each other less then 2 meters.

The script has 3 stages:

add azimuth field and calculate azimuth.
iterate and find close lines.
If lines have almost the same azimuth - select them.

I've used these lines here for first stage and need help with the loop that will go thru all the line and find close pairs.


texts:
I'm attempting to convert a USGS US Topo GeoPDF in Windows.  I'm using GDAL 1.11.1.  My initial attempts were with this Python script, which automates the process by feeding a command to the Windows Command Prompt.  After some initial problems (I couldn't get past line 10), I decided to try directly with the Command Prompt.  I used this USGS resource as a guide, which offers step-by-step instructions for using GDAL in the Command Prompt.  To simplify the process, I have been following the guide literally and using the same Santa Fe topo map.  When I enter the following command:

gdalinfo NM_Santa_Fe_20131108_TM_geo.pdf –mdd LAYERS

I get this error:

gdalinfo failed - unable to open 'NM_Santa_Fe_20131108_TM_geo.pdf'

I found a few different resources that said that this error could originate from the GeoPDF file format not being compiled into the build, and suggested the following command:

gdalinfo --formats

When I run this, the GeoPDF format is indeed listed:

PDF (rw): Geospatial PDF

GDAL functions properly if I read in other formats (PNG, TIF).  Does anyone have experience with this issue?


texts: Postgresql - problem updating row with complex calculation

texts:
A calculation for the replacement value of a tree exists:

Variables:

index_price for species;
circumference of tree;
condition measures: vitality, root_injuries, stem_injuries,
  crown_injuries
(A condition measurement can be between 0 and 4, 0 worst, 4 best).

Calculated inputs:

price_per_cm2 for small tree: index_price/13.45
area_of_tree: (circumference*circumference)/(4 * pi)
condition_adjustment: (vitality + root_injuries + stem_injuries +
  crown_injuries)/16
establishment cost: (area * 70 SEK) + 10000 SEK, or 75000 SEK, whichever
  is smaller

The calculation is:

((Price per cm2 * area) * condition adjustments) + establishment cost.

So, with javascript, this works:

circumference = diamter * pi;
establishment_cost = (area * 70) + 10000 > 75000 ? 75000 : (area * 70) +
  10000;
replacement_price =
    ((index_price/13.45)*((circumference*circumference)/(4*pi)) *
    ((vitality + root_injuries + stem_injuries + crown_injuries)/16)) +
    establishment_cost

With PostGreSQL psql I try:

update trees
set replacement_price =
  ((index_price/13.45)*
  (((diameter*3.141592654)*(diameter*3.141592654))/(4*3.141592654)) *
  ((vitality + root_injuries + stem_injuries + crown_injuries)/16)) +
  establishment_cost;

The query only works for one row, where the diameter is 4cm, in all other rows (diameters of 20cm and greater) replacement_price = establishment_cost, suggesting the start of the calculation is returning 0.

Can any talented mathematicians or SQLers help me find my error?

Eventually, I hope to place the code in a stored procedure, but thought I'd start with manual updates...


texts: Determining if lat/lon coordinates in OSGB36 or WGS84?

texts: If my DEM data's spatial reference is the British National Grid, when I use the xy tool in ArcGIS 10.2, are the lat/lon coordinates in OSGB36 or WGS84?


texts:
I am new in pygeoprocessing and I am trying to do a really simple raster calculation. I would like to calculate element-wise average over several rasters.

I first though pygeoprocessing could handle it the same way R's Raster package does it (stack of rasters and then apply a function) but couldn't find a way to do it so I am trying with the raster calculator.

I have 2 ratsers R1.tif and R2.tif and would like to calculate the raster Rmean.tif with same extent as R1 and R2. Here is what I tried:

import pygeoprocessing.geoprocessing as geop

listraster_uri = [('R1.tif',1),('R2.tif',1)]
rasterout_uri='R.tif'

def avg(r1,r2):
   result_avg=(r1+r2)/2.
   return result_avg

geop.raster_calculator(base_raster_path_band_list=listraster_uri,
            local_op=avg,
            target_raster_path=rasterout_uri,
            datatype_target(np.float32,1),
            nodata_target=-1,
            calc_raster_stats=False)

but I get this error :

TypeError: in method 'Driver_Create', argument 6 of type 'GDALDataType'

because I don't know how to specify the "desired GDAL output type of the target raster"m beside the fact that I don't know if I am doing the right thing here...

Is there a way to stack the ratsers and apply a numpy function like pandas does over dataframes ?


texts: I am trying to import .shp to PgAdmin III but is error because of "dbf file can not be opened" and I dont understand why.


texts: Creating dynamic maps based on updated attributes in SHP file

texts:
I am new to GIS. I can create maps and do basic transformation of the SHP file using QGIS. I also succesfully applied formatting based on attributes. I am able to add new attributes to SHP file in QGIS, and I have published my map with styles to my instance of Geoserver (SLD styling compatibility is painful). It's all working fine.
I've also familarized myself with kartograph.py
note: I am only using free software

My workflow is:
I get SHP file for my country with counties. I then add new CSV layer containing KEY matching the attributes in CSV file with the attributes in SHP file. ( I do this to add some numeric data to Counties). Then I will add some styling to counties based on the new attributes and export to new SHP file.
I then publish this to GeoServer, correct the styling and I am done - I can see my data on the map via webservice.

The question:
I am able to generate CSV with new data regularly (e.g. daily) but adding the updated attributes to shapefile and then publishing to geoserve, correcting the SLD style manualy everyday seems very troublesome. I feel there must be better way.

I am wondering how some of you Experts tackle this issue?

Is my workflow/approach incorrect?


texts: Count points within a given distance of a border

texts: I have a vector layer of counties in Florida and a layer of points. I want to count the number of points that lie within 10km of the boundary for each county pair. This is what this looks like in QGIS -- any ideas on an efficient way to do this?


texts: Remove fill / color but keep border lines

texts: I would like to remove the colors (fill) of each province (shapefile), but keep the border lines, but can't seem to find out to remove fill? I am using QGIS 3.4.


texts: zoom rubber band missing in IE11 with esri js api 3.10

texts:
I am no longer seeing the rubber band in the browser when a user does a zoom operation where they can draw their zoom extent. This started breaking when i added in a X-UA tag = 'Edge'.
I am testing with IE11 but need to support IE10 as well.
The reason i had to change that X-UA tag was because some other thing was breaking in IE.

Anyone having issues with this zoom rubber band?
All of the esri samples now use the +/- zoom buttons and not the variable zoom so i can't even test one of those.


texts: CSVLayer doesn't show HeatMap in Angular 6

texts:
I am using ArcGIS JavaScript 4.8 API and I have integrated the CSVLayer to an Angular 6 application but it is not showing up any heat map.
The legend shows that it is getting renderer data and heatmap is added as renderer. But no heatmap is getting shown.

But when I click on the map where I am expecting the heatmap, I can see the PopupTemplate showing my data. An image is given below:

When I zoom into my CSV layer, I can see the points on the map which I have created by replacing heatmaprenderer with simplerenderer on scale event. Screenshots are given below:

No Heatmap:

After Zooming in on CSV Layer (Top one):

Also, I converted an example from the ArcGIS Website to an Angular component and tried to show it on the website. But it is also not showing up any data.
I used 'Create a scale-dependent visualization text' example from here. Sandbox Code is here.

My code is on GitHub. (https://github.com/SiddharthMishraPersonal/arcgis-csvlayer-example)
In the above images, you can see that the second map is not showing any heat map.

Please let me know if I am doing something wrong.


texts:
I need to develop a prototype of a web application that speaks :

When you click on a city (point) it says the city name (Audio files are ready)

I want to use ArcGIS Javascript API. do you have any idea how can I implement this idea ?


texts: Trouble downloading Geodatabase wetlands from FWS.gov

texts:
Have tried using every browser (chrome, firefox, safari, ie, opera) to download wetland data from:
FWS.gov wetland data download with no luck. It either chokes around 400 MB or finishes up with a corrupted zip file which cannot be opened.
It is only an issue with the larger .gdb files (over 400 MB), like Texas which is 1.7 GB.

Any ideas on how to download these files without error?


texts: Atmospheric correction for pansharpened imagery?

texts: In converting DN values to radiance as well as reflectance, do we use the same equations in raw satellite imagery for pansharpened products? I saw this DigitalGlobe technical note for DN to radiance conversion where "the information applies to all Quickbird products EXCEPT pansharpened" imagery. So may I know then how pansharpened imagery are atmospherically corrected, i.e from DN to radiance to reflectance?


texts:
edit:

This question was about a free service that no longer exists for viewing Sentinel 2 data from AWS before the 'requester pays' condition was put in place.

If you want to find out more on how to access data yourself consider the answer to this question or go to Sentinel Hub forum by Sinergise (the company that ESA pays to manage this data)

https://forum.sentinel-hub.com/c/aws-sentinel

specifically this thread

https://forum.sentinel-hub.com/t/changes-of-the-access-rights-to-l1c-bucket-at-aws-public-datasets-requester-pays/172


texts:
I am very confused about the measurements of ST_Distance in PostGIS. If I measure along the equator I get longer measurements using the SPHEROID instead of the sphere, but only up to 179 degrees (and this makes sense since distance on a spheroid is obviously longer). From 180 degrees + this turns around and SPHERE surpisingly yields longer results

SELECT  ST_Distance(ST_Point(0,0), ST_Point(179,0), True)/1000 AS spheroid_dist,
ST_Distance(ST_Point(0,0), ST_Point(179,0), False)/1000 AS sphere_dist;

outputs

spheroid_distance;sphere_distance
19926.1888494779;19903.9192724991

while

SELECT  ST_Distance(ST_Point(0,0), ST_Point(180,0), True)/1000 AS spheroid_dist,
    ST_Distance(ST_Point(0,0), ST_Point(180,0), False)/1000 AS sphere_dist;

outputs

spheroid_distance;sphere_distance
19903.5933909347;20015.1143522337

To make it more confusing I get different results for the spheroid measurement for different PostGIS versions. Above are results for Version 2.1

Version 2.3 yields the following results (for the second SELECT statement, i.e. 180 degrees):

spheroid_distance;sphere_distance
20003.9314586255;20015.1143522337

So what numbers are to be trusted? Can anyone shed some light?


texts: Creating a GUI toolbox for a GRASS add-on

texts:
In Vaclav Petras' GRASS add-on guide it is stated that GUI toolboxes for GRASS add-ons are defined in the toolboxes.xml file, that must then be added to the main installation folder:

  You have to copy this file toolboxes.xml to the GRASS GIS configuration directory in your home directory (e.g., ~/.grass7/toolboxes/toolboxes.xml on Linux).

This mean that when a user installs an add-on using the g.extension tool a GUI is not made available automatically. As far as I understand, the user must update the toolboxes.xml file her/himself.

Is there any way to automatically register a GUI toolbox in the user's system?


texts:
I want to add database connection to SQL Server from ArcMap.

But, I don't know what to write as 'instance'. I've tried 'NT Service\MSSQLSERVER' but it doesn't work.


texts: ESRI Geometry Union slower than Built in dissolve which is slower than Shapely

texts:
I work with point datasets a lot.  We try to keep a running layer of points and update the dataset with new points while removing the old points.  So, we do a lot of checking the extents of the new points with the extents of old points to make sure we're not overwriting new data or keeping old data.

To do that, one of the methods I run is a buffer of the various points and then dissolve all of the buffers into a single polygon.  I'm wondering if anyone has any ideas how to speed up ESRI without third-party software (IT is hesitant to install any third-party "software" and they don't want to install the shapely, GDAL, and fiona libraries on our virtual server).

I did a recent test with one file that had a little over 11k points in it.  I ran ESRI's built-in Buffer tool under the Analysis toolbox and it took just about 4 minutes to buffer at 120m and dissolve all of the features.  My shapely/fiona implementation took 13 seconds.

So, I modified my code to try to use the ESRI geometry features.  Here's the essentials of my code (I also implemented allowing to group features but that is extraneous to this discussion):

tempFC = "ProjectedInput" #the code either copies a projected input or projects the input

fields = ["SHAPE@"]
inputCursor = arcpy.da.SearchCursor(tempFC, fields)
outputCursor = arcpy.da.InsertCursor(outputFC, fields)

geoms = []
for inputFeature in inputCursor:
    geoms.append(inputFeature[0].buffer(bufferDist))

mergedGeom = geoms[0]
for i in range(1, len(geoms)):
    mergedGeom = mergedGeom.union(geoms[i])

outputCursor.insertRow(mergedGeom)

The union function is incredibly slow.  I started the above buffer 22 minutes ago and it's only half-way through unioning all of the features.  (It turns out I selected SR_01 instead of SN_01 when I ran the above code.)  With the correct dataset, it took a little under 3 minutes.  So, it was faster than the built-in but still not fast enough.

So, it seems there's not going to be any way to get the shapely speeds out of the arcpy library.  Anyone have any ideas or anything that has worked in the past?


texts: ArcGIS -ArcMap Advanced SQL Expression Returning Blank Output

texts: Want to return the comments that only show "Pending Annexation - Pending City Limits", but is returning a blank value.


texts: Is it possible to undo and redo outside edit session in ArcObjects?

texts:
I would like to know if it is possible to undo and redo changes outside of an edit session?

If yes, how can I do it in ArcObjects using .NET?

I am working on an ArcSDE 10.1 versioned geodatabase?


texts: How to prototype GIS workflows in Sketchflow?

texts:
I've started to use Microsoft Expression Blend and SketchFlow to visualize senarios in our line-of-business application prototypes. Right now I use the Sketch Style to do a Rectangle and write a text Map or Layout inside to visualize a Map or a Layout.

Is there any way to get a sketchy looking map without adding dependencys to any ArcGIS APIs. I want it to be light-wight so it's easy to distribute the prototype without getting license problems.


texts:
omnivore.wkt.parse("${feature}").removeFrom(map);

I select features by this line
 omnivore.wkt.parse("${feature}").addTo(map);
but when I remove from map doesn't work

this code works but if I don't know delete feature

<c:if test="${not empty mapBean.layers}"> console.log('Kaydedilen featurelar:');
<c:forEach var="layer" items="${mapBean.layers}"> <c:if test="${layer.id==1}"> <c:if test="${not empty layer.features}"> <c:forEach var="feature" items="${layer.features}"> console.log("${feature}");
omnivore.wkt.parse("${feature}").addTo(map); </c:forEach> </c:if> </c:if> </c:forEach> </c:if>


texts: Choosing Projected Coordinate System for Western Australia?

texts:
I am using QGIS to observe and analyse spatial distribution of earthquakes in Western Australia. I want to measure polygon areas, and create buffers, for which I need to use a projected coordinate system.

What projected coordinate system should I use for Western Australia?


texts:
I would like to open MXD using ArcGIS Runtime WPF.

There is no example about this.

Anybody has a solution?


texts:
I've used older versions of QGIS at an ultra newbie level and never had this problem. I'm on Linux Mint and just installed 3.16, but can't even figure out how to get it to recognize my data folders, or even the drive that the folders are on. Even if I open QGIS from a project or layer on that drive and in those folders, when I go to save or add layers it can't even see the drive those folders are on. It keeps going to my desktop, which is the unfortunate place I first opened a file with it. The Project Home is in the folder I want to access but still QGIS can't see it. I'm stumped.
If I try to add or save a shapefile, the dialog box for browsing to my chosen location won't even show the drive that I need to access.


texts: Using Label Expression with Text Formatting Tags and Ampersand

texts:
I am using a Label Expression to label a parcel layer.  The expression is:

"<FNT name=""Arial"" size=""12"">" & [Owners_Nam] & "</FNT>"

This works correctly for some labels, but any [Owners_Name] field that has an "&" in it, such as between two owners names (Smith, John & Sue), does not show up correctly.  The parcel is labeled as:

<FNT name="Arial" size="12">Smith, John & Sue</FNT>

Is there any way to solve this without removing the "&" from the [Owner Name] field?

I know how to replace the & in Field Calculator; I would just like to preserve the original field text.  I've had this same problem in the past using the <BOL> tag too.


texts: Seeking high resolution real-color (RGB imagery) dataset?

texts: Specifically interested in the the South Florida region, USA. The aim is to pull up high quality aerials to detect manholes and wellheads from space! These things can typically be half a foot to a foot in diameter.


texts: ArcMap - Python: Point in Polygon - within/contains function

texts:
I have one layer with polygons (say, countries) and one layer with x and y coordinates of points (say, cities).

Now, for the city table, I want a new column giving me the id of the country it lies in. I believe the general approach would be to create a new field and then use the field calculator, either of the within or contains function and a loop, like this:

Is City A in Country A? No.

Is City A in Country B? No.

Is City A in Country C? Yes. --> Note down ID.

Is City B in Country A? And so forth.

However, I cannot figure out how to write the code. How do you use these functions if your inputs are from different layers/attribute tables? How to specify the point input since it consists of an x and a y variable? Does anyone have a Python code example for how s/he used the within or contains function?


texts: Showing dialog box even though QGIS is freezing

texts:
I have a long script in PyQGIS that contains three processing tasks after another.
I want to display a pop up box (like "QProgressDialog"-Box) for each task when it is finished.
Example:
When Task 1 (Merge) has finished, I want a box to pop up ("Task 1 (Merge) is finished). Then the Task 2 is executing and when it is finished i want another box to pop up...
Unfortunately only the first "QProgressDialog"-Box is shown and QGIS freezes as it is processing.
How can I achieve to display a pop up box (or even better a progress bar) when QGIS is freezing?

My code is:
bar = QProgressDialog("Operation in progress", "Cancel",0,100)
bar.show()
bar.setRange(0,0)

#the processing tasks are following below
...


texts: Clipping using multiple polygons from same layer in ArcGIS

texts: I have a polygone shapefile and now I want to select a couple of polygons which are islands in one big polygon and clip that big polygon with these selected islands. I know that I can use ERASE tool but that is not a good solution because I must to export features in separate layers and after that to paste it back. So,  I need some tool (clipper) which allow me to select more than one splitting polygons. I want it in ArcGIS 10.4,  but if someone know solution in qgis that is acceptable too.


texts: Calculating the fastest routes and some points on the same time

texts:
I have two drivers and each driver has a set of about 15 locations where they need to drive to on a single day. Calculating the shortest route for both drivers is not the problem (using a matrix routing api).

The drivers can have the same location in their route. And if they are the same, then they both need to be there at the same time. So I need to make software that calculates the fastest route but sometimes the drivers need to be at the same location at the same time.

My question: How can I make this software, and are there any libraries I can use?


texts: How to create a 50x50 m hexagonal grid with mmqgis?

texts:
I'm trying to use the mmqgis plugin to create a hexagonal grid, and I am really perplexed on how to fill in the variables required. I would like to have hexagons around 50m x 50m in size (real life), and my project is in long/lat. I hope that someone can help. Thanks

Sal


texts: Count number of pixels in different ranges for different polygons

texts:
I have two layers: one which contains few  seperated polygons, and one which has raster which I have uploaded (type feature collection). In addition, I have draw one for test (geometry).

My goal:

I want each polygon to get information based on its' location regard: how many pixels are in its' region that their value is lower than 0.2, between 0.2-0.4 ect.

What I have done:
I have used code that I had and tried to do this iteration for each polygon.

The result:
when running this with the GEE manual drawn polygon. calculation was done, but I got in the end imagecollection which contains this data in each feature properties. Also when I run this on my geometry (sometimes get error: geometry.geometry is not a function).

Code:

  var c0204 = image.gt(0.2).lt(0.4).selfMask().reduceRegion({
    reducer: ee.Reducer.count(),
    geometry: geometry,
    crs: 'EPSG:4326',
    //maxPixels:1310361348,
    scale: 20,
  }).values().get(0);

This is a link to my code, let me know if the imports are unavilable.

https://code.earthengine.google.com/24f5e8f9ef310b96d632ef89454c00eb

My end goal: To have the count pixels data in table so I can know for each polygon how many pixels were found in each count category. or at least to have it as featurecollection and not as imagecollection


texts: How to introduce table field parameter into field calculator formula in QGIS?

texts:
I am trying to build a very simple model in QGIS.
Simple intersection process (Density Zone-as intersection layer; watershed-as the other layer) followed by adding an "Area(ha)" field and finally the calculation of the population by multiplication of "Density Zone" field with "Area(ha)" field.
However, I want to add a table field parameter "Density Field" for example (density zone) which can be selected by the user from the "intersection layer".
I couldn't call the field parameter value into the field calculator algorithm.
I tried typing the functions:
Density Field * "Area(ha)"

"Density Field" * "Area(ha)"

'Density Field' * "Area(ha)"

I searched for "advanced python don't calculator". I don't know python code syntax so I couldn't manage to do so.
I also tried following texts in the advanced Python field calculator:
value= getvalue(<Density Zone>) * "Area(ha)"

value= getvalue(<Density Zone>) * (Area(ha))

All I wanted to do is to link the table field parameter into the functions' formula.

The problem is that a table field parameter cannot be introduced into a formula in order to create a new field.
Let's say there is a vector layer as an input. And another vector layer which will be used as an intersection layer. The first method is to intersect these two layers. The output of this intersection process will be used as a result layer but with further modifications. Secondly, in this processed layer, an area field will be added by the field calculator operator. Up to here, there is no problem.
Next, a second field will be created which is the multiplication of the area field previously created with the field which will be selected from the intersection layer (this field should be introduced as a "table field parameter" and must be selected by the user-user must know which field that is for this particular problem it is "density" field). The resulting output vector layer will have two new fields created which are "area" and "population" which is the multiplication of "area" and "density" fields.
I am struggling is how to type the multiplication formula by introducing the "table field parameter" (density) into the formula.


texts: What web mapping libraries are suitable for indoor maps?

texts:
I am looking for a custom mapping library but not to display an actual geographical area but to display a custom floor-plan. I have been looking for a while but most of the solutions that I found are for geographical areas.

I would like to customise a few zoom levels (probably 3 levels) and ideally I'm looking for a JavaScript, HTML5 (canvas, SVG) solution.

Does anyone have any suggestions?


texts: Calculate field in model builder using file name

texts:
I have browsed answers to similar questions but I still can't figure it out.  I have created a new field called "namemonyr" in all of 646 shapefiles.  I now want to populate that field using (preferably part of) the input file name.  My files are currently named e.g. "cl_Arthur_Apr_12.shp". I want the namemonyr field to be filled in with the string "Arthur_Apr_12". The name, month and year elements are different in every shapefile name.  I need to know when I output the tables as CSV files which subject the data came from and for what period - I need this info in the table, not just in the file name.  I have the iterator working fine.

I have tried using the field calculator in model building using the expression

!%Name%[3:]! to remove the cl_ element of the file name.

The error message reads:

ERROR 000539: Invalid field
cl_Arthur_Apr_12
failed to execute (Calculate Field)

I have tried using the Parse Path tool to get the file name as the Value that then goes into the Expression part of the field calculator.  But I don't really understand this tool.

I have tried variations on these.  Nothing is working!

I am very new to ArcGIS and Python and have never used Python coding.


texts: QGIS refactor field in processing modeler

texts:
I am trying to create a model to standardize data received from various clients into one single output format.
To do so, I use the refactor field and field mapper, following the answer given here: How to use fields mapper in qgis modeler？
One vector input is the client layer (unformatted)
Another vector input is our company standard layer (formatted)
Then the model becomes as follow.

The problem is that when I want to select in the list of "Source expression" my field from the client, I cannot find it. I only have the field from the formated layer available (see below).
When I do it using the algorithm directly from the processing toolbox, it is working well.
Am I missing something ?

I am using QGIS 3.14.16
EDIT:
Thank you MrXsquared. However this method doesn't allow the user to choose for each field of "client layer" which one is matching with the field of "company layer".
Sorry, I think my question was not developped enough. My goal is to reproduce the effect of Small ETL plugin (https://www.youtube.com/watch?v=yISFnyYCbqs&feature=youtu.be ) by using only the processing modeler (I want to share a project template to colleagues with the model saved in the template to avoid installing plugin).
I found a workaround: for each field, request the user to select it in the list and then rename it as the model, then refactor.
See the picture below:


texts: Finding Network Analyst tools in ArcGIS Pro?

texts:
I am trying to run some network analysis within 3D buildings to study handicap accessibility, best routes, emergency evacuation routes, etc. In ArcMap 10.3 there is a simple and easy to use toolbar that allows you to run the network analysis "on the fly" (directly to the map) and it is very quick to reset the start and end paths, restrictions, barriers, etc. With ArcGIS Pro so far I am able to run the network analysis, but it is very slow and requires individual inputs and outputs to be ran one at a time. With a large amount of starting and destination nodes this can be very tedious.

I was wondering if there is any way to get ArcGIS Pro to run the network analyst similar to the ArcMap 10.3 toolbar, or if there were any other options that can be used within ArcGIS Pro.


texts: How to set the map extent so as to show conccurently all the features included in different map layers

texts:
I want to load may diffent layers on my map using the ArcGis JS API.
I want to set my map's extent so as to show all the features of all the layers concurrently.
In other words I am looking for the extent that contains all my layers graphics and zoom to it.

How can I tell the map to find and zoom to a specific extent that will include all that features?
Is there any way other than iterating over all the features, like it is happening in the suggested solution of the this example ?


texts: Using MapInfo to create legend which shows points within named polygon?

texts:
I've created a map which shows some pupil post codes as points on the map and it also contains a layer which shows Wards within the city.

I would like to create a legend or table which shows how many pupils reside in each ward and the councillor responsible for the ward.

Is there a way to do this in MapInfo? I have already created a table in excel which shows this but i need the information to be on the map.


texts: I'm having an issue where when I save a QGIS project with a composer map made then reopen it all the map frames are blank with no graphic from the data. Is there something that would cause this? I've tried refreshing the render, changing the CRS, and the only thing that seems to work is adding another map frame. This only works until I reopen the file again.


texts: Calculating semi-major and semi-minor axis length of an ellipse in QGIS

texts:
I have a layer that contains ellipse-shaped geometries. The layer has projected CRS and there are also rotated ellipses in the layer. I couldn't find a tool or plugin related to semi-major and semi-minor of an ellipse.
Is there a way to calculate semi-major and semi-minor axis length?

I have added the blue points and green lines just for representation.


texts: Assigning projection to raster with PyQGIS?

texts:
Is there a processing algorithm in the pyqgis API for assigning a projection to a raster (GeoTIFF)?

I am clipping a GeoTIFF with a selected feature of a vector layer (shapefile).  Both are projected in GDA94/MGA Zone 55.  The QGIS project is set for the same projection. The result is defaulting to WGS84 for some reason, but not projecting to WGS84.  It is holding the GDA94/MGA Zone 55 coordinates.  I need to make sure the output is also in GDA94/MGA Zone 55 before progressing to the next step of the program.

I am running this in QGIS 2.14.  I ran it two years back in QGIS 2.8 and I got the same projection in the output as I had in the input.  I have tried manually running "Raster... Projections... Assign Projections" on the input raster before the clip, to make certain the input is georeferenced in a "QGIS friendly" manner.  No luck.  Here is the clipping algorithm I am running.

processing.runalg("saga:clipgridwithpolygon", veg_raster_for_clip, selected, outraster)


texts: Why can't I add the the qgis.org repository public key to Ubuntu 12.04

texts:
I'm trying to install QGIS in Ubuntu 12.04 but I can't import the public key:

mferreira@cibiows1:~$ gpg --recv-key 997D3880
gpg: requesting key 997D3880 from hkp server keys.gnupg.net
gpg: keyserver timed out
gpg: keyserver receive failed: keyserver error

Even if I try sudo it it gives an error message:

mferreira@cibiows1:~$ sudo gpg --recv-key 47765B75
[sudo] password for mferreira:
gpg: WARNING: unsafe ownership on configuration file `/home/mferreira/.gnupg/gpg.conf'
gpg: external program calls are disabled due to unsafe options file permissions
gpg: keyserver communications error: general error
gpg: keyserver receive failed: general error

Does anyone I have any idea how to work around this?


texts: Adding values from a table to an image band - Google Earth Engine

texts:
I'm trying to add values from a table to a cluster image. The cluster image has 50 cluster classes, and my table has 50 values for each cluster class. How do I add the values from that table to the image so that each pixel within that cluster class represents that value from the table? I've included an image for clarification.

I understand this might get flagged as too broad, but I've reached the limit of my ability, and I'm in need of some direction. Is something like this even possible in GEE?

Link to my code: https://code.earthengine.google.com/0d4094a289c9f7f139809b4e9e1f1d05


texts:
I am trying to calcualate the statistics of a TIFF. The script should ignore the NoData value, but doesn't. What causes this?

import gdal
from gdalconst import *

# register all of the drivers
gdal.AllRegister()

# open the image
ds = gdal.Open('test_slope.tif', GA_ReadOnly)

# get raster band 1
band = ds.GetRasterBand(1)

# set NoData value
band.SetNoDataValue(-3.4028230607371e+38)

# Calculate and print Statistics
stats = band.GetStatistics(0,1)
print stats

The results look like this

>>>
[-3.4028230607371e+38, 126.5991897583, -4.6432918375055e+37, 1.1680875238428e+38]
>>>


texts: Creating squared halo instead of circular in TileMill?

texts: Is it possible to style a halo around text that has a square shape instead of a circle? If not, what are the alternatives to make a "halo" around number angular, in a shape of a square? Software I am using is TILEMILL. On this picture below i posted, you can see blue circle around number which is halo. I don't want for a that blue halo to be a circle but a square instead.


texts: Creating multiple layer from single store in GeoServer?

texts:
I am new to GeoServer and using version 2.9.2.

I have a store based on a shapefile with multiple multipolygon features. And I want to create multiple layer based on CQL attribute query from that single store.

I try to add new layer and choose that store, but all I found is "Publish again" and "Create new feature type" link. I try the "Create..." link but it gives me error and I can't find any documentation about it.

Is it possible?


texts: Calculating Kernel Densities Using A Loop in Geospatial Modeling Environment

texts:
I'm currently trying to loop through my shapefile of telemetry locations to calculate kernel densities using the 'PLUGIN' bandwidth. I tried the following code in R to create a script that can be loaded into Geospatial Modeling Environment:

uniqueids <- unique(Locations_2001$Animal_ID)

for (i in 1:length(uniqueids))
  {cat(”kde(in=\”Locations_2001.shp\”,out=\”kde_”,uniqueids[i],”.img\”, bandwidth=\”PLUGIN\”,cellsize=30,where=\”Animal_ID=’”,uniqueids[i],”’\”);\n”, sep='')}

I literally built the code from the code listed in the Geospatial Modeling Environment PDF:

uniqueids <- unique(ANIMALID)
for (i in 1:length(uniqueids)){
cat(”kde(in=\”telemetry.shp\”, out=\”kde_”, uniqueids[i], ”.img\”, bandwidth=\”SCV\”, cellsize=20,
where=\”ANIMALID=’”, uniqueids[i], ”’\”);\n”, sep=””)}

My end goal is to get a script like the one in the PDF that can be saved and loaded into GME:

kde(in=”telemetry.shp”, out=”kde M01.img”, bandwidth=”SCV”, cellsize=20,
where=”ANIMALID=’M01’”);
kde(in=”telemetry.shp”, out=”kde F01.img”, bandwidth=”SCV”, cellsize=20,
where=”ANIMALID=’F01’”);
kde(in=”telemetry.shp”, out=”kde M02.img”, bandwidth=”SCV”, cellsize=20,
where=”ANIMALID=’M02’”);
kde(in=”telemetry.shp”, out=”kde F02.img”, bandwidth=”SCV”, cellsize=20,
where=”ANIMALID=’F02’”);
kde(in=”telemetry.shp”, out=”kde J01.img”, bandwidth=”SCV”, cellsize=20,
where=”ANIMALID=’J01’”);

However, I get the following error when I run my code above:

Error: unexpected input in:
"for (i in 1:length(uniqueids))
  {cat(”"
where=\”Animal_ID=’”,uniqueids[i],”’\”);\n”, sep='')}
Error: unexpected input in "where=\"

What am I doing wrong in the above code?


texts: Modifying the TIFF's attribute table value with R

texts:
I use raster::raster("test.tiff", RAT = T) to read the file with the format of TIFF.
The TIFF is derived from ArcGIS pro generating using Feature to Raster -> raster -> tiff.
ArcGIS TIFF's attribution table as below:

As you can see, the fields of OID, Value, Count are not allowed to edit in ArcGIS Pro. But I really need to make Value <- Name. Thus, though ArcGIS Pro doesn't allow I to modify Value field, I could use other platform like Python or R to edit the value , and finally I resave the modified data to the TIFF But I don't know how to resave the TIFF. And I learned another way as below:
db = read.dbf('test.tif.vat.dbf')
db$Value = factor(db$Name)
write.dbf(db,'test.tif.vat.dbf')

Even though I can save without any warning, but finally , when I open the TIFF in ArcGIS Pro , there's any change in attribution table .
The reason that I want to modify the value of attribute

Can you tell me what should I do ?


texts: Saving duplicate geometry layers into a geopackage?

texts: I have a bunch of variables to map with same geometry (Regional Boundaries), and wanted to create separate layers with differing attributes within a single geopackage file for ease of file management. However, when I save the duplicate layer into a geopackage I get an error "Export to vector file failed. Error: Cannot overwrite a OGR layer in place". This seems to be expected behavior according to this bug report. Is there a way around this? Creating separate gpkg files for each layer seems to be unnecessary/superfluous. Or am I misunderstanding the entire concept of geopackage as a an analogue to Esri's geodatabase?


texts: Writing attributes and associated feature class name from file geodatabase to CSV?

texts:
I need to write a csv file that lists every record from an ID field across all feature classes in a gdb. The code I have written below does that fine, but I also need to write in the second column of the csv which feature class the record was found in.

I have been trying to create the cursor so that it also contains the fc name, then when I use writerow(row) both fields will be printed, but I think this is not a viable approach. Can anyone point me in the right direction?
I'm using arcgis 10.2.2

arcpy.env.workspace = r'C:\Data.gdb'
output = r'C:\test.csv'

fcs = arcpy.ListFeatureClasses()
fieldnames = ['PropID', 'Feature Class']

with open(output, 'wb') as f:
    csvwriter = csv.writer(f)
    csvwriter.writerow(fieldnames)
    for fc in fcs:
         with arcpy.da.SearchCursor(fc,'PropID','PropID IS NOT NULL') as cursor:
            for row in cursor:
                csvwriter.writerow(row)
    del row, cursor


texts: Can you explain this WKT multipolygon format?

texts:
I am confused by the output I'm getting from ogr2ogr when converting a shapefile to WKT. What I'm expecting to see is something in the form below (from Wikipedia):

MULTIPOLYGON (((40 40, 20 45, 45 30, 40 40)),
((20 35, 45 20, 30 5, 10 10, 10 30, 20 35),
(30 20, 20 25, 20 15, 30 20)))

What I'm getting out of ogr2ogr -f CSV out.wkt source.shp -lco GEOMETRY=AS_WKT, in very abbreviated form, is this:

"MULTIPOLYGON (((-120.6 50.8, ...)", 33.9, 137.1, 2, 1,1,1

I do not understand how to interpret the final six numbers, which are not enclosed in brackets and fall outside the quotation marks. Can someone explain this or point me towards the appropriate documentation?


texts: Determining Min - Max within zip code (statistics by area)?

texts:
I'm new to ArcGIS Desktop. I have created a layer representing zip codes and overlayed it with a layer containing some data points. I want to find the min and max value of the data points in each zip code. Currently, I am manually checking the values if all the data points in each zip code but this is becoming very time-consuming and cumbersome.

Is there a faster way to do this?


texts: Publish existing table from postgis to geoserver

texts:
I would like to know if there is any way for publishing existing table from postgis to geoserver using Geoserver REST API or geoserver-manager java implementation.

The thing is that i'm actually working with java, and it's sending this error message:

ERROR [geoserver.rest] - :Resource named 'aerodromos' already exists in store: 'appdemo'

However running this command from terminal works fine:

curl -v -u admin:geoserver -XPOST -H "Content-type: text/xml" -d "<featureType><name>aerodromos</name></featureType>" localhost:8080/geoserver/rest/workspaces/appdemo/datastores/appdemo/featuretypes

But then, I need to setup the xml configuration, one thing that with geoserver-manager implementation works using this code:

public boolean createLayerOnPostgisDatastore() {
    GeoserverConfig conexionGeoserver = new GeoserverConfig();

    String workspace = parametros.get("workspaceGeoserver").toString();
    String storeName = parametros.get("datastoreGeoserver").toString();
    String proyeccion = "EPSG:4326";
    GSLayerEncoder layerEncoder = new GSLayerEncoder();
    layerEncoder.setEnabled(false);
    layerEncoder.setQueryable(true);

    if (parametros.get("tipoGeomCapa").toString().equals("POLYGON")) {
        layerEncoder.setDefaultStyle("polygon");
    } else if (parametros.get("tipoGeomCapa").toString().equals("LINESTRING")) {
        layerEncoder.setDefaultStyle("line");
    } else if (parametros.get("tipoGeomCapa").toString().equals("POINT")) {
        layerEncoder.setDefaultStyle("point");
    } else if (parametros.get("tipoGeomCapa").toString().equals("MULTIPOLYGON")) {
        layerEncoder.setDefaultStyle("polygon");
    } else if (parametros.get("tipoGeomCapa").toString().equals("MULTILINESTRING")) {
        layerEncoder.setDefaultStyle("line");
    } else if (parametros.get("tipoGeomCapa").toString().equals("MULTIPOINT")) {
        layerEncoder.setDefaultStyle("point");
    }

    GSFeatureTypeEncoder fte = new GSFeatureTypeEncoder();
    fte.setName(hlsNombreCapa);
    fte.setTitle(hlsNombreCapa);
    // fte.addKeyword("TODO");
    fte.setNativeCRS(proyeccion);
    fte.setDescription(parametros.get("descripcionFinal").toString());
    fte.setEnabled(false);
    fte.setLatLonBoundingBox(getShapeFile().getHeader().minX(), getShapeFile().getHeader().minY(), getShapeFile().getHeader().maxX(), getShapeFile().getHeader().maxY(), proyeccion);
    boolean ok = conexionGeoserver.getPublisher().publishDBLayer(workspace, storeName, fte, layerEncoder);
    return ok;
}

Thanks in advance.


texts: geting geograpic locations from an address list

texts: I have a table with 150 addresses. I would like to automaticly generate a point shapefile/PostGIS layer from all of them. Is that possible in any way or do I have to search for every single address to find its geographic location?


texts: When building a VRT in QGIS, which input raster will be assigned which band number?

texts: I don't understand how QGIS assigns the band numbers to the input rasters for a virtual raster catalogue. Unfortunately the bands in the compiled virtual raster catalogue will only have the names "Band 1", " Band 2", "Band 3" and so forth, but the name of the original input raster will not be preserved. This causes a problem when trying to create a false colour image (render type: "Multiband colour" option in QGIS) with Landsat 8 satellite imagery.


texts: Is r.cost equivalent to diffusion interpolation with barriers?

texts:
I want to perform interpolation using cost distances rather than euclidean distances similar to Greenberg et al. (2011). In that paper they use r.cost in GRASS to compute cost distances for each prediction location. Unfortunately I am more familiar with ArcGIS than GRASS. Can I use the diffusion/kernel interpolation with cost surfaces in ArcGIS rather than GRASS/r.cost? Would I get equivalent results?

Greenberg, J.A., C. Rueda, E.L. Hestir, M.J. Santos and S.L. Ustin. 2011. Least cost distance analysis for spatial interpolation. Computers & Geosciences 37:272-276. http://dx.doi.org/10.1016%2Fj.cageo.2010.05.012


texts: How to position LeafLet CircleMarkers on top of Markers

texts:
Question 1: With mapbox.js 2.3.0 (not sure what LeafLet version), are CircleMarkers supposed to go into the overlay pane by default, and Markers into the marker pane?

I basically have 3 sets of objects (labels, contours, points), which I'd like to display in the following z-index order (i.e. bottom to top): contours < labels < points.

My points are CircleMarkers and labels are Markers. If both were in the same pane, I could manipulate the z-index of each at instantiation. But since contours and points are in the overlay pane, and the labels in the marker pane, and LeafLet makes the z-index of the marker pane much higher than that of the overlay pane, the order I get is: contours < points < labels.

I saw this thread, whose best answer suggests to manipulate the z-index of the panes via CSS. The issue with that is that I end up with: labels < contours < points. It is better but not ideal, so...

Question 2: is there a more robust way of manipulating the z-index of CircleMarkers and Markers? It is irrelevant to me which panes they land on.

Below is how I instantiate points and labels, if it helps:

 // Labels:
 L.geoJson(file , {
       onEachFeature: function(feature, polygon) {
          var myIcon = L.divIcon({
              className: 'iconCSS',
              html: "
" + feature.properties.LAD13NM + "
"
          });
          var label = L.marker(polygon.getBounds().getCenter(), { icon: myIcon });
          map.addLayer(label);
        }
    });
// Points:
L.geoJson(file, {
  onEachFeature: function(feature, point) {
    var station = L.circleMarker(point.getLatLng());
    map.addLayer(station);
  }
});


texts: Creating minimalistic map in *.svg file format?

texts:
I am interested to obtain a .svg map file format in two colors, without labels or marks, like we have here or here for example.

what I've tried:

to use Google Styled Maps Wizard but unfortunately I'm unable to save the result into a decent file. I've said decent because, there is a possibility to use the Static Map button that will save a tiny image (not the desired result). Then,
to use Maperitive applying wireframe.mrules and Schwartzplan.mrules but this is not what I'm looking for; in this case the buildings are highlighted (not the streets on the land color background).

I find modifying the rules in Maperitive without seeing the layers in real time, pretty uncomfy; Having said that, is there any wysiwyg plugin for this app or a standalone editor (pre-Illustrator importing phase) that can support me to achieve the described result and also can save layered svg files?

or Is there anyone that already have and want to share a Maperitive .mrules file that can reproduce the above mentioned result?

@Mapperz I've tried to follow the below steps but I think that there is still a misunderstanding (I have to recognise that could be my lack of knowledge).

Download (Shapefiles) - Done! I've downloaded 828 KB coastline shapefile, 25 MB osm2pgsql shapefiles and 36 MB imposm shapefiles.

Load into QGIS - In the above described archives I've found/load vancouver.osm-roads and vancouver.osm-buildings. Is that enough?

and Change the map background to white The map background is already white if you reffer to Project properties/General/Background color (?) and Load the polygon and line work into QGIS already done that at the above step Load into QGIS, are you reffering to something else?

Set the layer black In the Print composer Should I create a new Print composer I guess but then, which layer should I make it black please?


texts: How to divide a buffer using geometry generator in QGIS 3.X

texts:
Some background and working PostGIS query and Python script available here.

How to split circles in 12 sections in PostGIS?

I don't want to create a new polygon layer out of points. I really like QGIS Geometry generator. My idea is to use QGIS Geometry generator and show this divide buffer as symbol.

Any suggestion?


texts: Appropriate resolution for regression analysis

texts:
I'm using ArcMap 10.3. I'm a student working on a research problem where I'm taking Landsat imagery (at 30-meters resolution) and regressing it against other datasets captured at different scales (up to about 1 km^2). I am then resampling everything down to 30 square meters prior to running the regression.

My ultimate goal is to understand vegetation dynamics, which is why I wanted to bring the whole analysis down to Landsat's 30 square meters. Or is this thinking incorrect, ie, should I run the analysis at the scale of the coarsest dataset, being 1 km^2?

From what I've read- in particular, Hengl 2006 "Finding the right pixel size", there is no "ideal" grid resolution and it is dependent on the question. Another paper I read, Tarnavsky 2007 "Multiscale geostatistical analysis of AVHRR, SPOT-VGT, and MODIS global NDVI products" recommended choosing a pixel size to minimize the spatial variability between sensors. Yet other papers suggest statistical corrections in order to use multiscale data.

I can do the regression. My question is whether there are any standards for how to approach an analysis where all the input data is at a different scale, in my case, ranging from 30 square meters up to 1 km^2.


texts: How to find out to which polygon each points belongs to?

texts:
I have searched for an answer to this, but couldn't find anything useful. I have got a shapefile consisting of a number of polygons, and another shapefile with a number of points.

I want to add a column in the points shapefile. For each point, I want the column to contain the name of the polygon that surrounds it.

In simple terms, I just want to find out to which polygon each points belongs to. Any tips?


texts:
My code is:
import arcpy
import os
from arcpy import env
import math
#arcpy.env.workspace = "CURRENT"

env.overwriteOutput = True

try:
    Nearest_line = arcpy.GetParameterAsText(0)
    Origin_locations = arcpy.GetParameterAsText(1)
    Destination_locations = arcpy.GetParameterAsText(2)

    origin = 'orig'
    destination = 'dest'
    origin_points = []
    destination_points = []

    #Create feature class
    arcpy.CreateFeatureclass_management(os.path.split(Nearest_line)[0], os.path.split(Nearest_line)[1], "polyline", spatial_reference=2039)
    arcpy.AddField_management(Nearest_line, 'Orig_Id', 'SHORT')
    arcpy.AddField_management(Nearest_line, 'Dest_Id', 'SHORT')
    arcpy.AddField_management(Nearest_line, 'Dist', 'FLOAT')

    #Loop each layer to create list points
    with arcpy.da.SearchCursor(origin, ['OBJECTID', 'SHAPE@XY']) as origins :
        for point in origins :
            origin_points.append(point)
    with arcpy.da.SearchCursor(destination, ['OBJECTID', 'SHAPE@XY']) as destinations :
        for point in destinations :
            destination_points.append(point)

    #Calculating the distance between two points from the layers above
    distan = []
    for i in origin_points :
        Dista = []
        x1 = i[1][0]
        y1 = i[1][1]
        Orig_Id = i[0]
        for j in destination_points :
            x2 = j[1][0]
            y2 = j[1][1]
            Dest_Id = j[0]
            Distance = math.sqrt((x2 - x1) * (x2 - x1) + (y2 - y1) * (y2 - y1))
            Dista.append([x1, y1, x2, y2, Orig_Id, Dest_Id, Distance])
        Distm = min(Dista, key=lambda l : l[6])
        distan.append(Distm)

        #Create list of fields to create the lines
        lines = []
        for L in distan :
            array = arcpy.Array([arcpy.Point(L[0], L[1]), arcpy.Point(L[2], L[3])])
            sr = arcpy.SpatialReference(2039)
            Orig_Id = L[4]
            Dest_Id = L[5]
            Dist = L[6]
            lines.append([Orig_Id, Dest_Id, Dist, arcpy.Polyline(array, sr)])

    #I scanned each line and inserted it to the table in the Feature Class output (Polyline)
    with arcpy.da.InsertCursor(Nearest_line, ['Orig_Id', 'Dest_Id', 'Dist', 'SHAPE@']) as cursor :
        for line in lines :
            cursor.insertRow(line)

except Exception as e:
    print("Error: " + e.args[0])

I get error:

Error: Failed to execute. Parameters are not valid. ERROR 000735:
Feature Class Location: Value is required ERROR 000735: Feature Class
Name: Value is required

Failed to execute (CreateFeatureclass).


texts: OpenLayers: Finding features in pixels radius

texts:
I'm using openLayers 3.19.1.

When clicking on the map  (map.on("click") event ) I retrieve the pixel I clicked on (event.pixel).

I need to get all features in an X radius on the pixel.

There is a function called "forEachFeatureIntersectingExtent"  I can execute on the ol.source.Vector, but it needs an "extent" parameter.

How can i convert my pixel to an extent?

If anyone has another idea for getting all features in a specific radius of a click on the map please share.


texts:
I want to process a raster with values 0-16. I want to set all values between 1 and 15 as 1, a value of 16 as 0 and cells with a value of 0 as 0.

How do I maintain a cell value of 0?

I have tried the following:

("raster_name@1" !=0 AND "raster_name@1" <= 15*1)+("raster_name@1" = 16*0)

Results in: Values of 16 as 0, values 1-15 as 1, values 0 also = 1

I also tried:

("raster_name@1" !=0 AND "raster_name@1" <= 15*1)+("raster_name@1" 16*0)+("raster_name@1" = 0*0)

Results in: Values of 16 as 0, values 1-15 as 1, values 0 = 2

Please, can you explain how to retain cells with a 0 value as 0?


texts: Creating tracks for moving points (while discarding records without movement)

texts:
I have a csv with columns:

id
latitude
longitude
timestamp

Many of the points (or "id"s) do not change "latitude"/"longitude" over time, others do.
How do I filter out those "id"s and "lat"/"long"s for the ones that did not move?
For the "id"s that do move, there's a path. How do I export or select only the paths or the "lat"/"long"s of those ids that moved?
Is this possible with MMQGIS or another time-series Plugin?


texts: Different length measures in QGIS vs. ArcMap

texts:
I noticed that QGIS and ArcMap return slightly different length measures when calculating geometries (length) on one and the same shapefile (see image).

In QGIS I used the $length command, in ArcMap the "Calculate Geometry" function from the attribute table.

Is there an explanation for the difference observed?

The data are in the Lambert Conical Conformal projection, with parameters used by the Mexican authorities (INEGI). Units are meters and I use the data’s coordinate system (the LCC)

And generally speaking: How exact can I expect the length-calculus to be? Do I have to consider a confidence interval?


texts: ArcGIS align grid of 10km squares with grid of 1km squares

texts:
I have a pre existing layer of 1km x 1km grid squares, which is aligned with the 10km squares of the British National Grid:

I want to create another 10km x 10km square grid on top to highlight the 10km grid lines. I have been trying to create a grid using the create fishnet tool, but I can't figure out how to align the resulting grid with my previous one. Does anyone know if there is a way to align a fishnet with a preexisting feature, or if there is any other way to achieve this aim? Basically I want the black gridlines on the green transparent grid layer below to line up with the blue gridlines on the basemap layer underneath. The projection of all the layers is British National Grid.


texts: Multiple GDAL operations on Raster Dataset without opening and closing files multiples times

texts:
I'm going to start writing a script, likely in python, where I need to open a raster file, reproject it, multiply it with a slope-shaded elevation model then save the full res outcome as well as multiple downsampled versions in specific resolutions (1m, 2m, 4m, 8m, 16m, 32m).

My question is, how do I go about opening a file and doing multiple operations on it, then save the outcome without the overhead of opening, exporting,  then opening exported files for the next operation, multiple files along the way using gdal libraries?
I'm basically trying to stay away from using each gdal terminal program individually as linking those programs together causes lots of disk reads and writes..

I'm hoping for a general example, even if it doesn't do the same functions that I'm looking for.


texts: Normalization of stack of features in Google Earth Engine (GEE)

texts:
How to normalize a stack of features in the google earth engine? Features are the image type you can see the Javascript code here. I'm going to normalize features between 0,1 values.
Can you provide suggestions and snippets?
var ***features*** = (feature 1

    .addBands(feature 2)
    .addBands(feature 3)
    .addBands(feature 4)
    .addBands(feature 5)
    .addBands(feature 6)
    .addBands(feature 7)
    .addBands(feature 8)
    .addBands(feature 9)
    .addBands(feature 10)
    .addBands(feature 11)
    .addBands(feature 12)
    .addBands(feature 13)
    )


texts: Connecting centroids of lines sharing same intersection in QGIS?

texts:
I have this street network of the center of Turin (IT), composed by lines representing streets, and black points representing intersections.

I have computed the centroids of each line and displayed these as red points. Now I would like to connect the centroids (red points) of each line that share the same intersection (black points) as in the "B" scheme of this picture:

The attribute of the intersections, that are the black points, are simple Ids:

Whereas in the attributes of the lines representing the streets I have the Ids, the start node (that is the Id of the black point where the line begin) and the end node (the Id of the black point where the line end.)

Do someone have an idea about how this can be done?

I've thought about pharsing the lines start points over the end points and vice-versa, but I really don't know how to do this.


texts:
I am using QGIS Essen on Windows 10. I have a data set (.txt file) that contains latitude, longitude, and black carbon (BC) values. I can style the layer using the values of BC - graduated style, color ramp. As such:

However, I would like to have the colors in log scale so I can still distinguish between the different low values (yellow dots). As you can see on the right side, all is just yellow. Is there a way to do this?

Update: I tried again the field calculator approach and got a nice plot. However, I would like the color ramp "logged" and not the values (please let me know if this makes sense). I'm trying to achieve something like this that I have done in Igor (but the graphics is better in QGIS):

As you can see, the color scale is also logged. Is this possible in QGIS?


texts: Adding shadedreliefarchive imagery to QGIS

texts:
I'm trying to add some relief imagery that is found at shadedreliefarchive.com

Specifically, I'm trying to add the SEAsia one found here.

I've tried making a custom CRS to match the details found on the page, but the units are way off and I'm not sure which units should be used anyway.  My definition is below.  Can anybody help?  Either with the correct already defined CRS or the user string for a custom CRS?

+proj=lcc +lat_1=10 +lat_2=21 +lat_0=17 +lon_0=103 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs


texts: Proxy for Secure ArcGIS Server Services. Secure?

texts:
To access secure services (token based or credentials), Esri recommends using proxy files  (.net github example). When routing requests through the proxy, you are able to request secure services on behalf of the client without exposing your credentials.  You can define a property called allowedReferers and assign a list of referring URL's that the proxy will work for.  Basically, the proxy won't make any requests for referring URL's that are not defined.  If set to '*', any referring request will be processed.

The problem is; the requesting header be spoofed easily by a hacker by just setting a false HTTP Referer property.  In this situation, they can access secure services by routing all of their requests through the proxy and setting the referer header to a valid address.

I am looking for recommendations on the best way to work around this issue.  Any recommendations?


texts: Generating a buffer of a vector line based on raster values where the vector intersects a raster

texts:
I am trying to figure out how to buffer a vector shapefile with a variable distance (in QGIS). However, the variable distance needs to be dependent on the value of a raster (DEM) where that vector overlies the raster.

Put differently, the yellow buffer extends out to a distance that depends on the height of the DEM above some known height at the point where the purple line intersects the raster.


texts:
I have a large collection of drone images that aren't georeferenced. I'd like to georectify them by pulling projection and coordinate data from the GeoTiff files, but I'm having difficulty.

My initial thought is to use GDAL in Python to pull lat/long data and use that to georectify the images, but I don't have much experience with GDAL or georectification through scripting.

Any advice to help me get started?


texts: How to edit an attribute while editing a feature?

texts:
I have several polygons to edit / modify and I want to modify some fields in the attribute table at the same time.

I have the table open, but I do not see where the attributes of the polygon I'm editing are right away.

To do this I can:

Hit the "Select Feature" tool
Select the polygon
Go to the table and select "Move selection to top"
Re-select the "Node" tool.
Modify the polygon; one or several nodes, or maybe nothing.
If needed, go to the table again and change the necessary fields.
Pass to another polygon, and do it again.

Or I can use the "Identify feature" and "Auto open form".

With these ways, I walk back and forth to the same tools. It becomes annoying and diverts attention from the information I need to control ...!

The problem:

When I´m using the "Node tool" I can´t see the attribute that I´m modifyng or open a form.
To open a form I need use the "Identify feature" or the "Select feature" tools and got back to "Node tool".

I need to focus on the "Node tool".

Is there any way I can access a form that opens automatically while I use the "Node tool"?

Or is it possible the attribute table move to the top the line I'm modifying, when using the "Node tool"?

(QGIS 2.18.14)


texts: Organizing kml created from shapefile in QGIS?

texts:
I have a shapefile that I would like to export into a kml to share it with my colleagues.

This shapefile contains a few hundred polygons. I have QGIS 2.2.0 on Windows 7.

I would like to create a kml that somehow keeps the same style as my shapefile. Indeed, in QGIS, I have a simple category style applied to a certain field ("Soil Type"). When I export the shapefile as a kml (after renaming the "Soil Type" field as "Name" to have it as a title of my kml polygons), the kml created doesn't have any sub-folders per Soil Type. As a result, if I want to set a color per Soil Type, I need to manually create folders in Google Earth, move every single polygons into the different new folders, then apply a style to the folders. It works, but it takes a long time when you have a lot of polygons.
Another solution I found is to create a different kml of each "Soil Type" in QGis. But again, it's not very efficient.

Is there any solution more efficient to organize the kml created from a shapefile in QGIS?


texts: Labeling one to many relationship class using ArcGIS for Desktop?

texts:
I’ve created a feature class with a one to many relationships to a table.

One point many classes.

However I want to label the point with the classes attributes.

Esri wrote VB solution but it's not supported at version 10.1.


texts: GEE- Supervised classification image is incomplete

texts:
Using Google Earth Engine, I'm trying to do a supervised classification over a semi-urban area by placing points with different numeric classes on areas of vegetation, impermeable surfaces (roads and buildings), and exposed soil. I have about 50 points of each classification type and have merged them into a new feature collection. When I try to import a landsat image over my area of interest, the image only covers about half of the points. In the tutorial "Hands-on Intermediate Training", at about 12:10, the same thing happens but he is able to correct the issue by moving one point:

https://developers.google.com/earth-engine/tutorials

I'm using essentially the same code as in the tutorial. I've moved about 60 points around and nothing changes.

Has anyone else had this problem and had any success getting the image to cover all the points?


texts:
I'm currently developing a slippy map. It is perfectly working with tile coordinates that are used in openstreet maps i.e. the zoom/x/y format. I have used the formula from the OSM wiki to map longitude and latitude coordinates to their respective tiles.
However I also want to implement a method to use a WMTS Tile Server. The problem is that the coordinate to tile mapping seems to be different here.

Is there a formula to get the WMTS tilenumber from a specific coordinate in longitude/latitude format with a given zoom level?

I am using geoserver to serve WMTS tiles.

I constructed an example to show my problem.

var lat = 50.5524; //Coordinates of Wetzlar, Germany
var lon = 8.5072;
var zoom = 12;
var tileLong = long2tile (lon, zoom + 1); //For some reason I need to increase the zoom level by 1 to get equivalent results in WMTS and OSM
var tileLat = lat2tile(lat, zoom +1);

tileLong -> 4289, this is correct as the OSM and the WMTS tile have the same value
tileLat -> 2758, correct in OSM, incorrect for WMTS as WMTS needs inversing the y-axis

var invertedTileLat = (1 << zoom) - tileLat - 1; //I have also tried zoom+1

invertedTileLat -> 1337, still incorrect for WMTS. The expected result is 897


texts: How to transform isolines to isopolygons with postgis?

texts:
I have a table a postgis table of isolines that's defined like this:

CREATE TABLE myisolines
(
  gid serial NOT NULL,
  isotime timestamp without timezone,
  val numeric(10,4),
  geom geometry(LineString,4326)
);

Visually this linestring objects look like this:

I know the spatial extent of my data, so I can add a Bbox, so the LineStrings can be sort of closed.

I want to create a table of isopolygons myisopolygons from myisolines table, with polygons, that won't overlap but create a continuous surface and have a column val with the lowest val of isolines, from which the polygon was formed. I understand it can be formed from self-closed isoline (island), or isoline closed with bbox, in that case the val should be taken from that particular isoline. Visually it should look like this:

I thought that I can create topology somehow and then transform faces to polygons, but I don't understand how to do it properly. How can this be done?

Another option would be to recursively use a difference function between bbox and each created polygon, but I guess that's not the right way to do it, and definitely not fast at all.


texts: How can I validate a SQL statement in ArcObjects?

texts:
I'm trying to validate the results of a ISQLQueryDialog2, but can't find any out of the box solutions.  This particular dialog has a built in validate button, but I would like to programmatically validate the SQL string against the feature class in question.

I've found ISQLCheck3, but the documentation for this is slim, as it's reserved for internal use only.

Is there a way to validate a SQL where clause?

EDIT: I've created a try/catch statement to attempt to use a QueryFilter, however this seems kind of hacky.

IQueryFilter qFilter = new QueryFilter();
qFilter.WhereClause = qDialog.SQLQuery; //result of ISQLQueryDialog2
try {
    fc.FeatureCount(qFilter);
} catch (Exception ex) {
    MessageBox.Show("Houston, we have a problem");
}


texts: The Buffer tool in ArcGIS 10.5 allows for the creation of geodesic buffers, however, I am unable to find that option for the Multiple Ring Buffer. I have 10 areas around which I need to create 6 geodesic ring buffers (huge areas so quite a difference regarding area) and I'd like to use the multiple ring buffer tool to avoid having to clip each ring. Any ideas how to best achieve this?


texts: Is it possible to set a default so that on opening a QGIS file it doesn't attempt to render?