texts: Dissolve tool is splitting features rather than dissolving them

texts:
I had UK area divided into six different zones.
I wanted to dissolve my feature class by zone just so I have six features in it.
I used Dissolve tool which worked fine for five zones. One of the zones rather than being dissolved, got split into something that looks tiles.

I used the tool in ModelBuilder and also independently. And makes not difference if my output is stored in geodatabase or in the folder.
I also run repair geometry on my input.

What would be the reason for that and how to make Dissolve tool work correctly?

Please see the image below.


texts:
The code sharing repl.it has a python interpreter, and proudly hails that any library is available for import.

https://repl.it/site/blog/python-import

This does include the GDAL library. However, the typical GDAL installation experience, to the best of my knowledge, involves
1. installing drivers, libgdal-dev
2. setting environment variables
3. AND THEN installing GDAL.

Sequences like this have been discussed in this related question:
Python GDAL package missing header file when installing via pip

Given these extra shell requirements, is it impossible to install GDAL from a python interpreter alone, as seems to be required for the available repl.it interface? Or is there a way to make this happen from within python?

https://repl.it/@ValerieAnne563/RunningWithGDAL


texts: Looping through Layers, Add Layer to MXD, Export Image, Repeat, using ArcPy?

texts:
ArcGIS 10.2.2

I am trying to loop through a list of layers in a folder, add each layer to an mxd and export an image of the added layer. I have added arcpy.mapping.RemoveLayer to the script to remove the added layer before adding the next layer, but it's not working the way I thought it would and is adding all the layers to the mxd and exporting after each one is added without removing the layer. I'm not sure if it's not nested/indented correctly, or if I'm not saving the map document at the right point(s). So, to recap, the script is successfully adding each layer to the mxd and exporting an image after each layer is added, but is not removing the layer after exporting.

It should be Add Layer 1>Export Image 1>Remove Layer 1; Add Layer 2>Export Image 2>Remove Layer 2; Continue Through List.

Here is what I have:

mxdPath = mxd path on disk
mxd = arcpy.mapping.MapDocument(mxdPath)
df = arcpy.mapping.ListDataFrames(mxd)[0]

imageFolder = path to image folder on disk
lyrFolder = path to folder of physical layer files on disk

arcpy.env.workspace = lyrFolder
rasterLayers = arcpy.ListFiles("*.lyr")
for layer in rasterLayers:
    image = imageFolder+"/"+layer[:-8]+".jpg"
    layerPath = lyrFolder+"/"+layer
    addlyr = arcpy.mapping.Layer(layerPath)
    arcpy.mapping.AddLayer(df,addlyr,"TOP")
    mxd.save()
    arcpy.mapping.ExportToJPEG(mxd,image)
    arcpy.mapping.RemoveLayer(df,addlyr)
    mxd.save()

I've also tried, to no avail:

for layer in rasterLayers:
    ...
    arcpy.mapping.AddLayer(df,addlyr,"TOP")
    mxd.save()
    arcpy.mapping.ExportToJPEG(mxd,image)
    # Delete and recreate map document object after saving to reflect added layer in ListLayers
    del mxd
    mxd = arcpy.mapping.MapDocument(mxdPath)
    for lyr in arcpy.mapping.ListLayers(mxd,"D_*"): #All layers in lyrFolder start with "D_"
        arcpy.mapping.RemoveLayer(df,lyr)
    mxd.save()

I'm sure I've missed something simple, but I'm stumped at this point.


texts: OpenLayers 3 : Problem with the tracker?

texts:
I try to make a tracker with OpenLayers. The localization point appears well but not in the right place.

<title>OpenLayers example</title>


texts: Seeking alternative to Google Maps Elevation API

texts:
The Google Maps Elevation API lets me pass a URL like

https://maps.googleapis.com/maps/api/elevation/json?locations=39.7391536,-104.984703

along with a developer key, and get a result like

{
   "results" : [
      {
         "elevation" : 1608.637939453125,
         "location" : {
            "lat" : 39.73915360,
            "lng" : -104.98470340
         },
         "resolution" : 4.771975994110107
      }
   ],
   "status" : "OK"
}

My company want to use something like this, but don't have the time for their legal department to go through the license.

Is there something similar with an MIT, Apache, GPL2, or similar license?

Basically, we want to input lat/long and get back elevation.


texts:
I a beginner with the ArcGIS platform.

I would like to create a new style with ArcMap. I try with the menu : Customize -> Style Manager

Style Manager : http://hpics.li/9471c73

My new style appears but it seems to be empty. I would like to add new symbols on it. Exactly 12 images (png).

First of all, I don’t know where can I put them in all the different categories. None of them seems to accept images…Is this the good way to create a new style with new symbols on it ?

Symbol Selector : http://hpics.li/a01f07a

In the second windows, with the selection of the style, I would like to have in the empty windows all of my 12 new symbols. I would like to click on it to attribute them to a certain entities.

How can I make a new style with new symbols ?


texts: get value from attribute table and change symbology with QGIS and Python

texts:
I'm newbie developing plugins for QGIS using Python, this is the first one and it is very frustrating. Ok, this is the situation, I have a 1 point layer in QGIS. Its attribute table has 4 fields and 1 row:

ID|Estation|Data|Estate

1     1     20  "Green"

What I'm trying to do with my plugin is to get the string data from "Estate" field and check the string data:

If string is "Green", the one point layer symbology will turn to green.

If string is not "Green", the one point layer symbology will turn to red.

The last part of this plugin is to refresh every X seconds to check "Estate" field, and change the symbology of this point layer if the "Estate" field has changed.

This is what I have done so far in my python script, I put from the run method:

# run method that performs all the real work
def run(self):

    # create and show the dialog
    dlg = PlayMonitoringDialog()
    # show the dialog
    dlg.show()
    result = dlg.exec_()
    # See if OK was pressed
    if result == 1:
      # do something useful (delete the line containing pass and
      # substitute with your code
      def h():

        self.iface.mapCanvas().refresh
        layer = self.iface.activeLayer()
        Estado = layer.attributeDisplayName(3)
        if Estado == "green":
          props = { 'color' : '0,255,0' }
          sl = QgsSymbolLayerV2Registry.instance().symbolLayerMetadata("SimpleMarker").createSymbolLayer(props)
          s = QgsMarkerSymbolV2([sl])
          layer.setRendererV2(QgsSingleSymbolRendererV2(s))
        else:
          props = { 'color' : '255,0,0' }
          sl = QgsSymbolLayerV2Registry.instance().symbolLayerMetadata("SimpleMarker").createSymbolLayer(props)
          s = QgsMarkerSymbolV2([sl])
          layer.setRendererV2(QgsSingleSymbolRendererV2(s))
        if hasattr(layer, "setCacheImage"): layer.setCacheImage(None)
        layer.triggerRepaint()
      lanza()

  def lanza():
    t=threading.Timer(3.0,h)
    t.start()
  lanza()

It doesn't work, so I hope someone can help me. I have seen the small examples in http://www.qgis.org/pyqgis-cookbook, I have checked http://qgis.org/api/index.html and I have looked for examples on the Internet, but I still don't know how to solve the problems in my script.


texts: Reading File Geodatabase into Postgres with ogr2ogr

texts:
I have looked at other Q&As on the subject and none have solved my problem.

Reading Esri File Geodatabase (FileGDB) with OGR?
How to Import ESRI Geodatabase format .gdb into PostGIS

My problem is similar to the first referenced Q&A where my error says 'unable to open datasource'. I downloaded the test dataset in the solution and still did not work. My set up of ogr2ogr looks like this:

ogr2ogr -f "PostgreSQL" PG:"host=pg.rambck.local dbname=projects user=david password=xxxxxx" \\source\gdbname.gdb\datasetname -nln v1_dataset

When I run ogrinfo on the test dataset I get error: Unable to open datasource 'C:\david\trecks.gdb'

When I run ogr2ogr on my set up I get error: Unable to open datasource 'host=pg.rambck.local dbname=projects.......'

OpenFileGDB is listed in the list of drivers. What am I doing wrong?


texts: How to limit a specific area in order to have a topographic 3D map using MatLab?

texts:
I am using MatLab, I want to build a topographic 3D map, I am following this guide, I have been doing the steps listed on the section Draping Geoid Heights over Topography, and I do not have any trouble getting the expected result, it looks like this.

Now, my trouble is that I want to have a specific area, just a country not the whole world... How can I limit a specific area using coordinates in order to have a topographic 3D map that looks like this one.


texts: If else if in ArcMap Field Calculator (Python)

texts:
Help!  I think I am missing something that must be so obvious but I cannot figure out what I am doing wrong.  Here's some background:

I have one short integer field (MajAspNum) that has values 1-10.  I have another text field (Direction) that needs to be filled.  I am trying to use the Field Calculator to calculate it but am not succeeding in any way.  I should know what I'm doing but I just can't seem to make it work.

Pre-logic Script Code:

def des(num):
if num = 1:
    return "Flat"
elif num = 2:
    return "North"
elif num = 3:
    return "Northeast"
elif num = 4:
    return "East"
elif num = 5:
    return "Southeast"
elif num = 6:
    return "South"
elif num = 7:
    return "Southwest"
elif num = 8:
    return "West"
elif num = 9:
    return "Northwest"
elif num = 10:
    return "North"
else:
    return "N/A"

Direction =
des (!MajAspNum!)

It gives me a GeoProcessing error of invalid syntax on line 2.
I have tried just selecting by attributes in the attributes table and then using field calculator more simply on just those but that gives me a very weird error that literally makes no sense.  I'm hoping my error is a really obvious fix to someone out there!


texts: Trying to get all Landscan Urban areas in the USA (from Natural World shp files)

texts:
I'm desperately trying to get a json file of all the Landscan urban areas in the USA.

The ne_10m_urban_areas_landscan.shp has a list containing polygons of all urban areas in the world with a NAME attribute that lists the name of the urban area (i.e. 'Aarhus')

But in order to work out whether that NAME attribute is in the USA, I also need to combine it with data from ne_10m_populated_places.shp which has the name of the city along with the country its in (i.e. 'Aarhus' and 'Denmark').

However, whatever I do I just cannot get the data from both files.

ogr2ogr -f GeoJSON
  -sql "SELECT ne_10m_populated_places.*, ne_10m_urban_areas_landscan.*
  FROM ne_10m_urban_areas_landscan
  LEFT JOIN 'ne_10m_populated_places.shp'.ne_10m_populated_places
  ON ne_10m_urban_areas_landscan.name_conve = ne_10m_populated_places.NAME"
us-urban_areas.json ne_10m_urban_areas_landscan.shp

The above command gets me all the urban areas in the world joined with the populated places data, but its not WHERE country = USA.
If I do put in a WHERE clause then I get an SQL error because I'm filtering based on a field in the populated_places file, not the urbanareas one.

More promisingly...

ogr2ogr -f GeoJSON
  -sql "SELECT ne_10m_populated_places.*, ne_10m_urban_areas_landscan.*
  FROM ne_10m_populated_places
  LEFT JOIN 'ne_10m_urban_areas_landscan.shp'.ne_10m_urban_areas_landscan
  ON ne_10m_populated_places.NAME = ne_10m_urban_areas_landscan.name_conve
  WHERE ne_10m_populated_places.ADM0_A3 IN ('USA')"
us-urban_areas.json ne_10m_populated_places.shp

The above command returns just the cities in the USA, with all the joined data, except the single Point Geometries from the populated places file are kept and the Polygon geometries from the urban areas file are lost.

I just don't know what to do. Does anyone have any ideas?


texts: Why does a geoprocessing service fail every other time it is run?

texts:
I have a ArcGIS Server 10 geoprocessing service right now that takes in a string and runs a parameterized query to SQL Server 2008R2 through pyodbc. The service works, but only every other time. I can run it on my machine locally every time and also on the machine this runs on (also every time), but I can't run this through ArcGIS Server but every other time, either through the web interface or the ArcGIS JavaScript API.

Every other time, I get this error: Error Executing Task. Error processing message... along with a 500 error.

Is it maybe a problem with pyodbc? Or am I overlooking something small?

Here's what I'm working with (snipped SQL and a little extra). There's two input parameters, SearchTerm and SearchType (only can be "standard" in this snippet), and one output parameter. All parameters are strings.

import arcpy
import pyodbc
import sys
import urllib

connection = None
cursor = None

def search(search_term, search_type):
    try:
        global connection
        search_term = urllib.unquote_plus(search_term)
        connection = pyodbc.connect(...)
        if search_type == "standard":
            return string_search(search_term) # snipped out other use case.
        else:
            return "[]"
    except:
        print_message(sys.exc_info(), True)
        return '{"error": "' + str(sys.exc_info()) + '"}'

def string_search(search_term):
    global cursor
    cursor = connection.cursor()
    sql = """SELECT * FROM Resources WHERE ResourceName LIKE ?"""
    search_term = "%" + search_term.upper() + "%"
    cursor.execute(sql, search_term)
    return return_results()

def return_results():
    global connection
    global cursor
    return_ids = ""
    for row in cursor:
        resource_name = ""
        if row.ResourceName is None:
            resource_name = "Not Available"
        else:
            resource_name = row.ResourceName
        return_ids = return_ids + ',{"ResourceId":' + str(row.ResourceId) + '}'
    return_ids = return_ids[1:]
    cursor.close()
    connection.close()
    return "[" + return_ids + "]"

def print_message(message, is_error=False):
    message = str(message)
    if is_error:
        print "ERROR: " + message
        arcpy.AddError(message)
    else:
        print message
        arcpy.AddMessage(message)

if __name__ == '__main__':
    if arcpy.GetArgumentCount() > 0:
        return_text = search(arcpy.GetParameterAsText(0), arcpy.GetParameterAsText(1))
        arcpy.SetParameterAsText(2, return_text)
    else:
        arguments = ["chicken", "standard"]
        argv = tuple(arguments)
        print search(*argv)


texts: Python error while trying to use openlayers_plugin as basemap

texts:
While trying to use openlayerPlugin in QGIS as base map. I am getting below error

warning:C:/Users/*User_Id*/.qgis2/python/plugins\openlayers_plugin\openlayers_plugin.py:174: DeprecationWarning: QgsMapLayer.setLayerName() is deprecated
          layer.setLayerName(layerType.displayName)

    traceback:  File "C:/Users/*User_Id*/.qgis2/python/plugins\openlayers_plugin\weblayers\weblayer.py", line 100, in addLayer
        self._addLayerCallback(self)

..I have tried running it with/without proxy unchecked and getting same error both time.

Note: I dont have admin rights to the workspace system.


texts: Converting Excel Table to Geodatabase with Repeating Numbers Results in all Zeroes

texts:
I'm trying to convert an Excel file to a DBF file. I understand how to do it correctly using the Table to Table tool. However I am having an issue where a column of repeating numbers (labeled ZoneNum) are converted into nothing but zeroes. There doesn't seem to be any other issues with the other columns such as street names, ID number, work done, etc. I believe it has to do with the repeating numbers. For example 16 streets might be under ZoneNum 1 so their attribute would be 1. The next 7 streets might be under ZoneNum 3. I've tried to store the number as text instead but then nothing showed up. I've formatted each column as numbers or text. I've tried converting to a csv file instead to then convert to DBF file with no results. Please help!

Here's 3 sample lines from my file (minus the linebreaks):

104,14,Yosemite Drive,Everglade Circle,South End,Cape Seal

105,14,Zion Avenue,Jefferson Street,Carlsbad Street,Slurry Seal

106,17,Aberdeen Drive,Venice Avenue,West End,"1.5"" Edge Grind & Overlay"


texts:
Follow up question to QGIS Model Builder Script - Possible to merge layers in a folder?

I think I'm missing something very obvious here but can't seem to find it. After merging shapefiles in a folder, I would like a grid to be created. However, I am unsure how to call the result of the merged file into the grid function.

## Test=name
## Select_Directory=folder
## cellsize=number 1000
## grid=output vector

import os
import glob

os.chdir(Select_Directory)
output = glob.glob('./*.shp')
Layer_0 = output[0]

processing.runalg("saga:mergeshapeslayers",  Layer_0, ";".join(output))
processing.runalg("qgis:creategrid", cellsize, cellsize, 24108, 18351, 258293, 665638, 1, 'EPSG:7405', grid)


texts: How to download all OpenStreetMap highways in a large area?

texts:
I need to download all highways as line features from OpenStreetMap in a large area (about 10.000km²). Therefore I wanted to use overpass-turbo with the following statement in the wizard:

 highway=*

in my chosen bounding box.

Unfortunately it looks like my bounding box is too big, because I get the error:

 runtime error: Query run out of memory using about 2048 MB of RAM.

Do I have to change my overpass turbo query or is it not possible to download large datasets in overpass turbo? Which alternative method could I use to download all highways as line features in a larger region. I need to have all attributes, like the type of the highway, "foot=yes/no/...", "access="yes/no/....", etc. and if I use something like geofabrik, these attributes do not appear.


texts: Time series coloring of polygons in browser (preferably Mapbox)

texts:
I have a problem where I need to display a choropleth map with a time slider.  I have seen many examples of time slider implementations with CartoDB Torque as well as with mapbox, but I can't find any method of doing this for larger polygon datasets.

I know that I can add a layer for each of the points in the time series, then switch them on and off based on the position of a slider bar.  This falls apart for the following reasons:

1) As the number of time steps increases, the number of layers increases.  This creates some overhead in the map.

2) We currently create tilesets for our shape information so we can display the appropriate number of points on the polygon boundary at different zoom levels.  Because of this, the time series data is transferred multiple times as the user pans and zooms across the map.

There are a number of nearly relevant posts I've reviewed but none lead me to a solution where I can transfer the "metadata" file to the client once and then attach it to the visible polygons for styling.


texts: Finding intersection point of point and polyline using ArcMap?

texts:
I have two shapefiles that looks something like this:

One point file and one polyline file.

I want to find the points where the points and line would intersect if you drew a perpendicular line between them (the shortest distance). The next step is to mark these points on the polyline or a rasterized version of the polyline, but first I need to find a way to define them.


texts:
Using QGIS 2.12.3, within the Field Calculator "Geometry" section, there are certain functions, such as "area", that also have a corresponding term beginning with a dollar sign, i.e. "$area".  I don't understand the difference between the two.

For example, in the case of "area", the Field Calculator's description pane states: "Returns the area of a geometry polygon object", while the "$area" description is "Returns the area of the current feature."  I confess that I don't understand the difference, nor how to use them in a formula.


texts: How do I add someone else's Tileset in Mapbox Studio?

texts:
I'm trying to see the Microsoft Building Footprints in my hometown with the Tileset Chetan has prepared but his Tileset page doesn't have the Map ID that my own Tilesets have and is required by Studio. When I'm in Studio, I enter the url for the tileset but that doesn't come up with anything.

Chetan says " If anyone wants to visualize the data in Mapbox Studio, you can use this tileset*." but there is no explanation of how to use his tileset.


texts:
I've been let loose in the workplace to learn python to do things in Arcmap 10.
so, I am learning python as I go and trying to remember the programming I have done.

Where I am in this project is converting a dbf, or csv, xls in a simple fashion.
from there, all the files will be copied together into one file.
I've got the all-in-one xls working, but I can't find an easy, simple dbf to xls solution.

I condensed code found here:
http://blog.gmane.org/gmane.comp.python.education/page=12

into:

from xlwt import Workbook
import dbfpy.dbf

input0 = '...file.dbf'
output = '...file.xls'

def test1():
    dbf = dbfpy.dbf.Dbf(input0)
    book = Workbook()
    sheet1 = book.add_sheet('Sheet 1')
    for row in range(len(dbf)):
        for col in range(2):#chop to two for my purposes, gm
            sheet1.row(row).write(col, dbf[row][col])
    book.save(output)

test1()

This works, minus the lack of field names.


texts: ESRI FIleGDB API dll failure: ogr2ogr can't recognize File GDB as driver type

texts:
I have been trying this for a day or two now, and I have followed the instructions listed here:
 How to add support for FileGDB (Esri file gdb API) driver in fiona

but for the life of me I can't seem to get ogr2ogr to recognize FileGDB as a driver. My end goal is to append shapefiles to an already exisiting .gdb file on my local machine. I have set all the environment variables to the seemingly appropriate locations (GDAL_DATA, GDAL_DRIVER_PATH, PROJ_LIB, etc.), and all my dll files seem to be in the correct location. I am using Python 3.7.4 via Anaconda3 command console. Furthermore, I am using Windows 10.


texts: How to install a QGIS plugin when offline?

texts:
Due to various IT policies at my workplace, QGIS is installed on a machine that is not connected to the internet. I wish to install a couple of QGIS plugins on this system.

I have downloaded the required plugins from http://pyqgis.org/repo/contributed. How do I install them in QGIS?


texts:
I'm trying to develop gis based web application. I'm using leaflet for frontend  and django for backend. My requirement is to plot raster data in leaflet.  I'm using rasterio for reading raster file. And here is my View to read raster file upload by the user
class UploadRasterView(View):
    def post(self, request):
        context = dict()
        raster = RasterFileUploadForm(request.POST, request.FILES)
        folder_name = 'raster_data'
        if raster.is_valid:
            file_name = ''
            for _file in self.request.FILES.getlist('file'):
                handle_uploaded_file(_file, folder_name)
                if _file.name.split('.')[-1] == 'tif':
                    file_name = _file.name
            file_location = settings.MEDIA_ROOT + '/' + folder_name + '/' + file_name
            data = rasterio.open(file_location)
        print(type(data))
        return JsonResponse(context)

Up to here everythings is working fine but i didn't know to serialize this data, os data i can plot raster data in leaflet.


texts: Aquisition date of digital globe imagery, Image Metadata File

texts:
It might sound strange but I can't seem to find the acquisition date of images I acquired from digital globe. In the metadata file I only find the date the images were sent to me. I am using QGIS.

Where do I need to look for the acquisition date?


texts:
My guess is that I'm doing something wrong with QGIS, but here goes.

I have multiple point data layers in WGS84, mostly of UK locations, but with point in Canada, South Africa and New Zealand. I also have line data generated from the point data with MMQGIS/Hub Lines. What I need are accurate distances (line lengths) for the Hub Lines. I have tried this in QGIS 2.12.2 in Windows 32-bit, and in 2.16.1 in Kubuntu 64-bit. The processor is an AMD A-series quad-core, and 8GB RAM.

First, I have remembered to change to an equidistant projection to correct the distortion in WGS84, and have tried using On the Fly re-projection, and simply starting a new project with minimal point and line data all converted to, for example, World Azimuthal Equidistant, but I am getting severe errors. To simplify my explanation I have also tried the QGIS tutorial on calculating line lengths here:
http://www.qgistutorials.com/en/docs/calculating_line_lengths.html

Following these tutorial instructions gives me a result for uident: 1506 (first line of attribute table) of 0.01 km; most other entries are 0.00 km. The correct result appears to be 637.07 km. If I calculate another column from this data in meters, without the division by 1000, the same entry shows the line to be 8.40 m in length. Further, if I load a Natural Earth country boundaries map to this project and re-project, the line measurement tool shows the whole of mainland USA (excluding Alaska) to be roughly 50 km east-west, and 25 km north-south. I am getting similarly wild results in my own data.

What have I got wrong?

EDIT: Turning Linux off and turning it back on again has cured the issue in the tutorial mentioned above; but why? And why do I still have the error in my own data? In the latter the distance from Bristol, UK to Auckland, New Zealand is calculated at 0.154 km, and the measuring tool shows the same.

I am still getting something wrong. If I now add the Natural Earth world coastline layer to the project with On the Fly project active, all my point and line data shrinks into a spot at coordinate 0,0.


texts: Krovak projection EPSG 5514: Correct data transfer between R and ArcGIS?

texts:
I have a problem which many ArcGIS, QGIS and R users from Czech Republic and Slovakia may encounter. In ArcGIS, I am working EPSG: 5514,  Projected Coordinate System, defined as:

S-JTSK_Krovak_East_North
WKID: 5514 Authority: EPSG

Projection: Krovak
False_Easting: 0.0
False_Northing: 0.0
Pseudo_Standard_Parallel_1: 78.5
Scale_Factor: 0.9999
Azimuth: 30.28813975277778
Longitude_Of_Center: 24.83333333333333
Latitude_Of_Center: 49.5
X_Scale: -1.0
Y_Scale: 1.0
XY_Plane_Rotation: 90.0
Linear Unit: Meter (1.0)

Geographic Coordinate System: GCS_S_JTSK
Angular Unit: Degree (0.0174532925199433)
Prime Meridian: Greenwich (0.0)
Datum: D_S_JTSK
  Spheroid: Bessel_1841
    Semimajor Axis: 6377397.155
    Semiminor Axis: 6356078.962818189
    Inverse Flattening: 299.1528128

When I load my file in R, my proj4string(data_in_5514) shows:

 "+proj=krovak +lat_0=49.5 +lon_0=24.83333333333333 +k=0.9999 +x_0=0 +y_0=0 +ellps=bessel +units=m +no_defs"

I've got the same information when I export my files from R using writeRaster() function from raster package. However, when I want to open my data modified in R again in ArcGIS, my ArcGIS crashes, and if I am lucky, I've got a warning message :

Warning, inconsistent extent!

which means (http://support.esri.com/em/knowledgebase/techarticles/detail/37732):

  This error message is returned if the data has an incorrect projection definition.

  Specifically, the data is in a projected coordinate system with units of feet, meters or other linear units, but the projection has been defined with a Geographic Coordinate System (GCS) that has units in decimal degrees.

Thus, obviously the projection of my data has changed, but when and how ???

How can I correctly transfer data between R and ArcGIS using Krovak EPSG 5514? I've also found, that this one PCS is also listed as EPSG: 102067 under ArcGIS, however, R doesn't recognize this code.

EDIT :

sample data: http://ulozto.sk/xNHrBv2D/r3-zip

as an output I was using several formats. I was able to opet Gtiff once, .img neidther one time, .grd seems not to be directly recognized by ArcGIS.

export data:

writeRaster(aaa, filename="r3.img", format="HFA", overwrite=TRUE)
# export data
writeRaster(aaa, filename="r3.tif", format="GTiff", overwrite=TRUE)
# export data
writeRaster(aaa, filename="r3", overwrite=TRUE)  # creates a .grd


texts: Determine Geometry of Multiple Layers in QGIS

texts:
I am sure this is a newbie question.

I am importing a large number of layers into QGIS from a complex .kml file and need to combine into a single layer and finally save as a .shp file.

I repeated run into the incompatible geometry error when I attempt to combine layers in Qgis.  I can check each layer to determine its geometry and change this, but given I am dealing with 100 or so layers, I am looking for a more efficient method (perhaps python query?) to identify all the layers'
geometries or other more automated method.


texts: How to replicate a polygon over several points with different locations in Arcgis?

texts:
I have two shapefiles, one is a polygon and the other is points.
I would like to create a new shapefile which has the shape of the polygon in all the places that there is points. I would like to replicate the polygon under the points.

1= the polygon that I want to replicate. 2= the points represents the location where the replicated polygons should be (one polygon per point).


texts: After DXF conversion, layer with same coordinate and projection not visible

texts:
I've converted a DXF file to shp (using another DXF importer/convertor plugin).
In QGIS this shp will open fine if standalone but is not visible when added to a project with other shapefiles.
The xy coordinates are comparable points are the same, as is the CRS.
I've tried various options with no success - any thoughts on the cause?


texts: Model Builder - Using Only Certain Fields

texts:
I have to use ModelBuilder in ArcGIS 10.1 to create a mix of two tools: buffer and select by location. However, I'm supposed to let the user only enter certain information, and have other fields be assumed. For example, I only want the user to fill in the first three fields, and have the Side Type, End Type, and Dissolve Type be the same values all the time:

How would I go about doing this is ModelBuilder?


texts: How to fill sinks in DEM with specific amount of water

texts:
I identified the sinks in my terrain based on the DEM (grid) by using the 'fill'-tool and substracting the original DEM from the filled DEM. I also identified watersheds as well as flow paths.

Now i want to model which water level each sink will have in case of a specific rain event. Can any of you recommend a understandable tutorial how to do that in ArcGIS 10.2 or explain the procedure for me?


texts: Heatmap of sold homes with constant cell size?

texts:
QGIS version 2.18.5

I'm working with annual house prices in England and Wales with excel columns as: A; Transfer Price B; Postcode C; Lat D; Long. There is approx 400,000 rows of data for 2011.

I'm playing around with the colour ramp editor with my points classified by house sale price (transfer price). The highest sale price for 2011 sits at £21,000,000. I've set my gradient colour ramp to 'discrete' with blue up to 99.9% with the remaining 0.1% as red. I've also set the maximum value to 21,000,000.

You would think that would single out the house worth 21mill, but I've attached a screenshot of my heat map below which apparently indicates there were multiple billionaires who sold their mansions last year. Can someone explain why this is happening?

I also notice when I zoom out, the red areas just get bigger! Why don't they stay the same size? Clearly I don't understand how a heatmap works...


texts: Flow length analysis only shows pixel value, how to get length?

texts: When I create a flow length analysis, and use identify tool on a cell in the layer, it just shows "stretched value" and "pixel value" I don't understand how it's supposed to show length. I'm using ARCGIS (ARCMAP) current version.


texts: Converting table drawn as annotation in AutoCAD to format suitable to link attributes to cables in ArcGIS for Desktop?

texts:
I am working on cables and each cable has attributes.

I want to link these attributes to the cables in ArcGIS for Desktop.

When I open these cables in AutoCAD it appears that their attributes are drawn as annotation. I would like these attributes to be in a table.

Later I need to export the data to an Excel spreadsheet.

There is an example of the Table open in AutoCAD below.


texts: Convert an image to an array or a list to train a supervised classifer in Google Earth Engine

texts:
I have two GeoTiff images that will represent training datasets classes that will be fed to a supervised classifier.

As opposed to this classification example where the geographic regions are extracted from a fusion table and the overlaid on the image, I don't have a fusion table so I'm going to use an image with its values directly as a training dataset.

My question is how can I transform an Image object into an array (or a FeatureCollection if I follow the example linked above)?

Edit:

I've managed to convert my two training datasets rasters into a KML vector, and from there I've uploaded them to Google Fusion Tables (See in the code below class1 and class2). The two raster images to train and classify (i.e. image1 and image2 in the source code) were uploaded as assets; while the first one was used to train the classifier, the second one has to be classified. The problem with the code below is that all the pixels in the second image were classified in the same class (class = 0):

// load vector masks from fusion tables
var class1 = ee.FeatureCollection('ft:xxx');
var class2 = ee.FeatureCollection('ft:yyy');
var classes = class1.merge(class2);

// load training image and the image to classify
var image1 = ee.Image('aaa').select('b1');
var image2 = ee.Image('bbb').select('b1');

// train the classifier
var training = image1.sampleRegions({
    collection: classes,
    properties: ['class'],
    scale: 1
});
var classifier = ee.Classifier.svm().train(training, 'class', ['b1']);
var classified_image2 = image2.classify(classifier);

// show classified image
Map.addLayer(classified_image2, {}, 'Classified image');


texts: Random pixel shift when adding WMS layer

texts:
I made a WMS service with MapServer to display a raster data. When I add it as an ArcGIS or a QGIS layer, there is a shift between the original layer and the WMS layer. This should not be a projection issue because the maps and the layer are all in the same coordinate system.

This shift is not static (if I slightly move the map, it will change in X and Y) and less than one pixel. Even if it is not a problem for many application, I need a very precise location of my pixels. I would therefore like to avoid this random shift. I've tried to remove the shift by setting (or not) the extent and the size of the map/layer in my mapfile, but nothing worked. The problem is exactly the same with ArcGIS and QGIS, therefore I am convinced that it comes from mapServer.

  MAP

  IMAGETYPE      png

  SIZE           120197 70397

  EXTENT      50184.875 26720.668 289726.486 167015.441 # EXTENT LB72 in meters, from gdalinfo

  PROJECTION

"init=epsg:31370"

  END

  LAYER # Raster layer

NAME         test

STATUS       ON

TYPE         RASTER

METADATA

  "wms_title"    "test"

END

PROJECTION

  "init=epsg:31370"

END

EXTENT    50184.875 26720.668 289726.486 167015.441 # EXTENT LB72 in meters, from gdalinfo


texts: PostGIS - Measuring the angle at the touching point of two lines

texts:
I want to calculate the angle of two lines of a routing network exactly at the "touching point" (the topological vertex).
My previous attempt was using the ST_ANGLE() function of PostGIS, which gave me a "wrong" result of 76° for the point 4060 between line [16543,4060] and line [4060,3382], since ST_ANGLE() takes the beginning and end node of the input lines to measure the angle betwen those two.

ST_LineInterpolatePoint to generate auxiliary points close to the node isn't that easy also due to the different orientations of the edges.


texts:
I am new to the ArcGIS.

Can any one please refer a beginners tutorial for this tool.

Referred the help, but couldn't able to get the basics to kick start.


texts: Is there a symbols level option/solution for SLD?

texts:
Is it possible to have an sld blend the symbology in the lines together just as QGIS and ArcGIS allow using symbol levels?

In QGIS I can use symbol levels on the line layer and have this appearance:

currently in sld it has this appearance:

I was hoping to have the appearance of the individual lines blending together as one. The intent is to apply the sld to a geoserver layer. Currently I am only using rules like the one listed below.

    <Rule>
      <MinScaleDenominator>2500</MinScaleDenominator>
      <MaxScaleDenominator>12500</MaxScaleDenominator>
      <LineSymbolizer>
        <Stroke>
          <CssParameter name="stroke">#ffe14c</CssParameter>
          <CssParameter name="stroke-width">3</CssParameter>
          <CssParameter name="stroke-linejoin">round</CssParameter>
          <CssParameter name="stroke-linecap">round</CssParameter>
        </Stroke>
      </LineSymbolizer>
    </Rule>


texts: Problem with interpolating a vertical profile in 3D Analyst (ArcMap 10.1)

texts:
I have a problem when I want to create a vertical profile in 3D Analyst from a DEM file. I want to interpolate the data and show points from 20 to 20 meters, as I specified in the option panel (see figure 1). After I press the button to show the Profile Graph, the X-axis do not show values between 20 meters as I specified in the option panel (see figure 2). What can I do to solve this issue?

Note: the Y-axis is displayed correctly in terms of altitude. Datum: WGS_1984


texts: Seeking command line tool for converting projection in .prj and GeoTiff to EPSG codes

texts:
Is there a command line tool to convert a shapefile projection in the .prj file into corresponding EPSG code?

I know that there is website http://prj2epsg.org/search for doing this online, but I wanted a command line tool to use in shell scripts (for things like shp2pgqsl) without involving internet.

Also, is there such a tool for extracting EPSG code from the auxiliary information in a GeoTiff file?


texts: I have coordinates in decimal degrees and need to calculate area for polygon objects. I tried to use AreaCalculator and set up criteria to find all these bigger than 100 sqm. But result is not as I expected. Should I convert first from decimal degree to metric?


texts: Analyse driving pattern from automotive dataset

texts:
I am trying to derive driving patterns from our automotive dataset. Our dataset has vehicle names along with their GPS points. One of the driving pattern I have read about is “hub and spoke” in which vehicles start and stop at the same location. Other vehicles like trucks will travel in the nights following a different pattern.

From my research, I understand that Kmeans clustering of data will not work on this data, as different vehicles might have different patterns. Just wondering whether any other ML algorithms would be effective in solving this problem? Because I have 1 million vehicles, if i consider all GPS points in a vehicle, that would be a massive dataset. Would it be better if i could restrict the GPS points based on some condition.


texts: Trying to build lines by cursing through points

texts:
I am trying to create a tool that creates lines from selected point features.  I am able to select the point features and create a search cursor to curse through them.  I can message box attributes from each point.  I am having trouble creating an array of points by getting the XY coordinates from each row in the search cursor.

Below is the code:

import arcpy
import pythonaddins

import arcpy.mapping as mapping
mxd = mapping.MapDocument("CURRENT")
layers = mapping.ListLayers(mxd)
# for lyr in layers:
    # print(lyr.name)

class ToolClass2(object):
    """Implementation for Fishnet_addin.tool (Tool)"""
    def __init__(self):
        self.enabled = True
        self.cursor = 3
        self.shape = 'Rectangle'
    def onMouseDown(self, x, y, button, shift):
        pass
    def onMouseDownMap(self, x, y, button, shift):
        pass
    def onMouseUp(self, x, y, button, shift):
        pass
    def onMouseUpMap(self, x, y, button, shift):
        pass
    def onMouseMove(self, x, y, button, shift):
        pass
    def onMouseMoveMap(self, x, y, button, shift):
        pass
    def onDblClick(self):
        pass
    def onKeyDown(self, keycode, shift):
        pass
    def onKeyUp(self, keycode, shift):
        pass
    def deactivate(self):
        pass
    def onCircle(self, circle_geometry):
        pass
    def onLine(self, line_geometry):
        pass
    def onRectangle(self, rectangle_geometry):
      """Occurs when the rectangle is drawn and the mouse button is released.
    The rectangle is a extent object."""

  extent = rectangle_geometry
  mxd = mapping.MapDocument("CURRENT")
  df = arcpy.mapping.ListDataFrames(mxd)[0]
  layer = mapping.ListLayers(mxd)[0]
  # Select the first layer in the project, the point layer.

  # Create a fishnet with 1 row and 1 column.
  if arcpy.Exists(r'in_memory\fishnet'):
      arcpy.Delete_management(r'in_memory\fishnet')
  fishnet = arcpy.CreateFishnet_management(r'in_memory\fishnet',
                          '%f %f' %(extent.XMin, extent.YMin),
                          '%f %f' %(extent.XMin, extent.YMax),
                          0, 0, 1, 1,
                          '%f %f' %(extent.XMax, extent.YMax),'NO_LABELS',
                          '%f %f %f %f' %(extent.XMin, extent.YMin, extent.XMax, extent.YMax), 'POLYGON')
  fishnet.visible = False
  # Select by location the point layer that intersects the fishnet
  arcpy.SelectLayerByLocation_management(layer, "INTERSECT", fishnet, "", "NEW_SELECTION", "NOT_INVERT")
  #  Count the number of points selected
  selectionCount = arcpy.GetCount_management(layer)
  count = int(selectionCount.getOutput(0))
  # if 2 points are selected create a search cursor for them
  if (count == 2):
      pointCursor = arcpy.SearchCursor(layer, ["SHAPE@XY", "NAVD88"])

      messageResult = pythonaddins.MessageBox(selectionCount,"Points Selected Event",3)
      array = []
      for row in pointCursor:
        NAVD = row.getValue("NAVD88")
        messageResult = pythonaddins.MessageBox(NAVD,row,3)
        # attempt to build an array of points from the points in the search cursor.  This is the part that does not work.  Comment this line and uncomment the array code below and it works.
        array.add(arcpy.Point(row.getPart(0)))

  # This is an empty line layer with one field called "slope"  I manually start editing on this layer before using the tool.
  fc = r"E:\ArcGIS\Working.gdb\testLines"
  #  Create an insert cursor to add the line feature to the line layer.
  cursor = arcpy.da.InsertCursor(fc, ["SHAPE@","slope"])

  #  The commented code below works.  What I am trying to do is to build this 'array' from the XY coordinates of the points in the search cursor.
  # array = arcpy.Array([arcpy.Point(-82.226344, 26.915444),
                 # arcpy.Point(-82.226499, 26.915005)])

  spatial_reference = arcpy.SpatialReference(4326)
  polyline = arcpy.Polyline(array, spatial_reference)

  theslope = 2.08

  cursor.insertRow([polyline,theslope])

  # arcpyRefreshTOC()
  arcpy.RefreshActiveView()
  return fishnet


texts: Vector Layer not displayed properly in 3D

texts: I am using Qgis2threejs plugin to display my map in 3D. But I am facing a problem in displaying vector layer properly. Please look at the following image.Buildings are not displayed according to DEM. it looks like they are floating in the air. Also Roads and Buildings are distorted from their original places. How to fix this? Can someone help?


texts:
I am new to Python and am having difficulties trying to loop multiple rasters into a conditional statement.  I want to extract all land-uses that have the value 21 from a folder with three different land-use layers.

When I run my code, I get an output that shows the land-use that equals value 21 for the first raster file in my folder, but it does not calculate an output for my other two raster files in my folder.  My code is:

    import arcpy
    from arcpy.sa import *
    arpy.CheckOutExtension('Spatial')
    arcpy.env.workspace=r'U:data\landuses'
    arcpy.env.scratchWorkspace='U\data\work.gdb'
    rasters=arcpy.ListRasters("*")
    for raster in rasters:
    ...     ras=Raster(raster)
    ...     outraster=Con(ras==21,21)
    ...     outraster.save('U:\data\work.gdb' + raster)

When I run this, I get:

  ERROR 010093: Output raster format UNKNOWN is unsupported.

My raster files are all "GRID" format.


texts: How to merge polygons from two different features without creating a hole in QGIS?

texts: How is it possible to merge the yellow polygon with the blue polygon without creating a hole and keeping the big yellow polygon intact. Each time I try to delete the small yellow polygon, the whole thing vanishes!


texts: How to compile QGIS master snapshot from a certain point in history?

texts:
I would like to find out which commit in QGIS code caused a certain change in behaviour. Since the archives at http://download2.osgeo.org/osgeo4w/release/qgis/qgis-dev/ are not complete, I think that the best would be to compile the source code and generate a snapshot of qgis-dev at certain commits that are suspected to have caused the change. How should I proceed, if I want to try e.g. to make a binary of QGIS based on the following snapshot tree: https://github.com/qgis/QGIS/tree/0b9d1dc49b?

Both, Linux and Windows suggestions, are welcome.


texts:
I am trying to access a VectorTile layers attributes in Openlayers3 like this example. The layer is being server by 'Chubbs'

The layer:

var observationLayer = new ol.layer.VectorTile({
      source: new ol.source.VectorTile({
        format: new ol.format.MVT(),
        tileGrid: ol.tilegrid.createXYZ({maxZoom: 22}),
        tilePixelRatio: 16,
        url: '/services/postgis/test_observations/geometry/vector-tiles/{z}/{x}/{y}.pbf'
      }),
      style: observationStyle
    });

Trying to get the attributes when the mouse hovered over the point:

function onMouseMove(browserEvent) {
        var coordinate = browserEvent.coordinate;
        var pixel = map.getPixelFromCoordinate(coordinate);

        map.forEachFeatureAtPixel(pixel, function(feature, layer) {
          // Small sub set of all the  things I have tried...
          console.log('feature:');
          console.log(feature); // returns the features!

          console.log('layer:');
          console.log(layer); // returns the layer!

          console.log('layer.getProperties():');
          console.log(layer.getProperties()); // returns the prope

          var source = layer.getSource();
          console.log('source:');
          console.log(source);

          var properties = source.getProperties();
          console.log('properties:');
          console.log(properties); // returns an empty object
        });
      }
      map.on('pointermove', onMouseMove);


texts: Changing KML markers/Icons in Google Earth Pro?

texts:
How to change KML markers/Icons in Google Earth Pro?

I have put in all of the paths that I could find on kml4earth.appspot.com/icons.html#notes?

I am not understanding the procedure.


texts: Python Widget Wrapper to extend basewidget class in ArcGIS Viewer for Flex?

texts:
I'd like to create a flex wrapper for python code run locally (utilizing arcgis data and outputting data that can be displayed on the flexviewer).  I'd like to use gdal or pyproj modules to process data locally with distributed scripts.

How would one go about creating a wrapper-widget that extends the basewidget class in the arcgis viewer for flex that can communicate with the local install of the python interpreter?

web.py rest services

in terms of the workflow:
user without arcgis server license (for geoprocessing tasks) wants to perform some spatial analysis on a graphic layer in the flex viewer (flash or air).

python script running as a web service (web.py?) accepts vector input ("public" viewshed geoprocessing task demonstrated by esri (http://sampleserver1.arcgisonline.com/ArcGIS/rest/services/Elevation/ESRI_Elevation_World/GPServer/Viewshed) creates a graphic polygon layer as output).  Using a graphic layer from the flex viewer as input to a "point in polygon" python script (http://geospatialpython.com/2011/01/point-in-polygon.html) would make the output from the viewshed (or other graphic generating tool) more useful.


texts:
I'm a complete novice to GIS and just began an Introductory course. I am having an issue when I try to access the Symbology tab on a layer. Specifically, when I click on the 'Categories' field and 'Unique Values' subfield, ArcMap crashes and comes up with a report window to send to ESRI. The message is

  ArcGIS for Desktop has encountered a serious application error and is
  unable to continue.

What could be the problem?


texts: Copy Paste of Feature in Vector Layer Throws Error in QGIS 2.18.10

texts:
I just upgraded to latest release of QGIS (QGIS version 2.18.10 64bit, QGIS code revision 59e0f78). When I try to copy/paste a feature in a vector layer, I get:

  Could not commit changes to layer reader

  Errors: ERROR: 1 feature(s) not added - geometry type is not
  compatible with the current layer.

Copy/paste worked fine with 2.18.9 and previous versions. Anyone else seeing this issue?


texts: Finding closest line given 2 points using PostGIS?

texts:
I have a table t that contains a column line_positions which is of type line.
Given 2 points I want to find the closest line that is close enough (less than 10km) and that does not pass too close to a point I want to avoid (20km minimum).
Currently I use

SELECT t.*
FROM path t
WHERE
  ST_DWithin(ST_GeographyFromText('Point(69.835 22.596)'), t.line_positions, 10000, FALSE)  AND
  ST_DWithin(ST_GeographyFromText('Point(69.856 22.519)'), t.line_positions, 10000, false) AND
  NOT ST_DWithin(ST_GeographyFromText('Point(-79.804 9.141)'), t.line_positions, 20000, false)
ORDER BY
  ST_Distance(ST_GeographyFromText('Point(69.835 22.576)'), t.line_positions, false) +
  ST_Distance(ST_GeographyFromText('Point(69.856 22.519)'), t.line_positions, false)
  ASC
LIMIT 1

There is a gist index ix_path_line_positions on the line_positions column.

It works but slow, between 3s and 30s for only 100000 rows in t.

explain analyze gives :

Limit  (cost=9.95..9.95 rows=1 width=1432) (actual time=21729.253..21729.254 rows=1 loops=1)
   ->  Sort  (cost=9.95..9.95 rows=1 width=1432) (actual time=21729.251..21729.251 rows=1 loops=1)
         Sort Key: ((_st_distance('0101000020E61000003D0AD7A370755140FA7E6ABC74933640'::geography, line_positions, '0'::double precision, false) + _st_distance('0101000020E6100000105839B4C8765140BE9F1A2FDD843640'::geography, line_positions, '0'::double precision, false)))
         Sort Method: top-N heapsort  Memory: 26kB"
         ->  Index Scan using ix_path_line_positions on path t  (cost=0.28..9.94 rows=1 width=1432) (actual time=93.490..21710.562 rows=690 loops=1)
           Index Cond: ((line_positions && '0101000020E61000003D0AD7A3707551407F6ABC7493983640'::geography) AND (line_positions && '0101000020E6100000105839B4C8765140BE9F1A2FDD843640'::geography))
           Filter: (('0101000020E61000003D0AD7A3707551407F6ABC7493983640'::geography && _st_expand(line_positions, '10000'::double precision)) AND ('0101000020E6100000105839B4C8765140BE9F1A2FDD843640'::geography && _st_expand(line_positions, '10000'::double precision)) AND _st_dwithin('0101000020E61000003D0AD7A3707551407F6ABC7493983640'::geography, line_positions, '10000'::double precision, false) AND _st_dwithin('0101000020E6100000105839B4C8765140BE9F1A2FDD843640'::geography, line_positions, '10000'::double precision, false) AND ((NOT ('0101000020E6100000FA7E6ABC74F353C0D578E92631482240'::geography && _st_expand(line_positions, '20000'::double precision))) OR (NOT (line_positions && '0101000020E6100000FA7E6ABC74F353C0D578E92631482240'::geography)) OR (NOT _st_dwithin('0101000020E6100000FA7E6ABC74F353C0D578E92631482240'::geography, line_positions, '20000'::double precision, false))))
           Rows Removed by Filter: 15365
Planning time: 0.491 ms
Execution time: 21729.321 ms

How could I improve it ? Using geometry calculation instead (but my track could span several thousands km, will the computed distances correct) ? Using <-> KNN operator (but since I order on the sum of 2 distances, it does not seem to use the gist index anyway) ?


texts: Connect to a specific WMS layer in QGIS (using Python or XML)

texts:
I would like to create a short list of useful WMS layers for non-GIS specialists colleagues. The idea is to give them a document allowing them to configure their connections easily - either a XML that they can charge using "add wms layer", or a short python script to copy/paste in the console. Is there a way to achieve that in QGIS?

So far, I've only been able to either connect to a WMS service using XML, without being able to choose only the layers that my colleagues will actually use :

<!DOCTYPE connections>
<qgsWMSConnections version="1.0">
    <wms ignoreGetMapURI="false" smoothPixmapTransform="false" dpiMode="7" password="" ignoreGetFeatureInfoURI="false" referer="" username="" url="http://wms.craig.fr/mnt" invertAxisOrientation="false" ignoreAxisOrientation="false" name="[CRAIG] Modèles Numériques de Terrain : ombrage, couleur, courbes de niveau"/>
</qgsWMSConnections>

or to open specific layers using Python (see here for example) - but I don't want the layers to actually open, I only want them to show in the browser panel...


texts:
I'm using a unique list to populate up to 10 fields for a fc with correct geometry, in an "overloading" manner as you would in excel. I noticed on my first run through I had an "index out of range" error, which was strange as values should not populate more than 10 fields per record. I did a bit of an error check with a "break" and print, and it seems from the final result that there are duplicate values being entered where there most definitely should not be. Problem code:

for x in uniqueList:
    query = "[BAND_NUM] = \"%s\"" % x
    if x <> "" and x <> " " and x is not None:
        try:
            gp.MakeFeatureLayer("FIRST", "lyr1")
            gp.MakeFeatureLayer("END", "lyr2")
            gp.SelectLayerByAttribute("lyr1", "NEW_SELECTION", query)
            gp.SelectLayerbyLocation("lyr2", "HAVE_THEIR_CENTER_IN", "lyr1", "", "NEW_SELECTION")

            cursor = gp.UpdateCursor("lyr2")
            urow = cursor.next()
            count = 0
            print x

            # update NULL or blank fields
            while urow:
                # this was my error check
                if count == 10:
                    print count
                    break

                if urow.getValue(fieldvalues[count]) is not None and urow.getValue(fieldvalues[count]) <> "":
                    count += 1
                # go to next feature only after an update
                else:

                    urow.setValue(fieldvalues[count], x)
                    cursor.updateRow(urow)
                    count = 0
                    urow = cursor.Next()

            gp.SelectLayerByAttribute("lyr1","CLEAR_SELECTION")
            gp.SelectLayerByAttribute("lyr2","CLEAR_SELECTION")
        except arcgisscripting.ExecuteError:
            gp.AddError(gp.GetMessages(2))
        finally:
            gp.Delete("lyr1")
            gp.Delete("lyr2")
            del urow, cursor

del gp

The list looks like - fieldvalues = ['FN1', 'FN2', 'FN3', 'FN4', 'FN5', 'FN6', 'FN7', 'FN8', 'FN9', 'FN10']

I believe the problem is in the way I'm trying to loop this - my first if statement has no cursor.next() because I want to stay on the same record, at a higher index. The logic is that if the record gets updated, go to the next one in the cursor. I had something similar before where I would append one field in a comma delimited manner, but it appears I won't be using a look-up table so 10 fields is the way I'm going.

Here's an example of the duplication -


texts: Is there ArcGIS (ideally) solution for spill modeling from particular pour point on crude oil pipeline?

texts:
wondering if anyone knows of a solution for spill modeling from a particular pour point on a crude oil pipeline? I currently have all my data within ArcGIS (DEM, major/moderate/minor watercourses, drains, embankments) so would ideally be looking for an ESRI solution, but I'm open to other options.

I did look into ArcHydro, but as far a I can tell wont let me input different variables such as spill volume and oil/soil density. I ideally would also need a model that calculates losses from oil adhesion to land over the oiled path, oil pooling in large depressions on the land surface, and oil evaporation to the atmosphere.

Thanks!


texts: Error when trying to save an Excel file from attribute table using QGIS

texts:
I keep getting this error (see below) when trying to save the attribute table to an Excel file using XY Tools. Can anyone help? I have deleted and re-installed xytools but this hasn't helped.

  An error has occured while executing Python code:

Traceback (most recent call last):
File "C:/Users/Doug/.qgis2/python/plugins\xytools\xytools.py", line 301, in excelSave
xlw.writeAttributeRow(rowNr, values)
File "C:/Users/Doug/.qgis2/python/plugins\xytools\providers\excel.py", line

  57, in writeAttributeRow
      self.ws.write(rowNr, colNr, cell)
      File "build\bdist.win-amd64\egg\xlwt\Worksheet.py", line 1030, in write
      self.row(r).write(c, label, style)
      File "build\bdist.win-amd64\egg\xlwt\Worksheet.py", line 1078, in row
      self.rows[indx] = self.Row(indx, self)
      File "build\bdist.win-amd64\egg\xlwt\Row.py", line 42, in __init
      raise ValueError("row index (%r) not an int in range(65536)" % rowx)
      ValueError: row index (65536) not an int in range(65536)

  Python version:
  2.7.5 (default, May 15 2013, 22:44:16) [MSC v.1500 64 bit (AMD64)]

  QGIS version:
  2.6.1-Brighton Brighton, e2a51df


texts: Counting points around a point conditional on time values

texts:
I want to count the number of points from one layer around each point from a different points layer. The thing is that each point (from both layers) has an "entrance date", so I want to count only those points that have an entry date prior to the one relating to.

For example:

I have a point with an entry date equal to August 2016, and a second layer of 2 points with an entry date equal to July 2016, and 2 points with an entry date equal to December 2016. All of the points are within 100m of the first point.
I want to be able to get the number of points, within 100m radius (in this example), that have an entry date prior to August 2016. In this example, the result would be 2.


texts: Earth Engine chart image series by region tripped up by geometry

texts:
ui.Chart.image seems to get tripped up by Landsat data. I've encountered the same error using Landsat 5, 7, and 8 and with several shapefiles so I'm confident it isn't a problem with the regions I'm providing. The scale I use (anywhere between 30m and 6km) doesn't make a difference either. The weird thing is this doesn't happen with MODIS data. See here and here for more on reducing MODIS data.

Both series1 and series2 plots fail with the following error:

Error generating chart: Collection.first: Error in map(ID=1):
Collection.geometry: Unable to perform this geometry operation. Please specify a non-zero error margin.

Series1 and 2 should show (for series1) individual band values over the course of the year, and (for series2) Band 1 (blue) values for individual counties over the course of the year.

Anyone have ideas on how to get past the error? Reproducible code below.

var cali = ee.FeatureCollection('ft:1S4EB6319wWW2sWQDPhDvmSBIVrD3iEmCLYB7nMM')
  .filter(ee.Filter.eq('StateName', 'California'));
print("cali shape",cali);

// This function clips images to the ROI feature collection
var clipToCol = function(image){
  return image.clipToCollection(cali);
};
print("clipToCol function",clipToCol);

// Load a Landsat collection, map the NDVI function over it.
var collection = ee.ImageCollection('LANDSAT/LT05/C01/T1_TOA')
  .filterDate('2001-01-01', '2001-12-31')
  .filterBounds(cali)
  .map(clipToCol);
print("TOA Col Preview, No Clouds",collection.first());

// Define a chart with one series in the region, averaged by DOY.
// Guide here: https://developers.google.com/earth-engine/charts_image_doy_series
var series1 = ui.Chart.image.doySeries(
    collection);
print("DOY series",series1);

var series2 = ui.Chart.image.doySeriesByRegion(
    collection, "B1", cali, ee.Reducer.mean(), 6000,
    ee.Reducer.mean(),
    'label');
print("County series", series2);


texts: In QGIS print composer, pdf in CMY/RGB, dosen't show correct colore after print

texts: By color setting in QGIS, I can only set RGB color but not CMY, so I've taken the cmy color from FreeHand, converted to RGB in Adobe Photoshop. Then exported the map as PDF (in QGIS). Before printing this PDF, I've controlled the color using of Adobe Photoshop, it shows me wrong cmy. Is this the problem with a printer or it is possible to get correct cmy or rgb color before print?


texts: Measuring areas of polygons in shapefile in the terrain

texts: How can I measure the actual areas of polygons in a shapefile in the terrain. Some of the polygons in the file is in flat terrain and some in steeper places. Is there an easy way to do this?


texts: Displaying Geodatabase Published Map On ArcGIS Server (JavaScript API)

texts:
I Google this a lot but unfortunately I didn't find any useful Tutorial for this.Can you please let me know how I can Put my published Geodatabase/ shapefile on the top of a ArcGIS JavaScript API map?

I have already published my geodatabase map as a Service now I do not know how to display that map on the top of a javaScript powered map in a page like this:

<!DOCTYPE html>

    <meta http-equiv="Content-Type" content="text/html; charset=utf-8">
    <meta name="viewport" content="initial-scale=1, maximum-scale=1,user-scalable=no"/>
    <title>Simple Map</title>
    <link rel="stylesheet" href="http://js.arcgis.com/3.8/js/esri/css/esri.css">


texts: QGIS 2.4 Processing (Zonal Statistics) error

texts:
In QGis 2.2 all my models (processing model builder) are working fine. When I start the same models in 2.4 I get errors like this:

  Prepare algorithm 5: Zonal Statistics
  Running Zonal Statistics [6/8]
  Parameters: RASTERLAYER_RASTER =C:\Temp\processing\779310258f5d4cc5a0d3f0c747b2627f\RESULT.tif, VECTORLAYER_GRENZEN =C:/Tempordner_GIS/ASL_EK/Bereinigt/018-01/DGM5_UTM32_6stell.shp
  Prepare algorithm 0: Zonal Statistics
  Running Zonal Statistics [1/1]
  Parameters: INPUT_RASTER =C:\Temp\processing\779310258f5d4cc5a0d3f0c747b2627f\RESULT.tif, RASTER_BAND =1, INPUT_VECTOR =C:/Tempordner_GIS/ASL_EK/Bereinigt/018-01/DGM5_UTM32_6stell.shp, COLUMN_PREFIX =_, GLOBAL_EXTENT =True
  Failed
  Failed
  Error executing algorithm 5 Error executing algorithm 0 'NoneType' object has no attribute 'GetGeoTransform' See log for more details

This is the error log:

  Uncaught error while executing algorithm
  Traceback (most recent call last):
  Traceback (most recent call last):
    File "C:/PROGRA~1/QGISCH~1/apps/qgis/./python/plugins\processing\core\GeoAlgorithm.py", line 217, in execute
      self.processAlgorithm(progress)
    File "C:/PROGRA~1/QGISCH~1/apps/qgis/./python/plugins\processing\algs\qgis\ZonalStatistics.py", line 78, in processAlgorithm
      geoTransform = rasterDS.GetGeoTransform()
  AttributeError: 'NoneType' object has no attribute 'GetGeoTransform'

Has anyone a tip for me?


texts: Why do I get incomplete answers using google maps api in java?

texts:
I'm using google maps API in java to check wether some places are correct, but I don't get all possible responses:

GeocodingApiRequest req = GeocodingApi.newRequest(context).address("Sydney");
try {
    GeocodingResult[] result = req.await();
    System.out.println("number of results: "+result.length);
    for(int i = 0;i < result.length;i++){
        System.out.println(result[i].formattedAddress);
        AddressComponent[] ac = result[i].addressComponents;
        for(AddressComponent acSingle : ac){
            System.out.println(acSingle.longName + " " + acSingle.shortName);
        }
    }
    // Handle successful request.
} catch (Exception e) {
    // Handle error
}

I get the following output:

number of results: 2
Sydney NSW, Australia
Sydney Sydney
New South Wales NSW
Australia AU
Sydney NSW, Australia
Sydney Sydney
New South Wales NSW
Australia AU

While, if you see what you can get by inserting the following url request, the results are many more:

Sidney

In the link above, you can see that the results are many more than 2, indeed there is Sidney also in Ohio, for instance.


texts: How to position regularly-spaced arrows along line with Geoserver?

texts:
I'm using Geoserver 2.6.2. I would like the arrows I've positionned on linear features to be regularly spaced. The fact is that what appears as a line is in fact composed of several lines of different length.

When I position an arrow mark on the centroid of each line feature, the marks are not regularly spaced and appearence is not good (see below). Is there a way to achieve that?


texts:
I am new to MapServer and I wanted to host offline map server on Windows, so after following this link https://github.com/mapserver/mapserver/wiki/RenderingOsmDataWindows#osm-import
I have setup the map for only India.

and with the below URL I can get a map of India, but I can't zoom in. How can I get to a particular lat long at center with some zoom value using URL parameters?

http://localhost/cgi-bin/mapserv.exe?mode=browse&template=openlayers&layers=all&map=C:/ms4w/apps/osm/basemaps/osm-google.map


texts: How to compute a map-to-screen transformation?

texts:
Every projection calculation I have come across dealing with the Lambert Conic projection
converts the given Lat & Long to meters.
So now the magical question: how on earth do you take the meters and convert to X, Y in pixels?
I have a map of the US. The following is the projection info:

ullat=53.3638
ullon=-130.964
lrlat=18.2977
lrlon=-69.1058

sp1=0
sp2=39
rlon=-95.5
rlat=36.5
width=640
height=480

I converted this http://svn.onelab.eu/deprecated.PLCWWW/imports/import-001-head-04-24/plot-latlong/plot-latlong over to php, but get returned meters.
I am working on an open source project.  Most of the users will not be able to install Proj4, so that is not an option (even there the return is in meters).

Anybody?


texts: All layers missing in QGIS project - not 'Bad Layers' but simply gone

texts:
I have a QGIS project I have been using for several years for a large archaeological project. Today, on opening, all layers are missing - though the groups they were in are all present and the map composer is empty.

I do not get a Bad Layers prompt, they simply aren't present in the groups. Re-adding them works but I have several hundred with a whole host of layer styles, and dozens of maps, that I would need to work through.


texts: Arcgis/Arcgis Online Domain coding table

texts:
I'm having an issue with coding domains for ArcGIS online and am hoping someone has an idea. Maybe this can be done or not.

I have to edit a .GDB to add a couple things into a feature service for some inspections. Easy enough. The issue I'm having is the new fields in the feature class properties goes straight to the bottom below globalID, Creation Date, ETC... I can move this around in the attribute table but when I'm using the collector app or am on AGOL these new fields appear at the very end. Is there anyway to adjust this?

edit: the issue is in the .GDB when I try to adjust the "fields" in the feature class properties.


texts: Error 000725 message from output workspace when using Python (ArcPy) script tool?

texts:
I am having trouble running a script tool that I just created in Arc 10.1.  The script runs fine in PythonWin, but when I create a script tool with it, the tool seems to not want to accept folders or geodatabases as its output workspace.  I have the output parameter (argument) set as both an output and a workspace in the script tool, but it still keeps giving me an error that the "Dataset.....already exists"

I'm a new user, so I can't post a picture, but the error is 000725

Here is my script, if that helps out.

import arcpy
#allow for overwrites
arcpy.env.overwriteOutput = True

#set the workspace
inWorkspace = arcpy.GetParameterAsText (0)

#set the erase feature
eraseFeature = arcpy.GetParameterAsText (1)

#set the output workspace
outWorkspace = arcpy.GetParameterAsText (2)

#get a list of all the features in the workspace
arcpy.env.workspace = inWorkspace
featureClassList = arcpy.ListFeatureClasses()

try:

    #loop through all of the features in the workspace
    for featureClass in featureClassList:

        #construct the output path
        outEraseFeature = outWorkspace + "\\erase_" + featureClass

        #perform erase
        arcpy.Erase_analysis(featureClass, eraseFeature, outEraseFeature)

        arcpy.AddMessage("Wrote clipped file " + outEraseFeature + ". ")
        print "Wrote clipped file " + outEraseFeature + ". "

except:

    # Report if there was an error
    arcpy.AddError("Could not erase feature classes")
    print "Could not erase feature classes"
    print arcpy.GetMessages()


texts: Update attributes during edit session using ArcObjects

texts:
I have made an attempt to create an ArcObjects addin to help with some workflows here in the office.  I apologize in advance if I have something fundamentally wrong with the code.  This is my first ArcObjects endeavor. Essentially I have geodatabases with a large number of feature classes.  Each feature class has the fields "DateModified" and "ModifiedBy" and I am hoping to populate them as edits are made.  I am aware of editor tracking currently in ArcGIS, but the powers that be would like to track edits without changing current data models.

It appears that the code below crashes ArcMap when it attempts to store the feature.

So my questions are 1) Can someone tell me why my code causes ArcMap to crash? 2) Are there easier ways to accomplish what I have been tasked with?  Thanks in advance.

I have the following code:

Public Class TrackEditsExtension
    Inherits ESRI.ArcGIS.Desktop.AddIns.Extension
    Public Sub New()
    End Sub
    'Invoked when the Editor Extension is loaded
    Protected Overrides Sub OnStartup()
        AddHandler Events.OnStartEditing, AddressOf Events_OnStartEditing
        AddHandler Events.OnStopEditing, AddressOf Events_OnStopEditing
    End Sub
    'Invoked at the start of the Editor Session
    Private Sub Events_OnStartEditing()
        AddHandler Events.OnCreateFeature, AddressOf Events_OnCreateChangeFeature
        AddHandler Events.OnChangeFeature, AddressOf Events_OnCreateChangeFeature
        System.Windows.Forms.MessageBox.Show("Editor tracking extension enabled.")
    End Sub
    'Invoked at the end of the Edit session
    Private Sub Events_OnStopEditing()
        RemoveHandler Events.OnCreateFeature, AddressOf Events_OnCreateChangeFeature
        RemoveHandler Events.OnChangeFeature, AddressOf Events_OnCreateChangeFeature
    End Sub
    'Invoked when a feature is created or modified
    Private Sub Events_OnCreateChangeFeature(ByVal obj As ESRI.ArcGIS.Geodatabase.IObject)
        Dim feature As IFeature = CType(obj, IFeature)
        Dim dateIndex As Integer = feature.Fields.FindField("DateModified")
        Dim nameIndex As Integer = feature.Fields.FindField("ModifiedBy")
        Dim changeMade As Boolean = False
        If dateIndex > 0 Then
            feature.Value(dateIndex) = DateTime.Now
            changeMade = True
        End If
        If nameIndex > 0 Then
            feature.Value(nameIndex) = Environ("USERNAME")
            changeMade = True
        End If
        If changeMade = True Then
            feature.Store()
        End If

    End Sub
    Protected Overrides Sub OnShutdown()
    End Sub
End Class


texts: Converting coordinate results of arcpy.AddXY to Degrees Minutes Seconds Format?

texts:
I ran the AddXY Coordinates tool on a feature class that did not have any coordinates fields. I now have two fields "POINT_X" and "POINT_Y". My objective is to convert the value in those two fields into degrees/minutes/seconds. I have tried the Convert Coordinate Notation tool but I always get empty fields in for the coordinates. Here is an example of what I am running:

arcpy.ConvertCoordinateNotation_management ("Orchids_Layer", "Orchids.shp", "POINT_X", "POINT_Y", "DD_2", "DMS_2","","",sr)

The image shows you the results I get from trying to convert "POINT_X" and "POINT_Y" using the above code. My input file's spatial reference is 'GDA_1994_MGA_Zone_56'.

What am I doing wrong?


texts: Highlight, count and delete overlapping point locations from different CSV layers in QGIS?

texts:
I have many CSV file layers in my project in QGIS, and quite a few layers have points which share the same location. Is there a way I can compare, count and highlight those "duplicates" and then delete directly by selection. I could potentially do a workaround by first creating polygon layers first via buffering each point and then intersect the point layers.

I just wonder, as I remember there is a quite easy option in ArcGIS's spatial analyst where you can do so, I think.

Any suggestions?


texts: Where is the OpenLayers plugin in QGIS 2.6?

texts:
I have just installed QGIS 2.6 and have installed the OpenLayers plugin. It lists the OpenLayers plugin in the installed tab, however, I can't add a layer to my project as it does not give the option to add an OpenLayers layer as previous versions do.

Does anyone know what is going wrong here?


texts: Passing multiple tab files into QuickImport via ArcPy only outputs file geodatabase with one feature class?

texts:
I'm trying to convert tab files to shape files. Issue is my code below converts only the last tab file in the os.walk list and not all the others. No python errors are reported and output prints the list of tab file paths below. Can anybody see something wrong as to why Quickimport converts one tab file and not the others?

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

import arcpy
import os

path = "C:/temp1"
for root, dirs, files in os.walk(path,topdown=False):
    for file in files:
        if file.endswith(".TAB"):
            table = os.path.join(root, file)
            print(table)
            output = "C:/temp1/results.gdb"
            arcpy.QuickImport_interop(table, output)


texts: saga 2.0.8 and gqis 2.0.1 on ubuntu 13.04

texts:
is it possible to have saga 2.0.8 and qgis 2.0.1 installed  on ubuntu 13.0 4at the same time?
I have the ubuntugis-repositoy activated and I got no problem with qgis and grass, but it seems impossible to install saga.

The problem is the libgdal1 library: during the qgis installation, libgdal1 h is installed too. But saga want the libgdal1 library and then the conflict.

Some idea or suggestions?

BTW: the same error occurs with the johanvdw repository (for the 2.1 version of saga).


texts: Show markers only at certain zoom levels

texts:
I'm building a cycle touring map in leaflet.js, with subjective information overlaid as markers, along the lines of "this is a great road", "this town is quite nice" etc. To avoid information overload, I want to only show the most prominent markers at a given zoom level (eg, show awesome roads highlighted at zoom 9, good roads at 11, and all road information, including negative comments, at 13+)

I have two questions here:

Is there a recognised idiom or design pattern for this kind of thing? I'm quite happy to manually decide the prominence of each marker, so I don't need any kind of clustering algorithm. There'll be maybe a couple of hundred for the whole state, so it's manageable.
How can I implement this in Leaflet? My obvious thought is to attach a handler to map.zoomend() and simply show or hide the relevant markers. Is there a better way?


texts: How to create WKT file readable by QGIS 2.4.0

texts:
I am trying to create a well-known text (wkt) file that can be imported into QGIS ver. 2.4.0

Here's a minimal working example of the geometry I wish to represent.

POLYGON((-125 52, -125 21, -98 21, -98 52, -125 52),(-120 50, -120 35, -100 35, -100 50, -120 50))

I have created a text file named polygon.wkt containing only the above wkt representation of a polygon. When I try to import it into QGIS (Layer/Add Delimited Text Layer) the import dialog responds "no data found in file."

How do I need to modify the file (a header, perhaps) to make it readable by QGIS?

Solution:
I needed to think of the file as delimited text with a header line. For some reason, I needed more than a single field to make it work. See the code block below for the working solution to the minimal example I provided.

name, geom
"my polygon", "POLYGON((-125 52, -125 21, -98 21, -98 52, -125 52),(-120 50, -120 35, -100 35, -100 50, -120 50))"


texts: How to write a new layer to PostgreSQL from QGIS plugin

texts:
The PyQGIS cookbook states that QgsVectorFileWriter can write vector layers to any format supported by OGR:

https://www.qgis.org/en/docs/pyqgis_developer_cookbook/vector.html#writing-vector-layers

However when I try this I get an error (based on the .source() of an existing layer):

error = QgsVectorFileWriter.writeAsVectorFormat(layer, u'dbname=\'snorf\' host=localhost port=5432 user=\'snorf\' sslmode=disable key=\'id\' srid=4326 type=POINT table="public"."mynewlayer" (the_geom) sql=', 'UTF-8', layer.crs(), 'PostgreSQL')

This returns QgsVectorFileWriter.ErrCreateDataSource.

I don't have any problems modifying existing PostgreSQL layers, and I can write a new layer to the database using ogr2ogr without issue.

How do I create a new empty postgres layer?
How do I save an existing layer into a postgres database?


texts: QGIS 2.18.3 path and sip module cannot be load in ubuntu 16.04

texts:
I'm currently installing QGIS 2.18.3 it works fine at first, after upgrading my ubuntu the problems comes out every time I run my QGIS; it always shows "Couldn't load SIP module. Python support will be disabled"

I figure out that this is happening because my sip python module installation is not in QGIS path. Then invoke my QGIS with path to my sip module where it is installed:

$ qgis --configpath /usr/lib/python2.7/dist-packages/

And now the python feature in QGIS works normally.
My question is how can I make this path permanent every time I run QGIS?


texts:
I have LineString in openlayers-3.  I want to set 2 label for this feature In the start and end of the feature. An example is as follow:

How Can I do this labels?


texts: Can I add Mapbox data (markers) to CartoDB basemap

texts: I have prepared a basemap in Mapbox which includes imported data (basically, marker points). I was expecting these markers to show when using that project as a basemap in CartoDB, but it only shows the basic tiles. Is anyone aware of this issue or how I could go about including the markers from Mapbox in CartoDB?


texts: Getting ERROR 000260 item with same path name already while creating MPK data packet

texts:
I am using ArcMap 10.7.
I have an MXD which sources data from a single ArcGIS Enterprise MSSQL Database. Every layer is referenced only once and there are no other mxd documents open. The MXD has 20 Layers and when I try to create an MPK file using ArcToolbox > "create map Package" Tool I always get this error. I need to set the "Engine Runtime" parameter which is why I use the tool from the toolbox. I have been able to successfully create 8 other packages using other data from the same source.
As a workaround I copied the layers into a file geodatabase and was able to create an MPK but this is not a solution, because the data is being updated every week, and this is only one of 11 mxds which have to be copied.
Does anyone have a solution?


texts: Does OpenLayers mergenewparams repeatedly hit a WMS server when you change the layer name?

texts:
Does mergenewparams repeatedly hit a WMS server if you change the layer name?

For example, you call the server with layer 'a' then you use mergewnewparams and call layer 'b', then you use mergewnewparams and call layer 'a' again. Does it cache layer 'a' or request a new layer 'a' ?

I'm actually using a tilecache, but I assume it works the same for either.


texts: Cumulative distance along a raster stream network

texts:
I am looking to calculate the cells along a stream network from starting points to a maximum distance and stopping points. I need something like r.cost in GRASS (https://grass.osgeo.org/grass75/manuals/r.cost.html), the only difference is that I need an algorithm capable to follow the stream network topology (without distinguishing between upstream and downstream). The problem is that r.cost only consider contiguous cells, but with streams in raster format, it is possible that contiguous cells are not really connected due to resolution issues, specially at river junctions.

In the picture, there is an example of issues with r.cost. The red point is a starting point, black and gray cells correspond to two different rivers originally connected only at the junction, but spatially contiguous in more points in raster format. Any solution in GRASS or SAGA?


texts:
I have code that suppose to calculte NDVI and then to cound the number of pixels in the raster that are higher than 0.3 There is part of the code where there is function that suppose to select the band NDVI and then calculte number of pixels. for some reason now, I have gotten the following error:

  Image.select: Pattern 'NDVI_mean' did not match any bands.

I have checked many times and this is the name of the band so I don't know why it says it doesn't match.
Also, when I check the images in the band they do have this name:

this is the fucntion:

var listOfImages =(byYear.toList(byYear.size()));
print(listOfImages, 'list of images');

var stats = listOfImages.map(function(img){
  var image = ee.Image(img).select('NDVI_mean');

  var c03 = image.gt(0.3).selfMask().reduceRegion({
    reducer: ee.Reducer.count(),
    geometry: geometry,
    crs: 'EPSG:4326',
    //maxPixels:1310361348,
    scale: 30,
  }).values().get(0);

  var c05 = image.gt(0.5).selfMask().reduceRegion({
    reducer: ee.Reducer.count(),
    geometry: geometry,
    crs: 'EPSG:4326',
    //maxPixels:1310361348,
    scale: 30,
  }).values().get(0);

    var c07 = image.gt(0.7).selfMask().reduceRegion({
    reducer: ee.Reducer.count(),
    geometry: geometry,
    crs: 'EPSG:4326',
    //maxPixels:1310361348,
    scale: 30,
  }).values().get(0);

  return image.set('c03', c03).set('c05', c05).set('c07', c07);
});

print(stats);

This is the full code, it's a bit long , I have gone over all of it but couldn't find the problem with the band name:
https://code.earthengine.google.com/a6726d65d9cddb10aa23f7683464a386

the polygon was the city abidjan

My end goal is to be able to calculte the number of pixels in each image that are higher than 0.3


texts: Adjusting data digitized from erroneously georeferenced basemap?

texts:
I have made quite a big georeferencing error, by attributing one of the control points of the map with a value that is 2 km less than the actual value, on the horizontal axis. The georeferenced map looked more like a right-angle trapezoid. I realized this only after digitizing.

Can I tweak/adjust somehow the digitized lines so they become correct according to the new shape of the map?

I have imagined rescaling around a point that is different than the center of the map, in Global Mapper, but I seem to be long way off achieving anything.


texts: Band rerendering netcdf file by variable not time. Using gdal, nco, cdo

texts: I want to average netcdf data (band1: precipitation) in bins (i.e. 0.01,0.10,0.25,0.5,0.75,1.0, etc.). What nco, gdal, or nco tools can I use to regrid the netcdf file this way?


texts:
I'm trying to convert coordinates from WGS84 to EPSG:3067 with proj4js and Openlayers, but can't get the conversion right. When I run proj4('EPSG:3067',[-24.9993393, -60.184716]); it gives
[-2143537.623518457, -7831430.649243464]
When clicked on the map the point of coordinates it gives
[-2182546.021603215, -12759115.371540245] which is accurate if I put the center of the map to the clicked coordinates. Below is the code where the EPSG:3067 is defined.

proj4.defs("EPSG:3067", "+proj=utm +zone=35 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs");
extent = [-3670733.46, 4601971.85, 642319.78, 9362767.00];
ol.proj.get('EPSG:3067').setExtent(extent);


texts:
I have a series of circles and lines representing some information about a geographical location. Is it possible to reuse that (complex) geometry by simply attaching it to a single point?

In other words, I would like to use that set (of geometries) as a prototype geometry for other locations, is that possible?


texts: Using a bounding polygon, instead of a bounding box

texts:
I've got a large set of LAT/LON coordinates. Up until now, when I needed to select out those that fell within a certain region, I'd just create a bounding box using some MINx, MAXx, MINy, MAXy values and try to get it as close as possible to the area I was interested in.

I'd like to get more precise by using more than 4 points as the extent of my search area, in this specific case, 14 points.

I'll probably be using Python or SQL for this, so the specific syntax on any solution isn't that important, I'm just interested in the formula for figuring this one out.


texts:
We have some jp2 or img satellite images that we compress to .ecw using ERDAS Imagine. We have realized that with some images we are losing too much information while doing that.
What i wanted to ask is about preferences on ERDAS Imagine.

Target ECW Compression (Greyscale): Value

Target ECW Compression (RGB):  Value

Target ECW Compression (Multispectral):  Value

Do I need to adjust them each time I compress an image or is there an optimal setting that i can use?

(We are using mostly 16 bit satellite images as input)

Or where else could the problem be that the ecw files are so foggy and unsharp compared to the input files.

img:

ecw:

Thanks in advance,

Yasmin


texts: How to show Vector Layer created by v.clean GRASS tool in QGIS

texts: I have created vector layer using GRASS tool v.snap, which resulted in a vector layer as output. I want to display the vector layer in QGIS. How can i view the vector layer created by snap tool in QGIS that is created by GRASS tool. After search on internet i found that there is a button used to Add GRASS vector layer in QGIS plugin. But in my current version that is 2.18. I did not find a option for this.


texts:
I tried to save my project under a new name ("Save as"), as I often do. However, the screen just sort of blinks; then stays the same and does not allow me to enter the new name, etc.. I have upgraded to the latest version of QGIS. I tried doing this with an earlier version of my project (saved under a different name) and it worked ok. Any idea what's the matter with this file, or what I need to do?

Here is a link to it.

Here's a link to the previous version, which did not have this problem.


texts: How to use ST_VoronoiPolygons() in PostGIS?

texts:
Sorry folks, I know it sounds like a lame question - and, yes I did read the manual... :-)

I'm trying to make Voronoi diagrams around points I have in a PostgreSQL table.

It should be rather easy - but I can't find a single example on using ST_VoronoiPolygons() on real data (i.e. a table of points).

All I keep seeing is examples similar to that in the docs. like this:

SELECT
    ST_VoronoiPolygons(geom) AS geom
    FROM (SELECT 'MULTIPOINT (50 30, 60 30, 100 100,10 150, 110 120)'::geometry AS geom ) AS g
;

That works fine on my computer. But I want to make Voronoi polygons around my real point, not just example-points hard-coded in the SQL.

Can anybody point me to a working example of SQL that shows the use of ST_VoronoiPolygons() with an existing PostgreSQL Point table as input.


texts: Transform from Lambert Azimuthal "old" to ETRS89_LAEA

texts:
I'm currently struggling to convert a small .asc dataset with no defined projection.
The projection that has to be first applied is Lambert Azimuthal "old" projection:

In PROJ.4

+proj=laea +lat_0=48 +lon_0=9 +x_0=0 +y_0=0 +R=6378388 +units=m

or
In ArcGIS

Radius of sphere of reference  6378388m
Longitude of center of projection 9 0 0 (Degrees Minutes Seconds)
Latitude of center of projection 48 0 0 (Degrees Minutes Seconds)
False Easting 0 m
False Northing 0 m

projection lambert_azimuthal
units meters
spheroid sphere

The target projection is EPSG:3035 with

PROJ.4

+proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs.

or in ArcGIS

Projected Coordinate System:    ETRS_1989_LAEA
Projection: Lambert_Azimuthal_Equal_Area
False_Easting:  4321000.00000000
False_Northing: 3210000.00000000
Central_Meridian:   10.00000000
Latitude_Of_Origin: 52.00000000
Linear Unit:    Meter

Geographic Coordinate System:   GCS_ETRS_1989
Datum:  D_ETRS_1989
Prime Meridian:     Greenwich
Angular Unit:   Degree

I'm using the Warp (QGis tool),  and this is what I got.

Output driver **`AAIGrid' not recognised or does not support
direct output file creation. The following format drivers are configured
and support direct output**:
GeoRaster: Oracle Spatial GeoRaster
VRT: Virtual Raster
GTiff: GeoTIFF
NITF: National Imagery Transmission Format
HFA: Erdas Imagine Images (.img)
ELAS: ELAS
MEM: In Memory Raster
BMP: MS Windows Device Independent Bitmap
PCIDSK: PCIDSK Database File
ILWIS: ILWIS Raster Map
SGI: SGI Image File Format 1.0
Leveller: Leveller heightfield
Terragen: Terragen heightfield
netCDF: Network Common Data Format
HDF4Image: HDF4 Dataset
ISIS2: USGS Astrogeology ISIS cube (Version 2)
ERS: ERMapper .ers Labelled
RMF: Raster Matrix Format
RST: Idrisi Raster A.1
INGR: Intergraph Raster
GSBG: Golden Software Binary Grid (.grd)
GS7BG: Golden Software 7 Binary Grid (.grd)
PNM: Portable Pixmap Format (netpbm)
ENVI: ENVI .hdr Labelled
EHdr: ESRI .hdr Labelled
PAux: PCI .aux Labelled
MFF: Vexcel MFF Raster
MFF2: Vexcel MFF2 (HKV) Raster
BT: VTP .bt (Binary Terrain) 1.3 Format
LAN: Erdas .LAN/.GIS
IDA: Image Data and Analysis
GTX: NOAA Vertical Datum .GTX
NTv2: NTv2 Datum Grid Shift
CTable2: CTable2 Datum Grid Shift
KRO: KOLOR Raw
ADRG: ARC Digitized Raster Graphics
SAGA: SAGA GIS Binary Grid (.sdat)

What I really don't know is why I cannot find or define a new projection for the .asc files and then hopefully to transform them into the ETRS_1989_LAEA.
So these probably has to be done like this (?)

Define the default projection for the dataset
Transform the projection to the ETRS_1989_LAEA

My question mainly focuses on (1).Any idea about how it can be done in ArcMap or even in QGIS?


texts: QGIS layer: import CSV and set 'Display HTML Map Tip' using python

texts:
My objective is to generate a python plug-in for QGIS that will open a CSV and define a field within the table to display as a map tip.

I'm trying to achieve this task and have got stuck on trying to set the layer properties using python, rather than opening the 'layer properties' dialog, selecting the Display tab and entering HTML expression...

the HTML expression I'm trying to insert into the layer properties is:

field1= '\"title\"
\n'+'[%CONCAT(\'\')%]'

where "title" and "file" are fields within the CSV, giving rise to:

"title"

[%CONCAT('')%]

This gives me the right HTML to slot into a relevant place - and it works when I use this manually, but I'm now lost trying to find how to change the layer properties using python...

There seems to be a function called setDisplayAttributes but I can't really follow the API description well enough to understand whether it is the correct function and I can't find any examples using this function.

I'm a toddler with Python but probably pre-pubescent with QGIS (Version 2.8 on Windows 10).


texts: Automating removal of all demo layers from GeoServer?

texts:
Is there an easy way to empty out an install of a standalone Geoserver jar (i.e., delete all the examples) or do I have to do it through either the admin GUI or the RESTful interface?

Can I just delete or munge some files in the work directory?

Or is there a command line switch to tell it to never unpack the examples in the first place?


texts: QGIS 2.16.3 crashes when trying to dissolve any shapefile

texts: When I try and use the dissolve tool (vector>geoprocessing tools>dissolve) QGIS 2.16.3 crashes and produces a minidump file. This is with any shapefile that I use. I have also tried it on a friends 2.14.3 and the dissolve process works fine for the same files that I tried.


texts: How do I clip/mask a satellite photo to a polygon?

texts:
Given a shapefile composed of contiguous polygons, I'd like to clip (or mask) regions of satellite photo to the polygon outlines, and export each resulting raster in medium resolution for further processing.

I've gotten as far as overlaying the shapefile on Google Maps imagery (in QGIS, with Openlayers plugin), but I can't find any operation to intersect vector data with raster data (presuming that the Openlayers layer is indeed treated as a raster layer).

A simple solution would be great, and suggestions of other free OSX-compatible software are fine.  A programmatic solution would be ideal (though not requisite) as I'd like to export many such images, each named by an attribute field in the shapefile.

I hope I haven't over-explained a simple problem, but alas, I'm rather new to GIS software.  Trying to move forward with this project as quickly as possible, but stymied by data acquisition... Thanks a ton!


texts: How to write expression with vb.net in mapwingis (ShapefileCategory.Expression)?

texts:
I would like to give color to the map with new category. to make must make previous expression. I've tried but still failed. What is wrong ? Please help me.

 Dim category As MapWinGIS.ShapefileCategory = sf.Categories.Add("Manual")
Dim ad As String = "Alaska"
Dim UniqueShape As New MapWinGIS.Shapefile
UniqueShape = FormMain.AxMapMain.get_GetObject(idx)

Dim af As String = UniqueShape.Table.Field(1).Name
Const cote As String = """"

'failed
' category.Expression = "[STATE_NAME] = """ + ad + """"
'failed
' category.Expression = "[STATE_NAME] <> """""
'failed
category.Expression = "[" + af + "] = " + ad
'failed
category.Expression = "[STATE_NAME] = " & cote & "Alaska" & cote
'failed
'category.Expression = "[STATE_NAME] = ""Alaska"""

Dim utils As New MapWinGIS.Utils
category.DrawingOptions.FillColor = utils.ColorByName(MapWinGIS.tkMapColor.Magenta)
sf.Categories.ApplyExpression(0)

sf.DefaultDrawingOptions.Visible = False
' FormMain.AxMapMain.AddLayer(sf, True)
FormMain.AxMapMain.Redraw()
FormMain.AxMapMain.Refresh()


texts:
I was under the impression that L.circle changes size with zoom, but L.circleMarker should remain the same radius as you initially set?

I've setup an example here, and whichever I choose the circle always changes size on zoom. What am I missing?

http://jsfiddle.net/qmL8ej19/

(double click to add circle marker)


texts: Is there a Cluster convex hull implementation in OpenLayers?

texts:
In this leaflet example. When hovered over a cluster marker, it shows the bounding area of the clustered features which is known as convex hull.

Is there a way to reproduce it with OpenLayers?


texts: Command line-tool to convert line-delimited geojson file to regular geojson

texts:
Now I need the opposite tool from this one: something to convert a line-delimited GeoJSON file to a single GeoJSON file.

In other words:

Write { "type": "FeatureCollection", "features": [ on the first line.
Append all the lines from the input file, with a comma at the end of all of them except the last.
Then write ]}

Is there such a tool?


texts: What Proj.4 parameters have to be modified if the Y coordinates are shifted?

texts:
I would like to transform Barcelona map coordinates to Latitude, Longitude with Proj.4. I'm using EPSG number 23031 and the parameters are:
+proj=utm +zone=31 +ellps=intl +towgs84=-87,-98,-121,0,0,0,0 +units=m +no_defs

However, the Y coordinates in the original maps are shifted by -4000000.

For example:
x=438120.5406, y=568084.9572 instead of x=438120.5406, y=4568084.9572.

I realized that I have to modify the +towgs84 values, but I can't find the exact values.

What is the correct method?


texts: ArcGIS Trace Geometric Network Tool Trace using Find Common Ancestors gives error

texts:
I have a geometric network stored in an ArcGIS 10.5 geodatabase. These lines represent all the river paths in a watershed.  The lines associated with this network have the required flow direction assigned based on line digitized directions.  I can manually select an end node of one stream and use the trace geometric network tool in ArcToolbox to select a single upstream path without error. My network branches and I want to select an end node and return all the connected streams that flow to that end node. If I change the parameter to Find Common Ancestors the tool returns an error 001191: No flag found from the flag junction feature class.  Ultimately I need to generate a Python script to select each nodes and select all upstream line segments but I cannot make the tool work.  The  inputted river lines, and the end nodes from those lines each share a CRS and reside in the same feature database feature class items.

Here is a screenshot of a trace successful find uphill path result using these data. I want to select all connected lines though.

If I use the Utilities Network Toolbar I can place a flag at an end and trace a single path like the image above.  The same tools that fail using the Trace Geometric Network tool in toolbox are also grayed out in the Utilities Network Toolbar.

Is there a bug in the Find Common Ancestors parameter of the Trace Geometric Network tool?


texts: Automated removal of dangles / loose ends in line with QGIS?

texts:
Here's what I want to do:

The original line is a dissolved osm extract of ways, which I want to simplify to one continous line without dangles / loose ends (red x's in picture).

I found the network generalization tool in GRASS, but possibly there is a QGIS way ro remove loose ends.

If there is no QGIS-only approach to do it, what would the parameterization in GRASS would have to look like to yield the desired result?