texts: Combining static text and dynamic text into single field?

texts:
I have a print template that is used with the Export Web Map task. This print  template contains place-holder fields that are replaced with actual values when the web map JSON is sent to the task for PDF generation. I have the place-holder fields being replaced with the real values but I need to add some static text to identify what the value represents. For example, instead of just throwing "123" on the top of the map it would be better to display "Lift Station ST-123".

The place-holder fields appear to be plain old text elements (Insert > Text) with an Element Name. I copied the place-holders from one map and pasted in the new print template so I might be missing something about how they SHOULD be created. I used Insert > Text to place a piece of text and, when compared to the Properties dialog of the place-holder, they look to be the same.

How do I add static text to this place-holder?

Along with adding some context to the displayed value, a single element with static and dynamic text is easier to move than two independent entities.

Update: I have the following text elements on my map.

text property = dp id : <dyn type="document" property="XXX"/>. On the Size and Position tab, the Element Name is blank. When printed, the static text remains but the XXX property is not inserted.

text property = pp Well No. : <dyn type="page" property="XXX"/>. On the Size and Position tab, the Element Name is blank. When printed, the static text remains but the XXX property is not inserted.

text property = XXX. On the Size and Position tab, the Element Name is XXX, exactly matching the CustomTextElement in the WebMapAsJson submitted to the ExportWebMap task. When printed, the XXX text is replaced by the value submitted in the JSON. This is the expected behavior but now I want to prepend static text to the dynamic value provided in JSON.

text property = dp w/ element name: <dyn type="document" property="XXX" />. On the Size and Position tab, the Element Name is XXX, exactly matching the CustomTextElement in the WebMapAsJson submitted to the ExportWebMap task. When printed, the entire value of the text property is removed and replaced by the value submitted in the JSON. The text displays the the dynamic value provided in the JSON (good) but strips the static text (bad).


texts: Updating Values in FIelds of ArcGIS editing only portion?

texts:
I have total 96 records which are having duplicate values like 999999.

What i want is to convert these values into unique ones like 999901,999902,999903 and so on, through the 96 records.

I want only 2 digits to be replaced as rowids but there may also be cases where i have to change three digits too.


texts: I am using postGIS. I am using ST_AsGeoJSON(geom) function to get GeoJSON data of Geom field. Its returning a long lat value that have 13 digits after decimal point. Is there any way to reduce or limit the decimal digits?


texts:
I'm trying to import vector data into SAGA GIS from a SQLExpress 2008 server.
Firstly is this possible?

I would have thought using an OGR VRT layer would have been the way to do it, I've attached a sample of my vrt file, any help would be great.

<OGRVRTDataSource>
<OGRVRTLayer name="tempGeom">
    <SrcDataSource>MSSQL:server=(local)\SQLEXPRESS;database=test;trusted_connection=yes</SrcDataSource>
    <SrcSQL>select *, geometry::STGeomFromText('POINT (' + CAST(x AS VARCHAR) + ' ' + CAST(z AS VARCHAR) +' '+ CAST(y AS VARCHAR) +')',0) as geom  from dbo.testTable</SrcSQL>
</OGRVRTLayer>


texts: SQL server views in GeoServer - problem with SQL

texts:
I try to extract non-spatial data from mssql as spatial data to my GeoServer 2.1.0 installation from a SQL statement looking like this:

SELECT  database.table1.id,
    geometry::STPointFromText('POINT(' + CAST(database.table1.YCOORD AS VARCHAR(6)) + ' ' + CAST(database.table1.XCOORD AS VARCHAR(7)) + ')', 3008) AS geometry
FROM database.table1

The above given statement returns good spatial info in mssql management studio but I'm not sure how to get it working in GeoServer (I have to CAST the data as it is stored as integer). Setting up mssql as a store works fine and I know I can access its tables. I then try to execute the statement as a SQL view. That does not work :(

The geometry attribute shows up as byte[] and I can't change that.

I can't write anything in native SRS and filling in epsg:3008 in declared SRS doesn't seem to do anything when calculating the bounding box, giving the values MIN X & Y as -1 and the MAX X & Y as 0 in native bounds.

@jdeolive - not an option to try trunk at this stage, I have to convince my boss who is a die-hard ESRI-fan and stability is a must.

@unicoletti - I guess I have to read my statements better not to mix geometry and geography objects :-) thank you. It makes sense when calculating bonds, will try it tonight. Creating a view of it is what I have in mind later on. The only thing is that I'm taking my information from another view but as it is only about 6500 points I don't think it will give any larger performance issues.


texts: Erratic results from filters.normal in PDAL?

texts:
Given this roof structure:

I get the following output for the normalx, after using:

{
         "type":"filters.normal",
         "knn":8
},

(orange to purple : -1 to 1)

Comparable situation for normaly and normalz.

Where I was expecting homogeneous colored roofplanes, I get a mixed-up situation, where the values appear to be sometimes flipped to the negative.

It seems to me there are some artifacts here based on the laser scan-line, considering the flipping of the normal on the same roof plane.

I cannot wrap my head around what happens here. Ideas?

The purpose of this work is to get the slope and aspect from the roofplanes. I can get the slope correct by taking the absolute values of normalz (slope will never pass the 90 degrees) but this is not an option for calculating the aspect it seems.

The resulting point set can be downloaded here.


texts: Geoserver : add extra style (SLD) with Curl

texts:
In my geoserver I have a layer with a default style (SLD) and i want to add another style to this layer with Curl.

I tried command below with <style></style> instead of  <defaultStyle></defaultStyle> but without success :

curl -u admin:pass -XPUT -H 'Content-type: text/xml' -d '<layer><style><name>sld_name</name><workspace>workspace_name</workspace></style></layer>' http://localhost:8080/geoserver/rest/layers/workspace_name:layer_name

Unfortunately I don't have found help in geoserver cookbook : http://docs.geoserver.org/latest/en/user/rest/examples/curl.html.

Anybody know the answer ?


texts:
I encounter edge issues when trying to mosaic USGS NED DEMs in ArcGIS 10.0.0. After projection and aggregation, I attempt arcpy.MosaicToNewRaster_management(). I've reproduced the same results using .tif, .img, and GRID formats. I've also reproduced the result using file geodatabases and file folders to store inputs and outputs. I've also tried using arcpy.Mosaic_management() with no luck. Note the unwelcome black border on the bottom which is problematic for analytical applications.

What can be done about the edge issues here?

This question was split from parallelization of raster mosaic.

 +  =


texts: High coloum and line number while using the Raster Calculator in QGIS

texts:
I downloaded the DEM25 from the EU DEM homepage and want to create in the end polygones for the exposition for an area of around 5000m² (like all mountainsites facing south are yellow etc.). So far, everything worked, except that in the Raster Calculator i get a huge number of coloums(aound 50000) and lines(45000) which are resulting in a big outputfile (over 10 Gig) after i used "actual layerdimension".
I tried to reduse the coloums (to 10000) and lines(9500). That worked, but i am not sure if i have overseen some big loss of data which i could regret later.
I had other DGM for a similiar region and in this i had just 2401 coloums and 2001 lines.
Has anybody an idea why the numbers are so high?


texts: Problems using REST API to generate Tile Packages from World Topo (For Export)

texts:
Using the python script below it seems that I can request the generation of tile packages from an ESRI ArcMap REST service, and then download the resulting tpk file generated by ESRI.  But when I try and open the package NOTHING shows up in ArcMap.  Can't tell what is happening - can't tell what is going wrong.  Any ideas?

import arcpy
import requests
import pandas as pd
import time
import json

def gen_token():
    arcpy.SignInToPortal_server("USERNAME", "PASSWORD", "")
    token = arcpy.GetSigninToken()
    return str(token['token'])

def check_status(jobID, token):
    req_status_url = 'http://tiledbasemaps.arcgis.com/arcgis/rest/services/World_Topo_Map/MapServer/jobs/' \
                     + jobID \
                     + '?f=json' \
                       '&token=' + token

    status_resp = requests.get(req_status_url)
    status_resp_cont = json.loads(status_resp.content)
    return str(status_resp_cont['jobStatus'])

def scrape_tiles(input_file, levels):

    token = gen_token()

    tpk_grid_df = pd.read_csv(input_file, header=0)

    for index,row in tpk_grid_df.iterrows():

        extent = '{"xmin":' + str(row['XMIN']) + \
                 ',"ymin":'+str(row['YMIN']) + \
                 ',"xmax":'+str(row['XMAX']) + \
                 ',"ymax":'+str(row['YMAX']) + \
                 ',"spatialReference" : {"wkid" : 4326}}'

        req_url = 'http://tiledbasemaps.arcgis.com/arcgis/rest/services/World_Topo_Map/MapServer/exportTiles?' \
                  'tilePackage=true' \
                  '&exportExtent=' + extent + \
                  '&optimizeTilesForSize=true' \
                  '&compressionQuality=90' \
                  '&exportBy=levelId' \
                  '&levels=' + str(levels) + \
                  '&f=json' \
                  '&token=' + token

        resp = requests.get(req_url)
        resp_cont_dic = json.loads(resp.content)
        jobID = str(resp_cont_dic['jobId'])
        jobStatus = str(resp_cont_dic['jobStatus'])
        print(jobStatus)

        while jobStatus == 'esriJobSubmitted' or jobStatus == 'esriJobExecuting':
            time.sleep(15)
            jobStatus = check_status(jobID, token)

        if jobStatus == 'esriJobFailed':
            print('esriJobFailed')
            break

        # GRABBING THE RESULTS
        res_url = 'http://tiledbasemaps.arcgis.com/arcgis/rest/services/World_Topo_Map/MapServer/jobs/' \
                  + str(jobID) \
                  + '/results/out_service_url?f=json' \
                    '&token=' + token

        result_resp = requests.get(res_url)
        result_resp_cont_dic = json.loads(result_resp.content)
        result_url = result_resp_cont_dic['value'] + '/Layers.tpk'

        file_name = 'NH_TP_' + str(levels) + '_' +  str(row['ID']) + '.tpk'

        r = requests.get(result_url)
        with open(file_name, 'wb') as f:
            for chunk in r.iter_content(chunk_size=1024):
                if chunk:  # filter out keep-alive new chunks
                    f.write(chunk)
                    f.flush() #commented by recommendation from J.F.Sebastian
        f.close()

if __name__ == "__main__":

    input_file = "TPK_Grid_100.csv"
    num_levels = 17
    scrape_tiles(input_file, num_levels)

First couple of lines from the input file:

ID  XMIN    YMIN    XMAX    YMAX
1   -72.573919  42.633528   -72.37901   42.909766
2   -72.37901   42.633528   -72.1841    42.909766


texts: Seeking GPS device for vehicle tracking?

texts:
I am a programmer with no connection to Geographic information, except that I work on application which uses Google Maps.

I need some device, which I will be able to set in a vehicle in order to keep track of it.

Every second (it would be better if this parameter can be configurable) a device should send information to the server:

GET/POST request to http://mysite/getLocation

with the following information:

lat: 23.1232
lng: 34.1232
deviceId: 321

Few important things:

information should be sent in real time (just storing info and
accessing it in the end of the day does not work)
vehicle will be operating in the city
device should work at least 24 hours till recharging

What device should I buy?

If possible, provide some links to the device, or where I can read about it.


texts: Google Earth Engine Coordinates from Table

texts:
I am trying to import Table data into GEE. The original CSV file is at TERN Star Transects. I have followed the steps at GEE Importing Table Data (including making sure to specify either the x and y columns as "ref_x" and "ref_y", or the primary geometry column as "geom"), however I'm encountering two issues.

Ingesting the Table (or seemingly any other table from a CSV file) increases the absolute value of the x and y coordinates of any point by around 7.
Points appear as (160.84, 34.56) when they should be (160.84, -34.56), as they are in the CSV. This missing minus sign results in the points being placed off the coast of Japan when they should only be over Australia.

The code I'm using to display the points is
var STAR_TRANSECTS = ee.FeatureCollection("your_file_location");
print('STAR_TRANSECTS',STAR_TRANSECTS);
Map.addLayer(STAR_TRANSECTS, {}, 'STAR_TRANSECTS');

Some things I've tried:

A fresh new CSV test file: fixes the problem with negatives not appearing but doesn't solve the weird scaling issue.
Using someone else's CSV file they had attached to a StackExchange question. This works fine, including when I add negative values to their table!? Adding a CSV with GPS coordinates points to Google Earth Engine
Pulling the CSV file into R. R is able to read the negative sign. Exporting as a CSV from R and trying to upload that file to GEE did not solve the issues.

More details to help solve this:

Previous thread on Tables and GEE:
Adding x,y coordinates from table in Google Earth Engine
The website I am retrieving the CSV file from (includes some more
info on its properties)
http://data.auscover.org.au/xwiki/bin/view/Product+pages/SLATS+Star+Transect+Data

At this point I'm not sure if I'm completely misunderstanding something about GEE or if there is a problem with the original data.


texts: Debugging AttributeError: Invalid destination path from ArcPy?

texts:
I'm doing some programming work for a client and we can't seem to get the script tool to work on her end. It works perfectly on my end but when she tries it, it grinds away for a long time and then gives her an AttributeError: Invalid destination path. I asked what version of ArcMap she's using and she says 10.1. I'm using 10.3. Could the different versions be the cause of the problem?

import arcpy

inputList = arcpy.GetParameterAsText(0)
bckGrndLyrs = arcpy.GetParameterAsText(1)
pageRng = arcpy.GetParameterAsText(2)
outLoc = arcpy.GetParameterAsText(3)
exportDPI = arcpy.GetParameterAsText(4)

mxd = arcpy.mapping.MapDocument("CURRENT")
df = arcpy.mapping.ListDataFrames(mxd)[0]

for lyrs in arcpy.mapping.ListLayers(mxd, "", df):
    if lyrs.name in inputList:
        lyrs.visible = False
    elif lyrs.name in bckGrndLyrs:
        lyrs.visible = True
    else:
        lyrs.visible = False

for inFC in arcpy.mapping.ListLayers(inputList, "", df):
    if inFC.name in inputList:
        arcpy.AddMessage("Activating " + inFC.name + " Layer")
        inFC.visible = True
        for titleItem in arcpy.mapping.ListLayoutElements(mxd, "TEXT_ELEMENT"):
            if titleItem.name == 'EditedTitle':
                arcpy.AddMessage("...Changing title to " + inFC.name)
                titleItem.text = inFC.name
                arcpy.AddMessage("......Exporting pages " + pageRng + " as PDF")
            for pageNum in range(1, mxd.dataDrivenPages.pageCount + 1):
                if str(pageNum) in pageRng:
                    mxd.dataDrivenPages.currentPageID = pageNum
                    blockName = mxd.dataDrivenPages.pageRow.Block
                    arcpy.RefreshTOC()
                    arcpy.RefreshActiveView()
                    mxd.dataDrivenPages.exportToPDF(outLoc + "\\" + inFC.name + ".pdf", page_range_type="RANGE", page_range_string=pageRng, multiple_files="PDF_MULTIPLE_FILES_PAGE_NAME", resolution=exportDPI)
            arcpy.AddMessage(".........Deactivating " + inFC.name + " Layer")
            inFC.visible = False

del mxd

AttributeError: Invalid Destination Path


texts: Equivalent of OpenLayers 2 ol.Layer.Text in ol3 (to show POIs)

texts:
I have this OpenLayers 2 code that I want to convert to OpenLayers 3:

var pois = new ol.Layer.Text(
    "My POIs",
    { location:"/media/pois.txt",
      projection: ...
    }
);
map.addLayers([
    pois
]);

The .txt file is a tab-separated file with the following columns:

lat lon title description icon iconSize iconOffset

What is the equivalent ol.source to achieve the same effect, and do I need to change anything in the input file?


texts: Export one attribute field from QGIS to a .TXT or .CSV file

texts:
I used this code in Python to export attributes from QGIS to a .TXT file :

##input_file=vector

from qgis.core import *
import processing
# Load the layer
layer = processing.getObject(input_file)

# Set the path for the output file
output_file = open('C:/.../cities.txt', 'w')

# Get the features and properly rewrite them as lines
for feat in layer.getFeatures():
    msgout = '%s, %s\n' % (feat["name"], feat["pop_max"])
    unicode_message = msgout.encode('utf-8')
    output_file.write(unicode_message)
output_file.close()

It worked, but I don't want to get all fields of the attribute table, but just one of them (ID field for example).

Will it be possible?


texts: Qgis. Layers not shown when trying to import gdb file. 1 combined layer added instead

texts:
i am needing a little help, as somewhat a beginner, to work with .gdb files.
When trying to import the clc data for landcover in germany (https://opendata-esri-de.opendata.arcgis.com/datasets/esri-de-content::corine-landbedeckung-2018), i run into the problem, that via "add layer"-> "vector layer", while using "directory"  and "OpenfileGDB", qgis does not provide me witha list of layers to chose from in the next step, but adds just 1 layer with all the categories in it instantly, without prompting the "chose from Layer"window. Is there any pre-set i am missing, to enable to let me chose only specific layers from a gdb when importing?
Version 3.12.2
Many thanks to whoever might answer. Tim


texts:
Recently I was learning how to set up qgis server.

So here's the trouble I've encountered so far:

I was using OSGEO Live 5 Distribution ISO' Installed in dual boot mode with windows 7
I got the GetCapabilities work after fresh install.
the response was as indicate in qgis server manual, the XML info about layer etc
After I execute update software system, then I tried to do the GetCapabilities returned 500 Internal Server Error.

I didn't do any change in the qgis system which was worked fine before software update.

What should I do to tackle this error

I'm very new to linux or foss, it just my work project have positioned me to like it or not to learn linux and foss.

so I followed up advice I open the error.log have information such as following:

[Sat Dec 24 22:33:54 2011] [error] mod_fcgid: process /usr/lib/cgi-bin/qgis_mapserv.fcgi(1862) exit(communication error), get unexpected signal 6
[Sat Dec 24 22:33:59 2011] [warn] [client 127.0.0.1] mod_fcgid: error reading data, FastCGI server closed connection
[Sat Dec 24 22:33:59 2011] [error] [client 127.0.0.1] Premature end of script headers: qgis_mapserv.fcgi
[Sat Dec 24 22:34:05 2011] [error] mod_fcgid: process /usr/lib/cgi-bin/bik/qgis_mapserv.fcgi(1864) exit(communication error), get unexpected signal 6
[Sat Dec 24 22:35:10 2011] [notice] caught SIGTERM, shutting down

what the information state and a follow up action to take to troubleshoot the error.

I also attached a link to the error.log file http://www.mediafire.com/?3zeuz48eeu8g8v1


texts: Geoserver WMS Circle SLD style return incomplete Circle with higher zoom levels

texts:
I'm using OpenGeo Suite Boundless v2.6

I created a SLD style for a points layer and create 2 rules:

Rule to draw a graphics icon from png file.
Rule to draw a circle with radius(size) from db table using (OGC Attribute) value in meters.

I use the following code to load the layer in OL3:

    var wmsSource = new ol.source.TileWMS ({
    url: 'http://localhost:8080/geoserver/wms',
    params: {'TILED': true, VERSION:'1.3.0',
       'LAYERS': 'myds:antenna',
    },
    serverType: 'geoserver'
});

new ol.layer.Tile({
         id:'Message',
         title: 'Antennas',
         visible: true ,
         source:wmsSource
     });

The layer loaded correctly on the map in some lower zoom levels (0-10) if i zoomed in to more levels down (11-17) the circle shape is trimmed by the tiles like the image below:

The strange thing when i set the 2nd rule in the SLD that to get a fixed radius value like (2000) meters - Not from db table -, it draw a complete circle in all zoom levels, but i need to set a different radius for each circle.

I tried many solutions to avoid this error by changing the metatiling factors to be 100x100 instead of 4x4 and change the Gutter to be 100px, it solve the problem but affect the loading speed badly.

So kindly help me solve this problem.


texts: Complete QGIS Uninstallation to fix TauDEM Issues

texts:
I was looking to completely uninstall QGIS 2.18.14 from my laptop as I have managed to do something with the settings that has stopped the TauDEM provider from creating any outputs from any of the processes I have attempted to run.

The steps i have taken are as follows-
- Run QGIS uninstaller
- Deleted QGIS2 folder from C:\Users\me.qgis2
- Run Revo Uninstaller for removing trace files and folders

The problem is that somewhere there is a setting or folder remembering my last working configuration as when i reinstall QGIS the plugins, WMS layers and Panels that I have had opened for my daily work, are still opening as they were before even after a fresh install.

Is there something I am missing to have a fresh install with all the default tools and pathways reset to as they were if I had never had a copy of QGIS on my machine before?

I believe the problem may have arisen when i tried to install a newer version of SAGA GIS onto my machine but i can't be 100% certain as I have only noticed this TauDEM problem over the weekend. No matter what I do it always says the 'The following layers have not been created'. I have removed the newer version of SAGA but it is still not allowing me to create anything using TauDEM provider. the version of TauDEM is 5.3.7 x64 for Windows. I am currently running Windows 10 on my laptop.


texts: Testing Oracle query produces valid results in QGIS Python Console

texts:
I have a query that searches an Oracle database using a where clause to get the geometry of an object and then zooms to the extent of that object. The value for the where clause is retrieved from a dialog that opens after a button press of a custom tool button.

I want to be able to test if the input in the where clause for the query matches what is in the database, and if not to produce an error message stating that the input was not recognised.

I have created a connection to Oracle and I can zoom to the extents of the object, but if the query does not match any results it zooms to the coordinates 0,0.

I have tried the following code and made several variations on it but the query always returns as valid, whether the input matches rows in the table or not.

uri = QgsDataSourceURI()
uri.setConnection("provider","port","database","user","pass")
uri.setDataSource("user","table_name","geometry_column","where_clause")
try:
    layer = QgsVectorLayer(uri.uri(),"layer_name","oracle")
        if layer.isValid():
            qgis.utils.iface.mapCanvas().setExtent(layer.extent())
            qgis.utils.iface.mapCanvas().refresh()
except:
    qgis.utils.iface.messageBar().pushMessage("Error","The input was not recognised.", level = QgsMessageBar.CRITICAL, duration = 15)

I tried the code from the accepted answer here but this gave the results of query.value(0) as None and the db parameters set were all empty, despite the values being explicity set. I don't know whether this is because it is Oracle and not PostgreSQL.

Is there a way to test that the database query matches rows from the table before acting upon the results?


texts: Saving openJUMP data to be displayed in leaflet

texts:
Using OpenJUMP, I've managed to create a road map layer split at every intersection of the lines (using the noder tool).

I need to take this data and display this using leaflet by using end coordinates of each of the split lines. How do I save the layer which has been split?


texts: Where is pyCSW metadata repository stored in GeoNode?

texts: I have a working instance of GeoNode 2.0 over Debian7 (with PostgreSQL) and I want to know where does pyCSW stores the metadata by default. And if it's possible, how to set it to the same PSQL server.


texts: Length difference between Google Earth and ArcGIS for Desktop?

texts:
I have drawn a complex path with multiple vertices in Google Earth (GE) Pro and its measurement property is shown as being 2,460 nm. It is located East of Australia.

If i import this KML into ArcGIS and try to measure it, i get different results:
Using the Calculate field and the formula !shape.length@nauticalmiles! using WGS84 GCS, i got 2,554 nm.

Thinking that i should better use an equidistant projected coordinate system to get the length of my path, I redid the measurement and got 1,966 nm with the World Equidistant Cylindrical PCS.

I also played with the Plate Carree PCS and got 2,554 nm.

Am i wrong to use an equidistant projection for this calculation as it is normally supposed to maintain the distance?

Can i also consider GE length as being the correct value? I used external tool such as KML path measurer and got 2,467 nm.

What would be the correct process to do such a basic calculation?


texts:
We are publishing raster layers (computed via kriging) to ArcGIS server 10.1. A significant feature of our local geography is a large body of water, in our case, Lake Michigan.

Although I can clip the raster layer to exclude the water on ArcGIS Desktop, I get an error message: 000026 Data frames with a clip shape are not supported.

Because of this, the raster layer "bleeds" into the water. The result is not pleasant.  How can I effect "clipping" and still publish to the ArcGIS server?

P.S. I considered some kind of a mask, to block out the lake. But it seems that the transparency setting is for the whole layer, and so might block out the lake, but also cover up the land.)


texts: Generating Points for each Shortest Path

texts:
Python code is:

import arcpy

arcpy.FeatureToLine_management("Closestroute", "Featureline", "0.0001 Meters", "ATTRIBUTES")

#Generate points along lines for every 1 km
arcpy.GeneratePointsAlongLines_management("Featureline", "splitPoints", "DISTANCE", "1000 Meters", "", "NO_END_POINTS")

#Split line at point for every 1 km roads,search radius = 20 meters
arcpy.SplitLineAtPoint_management("Featureline", "splitPoints", "split1kpoints", "20 Meters")

#I added SEE field into splitPoints, then I assigned name for each rows.(s1, s2, s3, ...)
cursor=arcpy.da.UpdateCursor("splitPoints","SEE")

shpFld = arcpy.Describe("splitPoints").shapeFieldName
arcpy.DeleteIdentical_management("splitPoints", shpFld)

i = 0

for row in cursor:
    row[0] = 's' + str(i)
    i += 1
    cursor.updateRow(row)

del row
del cursor

My aim:

I am trying to obtain which start points are used to arrive in destination by using shortest path.

My problem is:

In table below, I have OBJECTID for each points. But, when more than one routes intersect with other route like in the picture below (yellow arrow), although those points are used by two routes, in table, it assigned for just one shortest path. (In table, please check yellow arrow) ex: ORIGIN id 10078 is just assigned to location3-location2. (from location 3 to 2)

I just want to get which points are used for each shortest routes and see which points are used by more than once.

How can I do this ?


texts: Why is the order of my features changed when using readFeatures in OpenLayers v4.6.5?

texts:
I have a geojsonObject that is a featurecollection of X number of features.
see picture

I Then add this to a vector:

var Source = new ol.source.Vector({
        features: (new ol.format.GeoJSON()).readFeatures(geojsonObject)
    });

I add this to a layer, and call a style function:

var Layer = new ol.layer.Vector({
        source: Source,
        style: stylefcn(Source.getFeatures())
    });

function stylefcn(_features) {
console.log(_features)
//When i log my features here, the order of the features have changed.
}

Why does the order change?

E.g: Feature number 0 from the picture above, has become feature number 4 in the style function.

(I have a small suspicion that this might happen when the number of features pass 10)


texts: Getting feature class from File Geodatabase in ArcObjects using VB.NET?

texts:
How do you open up a File GDB and read out a simple feature class?

The snippet for retreiving a shapefile from a regular folder on disk works, and I suppose I can use the same code but with a different kind of workspace.

''' <summary>
''' Get the FeatureClass from a Shapefile on disk (hard drive).
''' </summary>
''' <param name="string_ShapefileDirectory">A System.String that is the directory where the shapefile is located. Example: "C:\data\USA"</param>
''' <param name="string_ShapefileName">A System.String that is the shapefile name. Note: the shapefile extension's (.shp, .shx, .dbf, etc.) is not provided! Example: "States"</param>
''' <returns>An IFeatureClass interface. Nothing (VB.NET) or null (C#) is returned if unsuccessful.</returns>
''' <remarks></remarks>
Public Function GetFeatureClassFromShapefileOnDisk(ByVal string_ShapefileDirectory As System.String, ByVal string_ShapefileName As System.String) As ESRI.ArcGIS.Geodatabase.IFeatureClass

  Dim directoryInfo_check As System.IO.DirectoryInfo = New System.IO.DirectoryInfo(string_ShapefileDirectory)
  If directoryInfo_check.Exists Then

    'We have a valid directory, proceed

    Dim fileInfo_check As System.IO.FileInfo = New System.IO.FileInfo(string_ShapefileDirectory + "\" + string_ShapefileName + ".shp")
    If fileInfo_check.Exists Then

      'We have a valid shapefile, proceed

      Dim workspaceFactory As ESRI.ArcGIS.Geodatabase.IWorkspaceFactory = New ESRI.ArcGIS.DataSourcesFile.ShapefileWorkspaceFactoryClass
      Dim workspace As ESRI.ArcGIS.Geodatabase.IWorkspace = workspaceFactory.OpenFromFile(string_ShapefileDirectory, 0)
      Dim featureWorkspace As ESRI.ArcGIS.Geodatabase.IFeatureWorkspace = CType(workspace, ESRI.ArcGIS.Geodatabase.IFeatureWorkspace) ' Explict Cast
      Dim featureClass As ESRI.ArcGIS.Geodatabase.IFeatureClass = featureWorkspace.OpenFeatureClass(string_ShapefileName)

      Return featureClass
    Else

      'Not valid shapefile
      Return Nothing
    End If

  Else

    ' Not valid directory
    Return Nothing

  End If

End Function


texts: Unexpected vertices from the ouput of ArcGIS union tool

texts:
I have found some unexpected behavior of ArcGIS union tool. My assumption was that this tool will never add any new vertex to the output layer, but this was not the case.

I have run the ArcGIS union tool on only one layer GADM (world admn01) and found that in output layer there are around 866 new vertices. I think this tool make tiles of data and then process. After processing there are extra vertices which are basically created because of tiles are not completely deleted.

Is there anybody else with the same problem/experience?

I am using ArcGIS 10.0 on Windows7 64bit


texts: Sum attributes from one table and add to another in QGIS

texts:
In QGIS 3, I have two layers:

Polygon layer representing properties including an attribute on the type of property: residential, commercial, industrial, open space etc.
point layer representing connections to watermains

I have used the NNJoin plugin to determine the closest point for each of the properties, and this reference is now in the property layer.

For each point, I now want to sum the number of properties of each type and add this to point layer.

For example:

Point  Residential Commercial Industrial etc
 1        3           0           4
 2        7           2           5
 ...
 5000     5           5           15

I was going to export it to Excel to undertake a series of "Sum Ifs", but didn't want to take it out of GIS unless I had to.

Is there a method to achieve this in QGIS?


texts: Global NetCDF datasets with station time series data?

texts:
I´m working on a web application which uses temperatur and precipitation data. Currently i´m using only globally interpolated, multidimensional rasterdata from terrestrial stations (queried with x,y,t), but i also want to use global datasets with values of the original time series (which means vector- instead of raster data). So far i didn´t found many of these, although the NetCDF data model seems to be well suited for this case of data storage.

I want to use NetCDF because i need a standardized way of querying the data. That means i don´t want to use datasets which i hvae to convert or edit before i can use them. That´s why i would prefer NetCDF over simple ASCII or csv files...

Does any one know of possible Websites/Portals/Sources where i maybe could find the data i need in this or similar formats following the common data model (e.g. HDF5)?


texts:
import arcpy
import csv
import glob, os
os.chdir("C:\\path\\folder")
fieldsource = "C:\\path\\sampleshp.shp"
fields = [r.name for r in arcpy.ListFields(fieldsource)]
csv_out = "C:\\path\\table.csv"
with open(csv_out, "wb") as f:
    wr = csv.writer(f)
    wr.writerow(fields)
    for file in glob.glob("*.shp"):
        in_shp = "C:\\path\\folder\\" + file
        with arcpy.da.SearchCursor(in_shp, fields) as cursor:
            for row in cursor:
                wr.writerow(row)

I have few thousand shapefiles with only a few rows in them. I want to export their table to a single csv file. The script works but it seems I have a character in some of the files that python(2.7) doesn't like:

  Traceback (most recent call last):
    File "C:\path\folder\table.py", line 20, in
      wr.writerow(row)
  UnicodeEncodeError: 'ascii' codec can't encode character u'\xe1' in position 2: ordinal not in range(128)

How can I exclude/change the character that creates the issue?


texts: Read only a crop or certain extent of a tif raster R

texts:
I am working with a huge amount of data in daily tif files (thousands of daily files). I am analyzing the mean of the raster in a shapefile area, repeated over potentially thousands of shape layers.

My current .tif files are for an entire country when in reality I only need a small area of the of each tif file for each shapefile layer. Each shape layer has a different set of the days to extract from, i.e. data frame like this:

df <- data.frame(shplayer=c("shp_layerbuffer1","shp_layerbuffer2", "shp_layerbuffer3"), start=c("2000_02_18", "2004_03_19", "2010_08_20"), end=c("2010_01_09", "2005_03_15", "2012_09_04"))

Is there a way, in R, to crop a .tif (or any raster type file) BEFORE reading the file? This way I could read just the cropped area of each of the tif files

I envisoned something like (repeating across the entire set of dates):

library(sf)
library(raster)
shp_layerbuffer1 <- st_read("myshpfolder", layer="shp_layerbuffer1", quiet=T)

###EXAMPLE BUT doesn't work to crop the raster as it comes in
tempraster <- raster::raster("mytif_2000_02_18.tif", ext=extent(shp_layerbuffer1))

Then the usual velox (or raster) extract, repeat.

Is this possible?


texts: How to produce high-resolution DEMs from new Sentinel-1 data?

texts:
This is a general question relating to the new Sentinel-1 data. I'm interested in the potential for using this data to create high resolution DEMs in GIS (I use QGIS), to investigate archaeological sites and landscapes and model their change.

The European Space Agency Copernicus program satellite Sentinel-1a has recently been commissioned, with data just beginning to become available on the ESA Data hub website (https://scihub.esa.int/). This data can be viewed and processed using specially prepared software called Sentinel-1 Toolbox (https://earth.esa.int/web/sentinel-tbx/sentinel-1-toolbox). The radar data will allow land surface observations and changes in land surfaces to be modelled using 'synthetic aperture radar interferometry', with most examples given so far detailing the impact of earthquakes or other large-scale natural landscape changes (https://earth.esa.int/web/sentinel/thematic-content/-/article/radar-vision-maps-napa-valley-earthquake). After searching the internet I have found that, other than the ESA's own website, there appears to be little information on the practical use of this new data.

Is there any literature on
how to produce high-resolution DEMs from this new Sentinel-1 data?


texts: How to get mean from image collection in numerical format?

texts:
I have some GEE code that works fine in the code editor and it plots a chart like this: It has this image collection called 'toplot' and you pass it to the UI chart function and it plots the chart for you. It has 45 bands, you pass the region of interest and means of all images in the collection are plotted in the chart. This is the code that does the job:

var annual_prof = ui.Chart.image.series({

imageCollection: toplot,

region: ROI,

reducer: ee.Reducer.mean(),

scale: scale,

xProperty: 'system:time_start',

});

annual_prof.setChartType("ColumnChart");

annual_prof.setOptions({

title: (Adm+'-Avg NDVI Anomalies in the 12 months\nto: '+lastdatetext),

hAxis: {format:'', title: '8 days period'},

vAxis: {title: 'Anomaly %'},

lineWidth: 2,

pointSize: 1,

legend: 'none',

series : {

0 : {color : 'darkgreen'}}

});

print(annual_prof)

I need to get those means in an array form so that I can use it to plot the chart on the webapp. The UI functions don't work in Python and JavaScript APIs. So after reading the documentation I tried this:

function get_mean(d){
  return d.reduceRegion(ee.Reducer.mean(), ROI);
}

var nums = toplot.map(get_mean);
print(nums)

But this gives me an error:

  Collection.map: A mapped algorithm must return a Feature or Image.

How can I get these numbers as an array or any other numerical form using GEE functions?


texts: Comparing values between two fields in QGIS Graphical Modeler

texts:
I'm building a graphical model in QGIS 3.4 to check if a point feature with a specific value in the A field is within a polygon with the same value in B field.
The image displays the draft of the model.

I set two vector layers as input (point layer and polygon layer), then I added the algorithm "Join attributes by location". After I tried to use "Select by attribute" to compare the attribute values between two fields and I continued adding "Extract selected feature" algorithm.
My issue is in "Select by attribute" part, I don't know what I need to write in the "Value" field or if there are other solutions to compare the value of the two fields.

I found this post that uses actions but I prefer Graphical Modeler.


texts: Add data to attribute problem in PyQGIS?

texts:
Based on xunilk code in previous post I improved it to add mean distances to "MeanDist" field, while I can print mean distance for each point id, but this data didn't write in attribute of layer. what is the problem for add data to the attribute?

from PyQt4.QtCore import *
import itertools
import numpy as np

layer = iface.activeLayer()

feats = [ feat for feat in layer.getFeatures() ]

n = len(feats)

comb = range(n)

distances = [ [] for i in range(n) ]
indexes = [ [] for i in range(n) ]

for i, j in itertools.combinations(comb, 2):

    dist = feats[i].geometry().distance(feats[j].geometry())

    if dist < 500:
        i_dist= distances[i].append(dist)
        i_index= indexes[i].append([i,j])
        j_dist= distances[j].append(dist)
        j_index= indexes[j].append([i,j])

prov = layer.dataProvider()
layer.startEditing()

if prov.fieldNameIndex("MeanDist") == -1:
    prov.addAttributes([QgsField("MeanDist", QVariant.Double, "double",    10, 2)])
else:
    pass

for feature in layer.getFeatures():
    attrName = "MeanDist"
    for i, group in enumerate(distances):
        if i == feature.id():
            print feature.id(), np.mean(group)   # print for test results, it must remove.
            prov.changeAttributeValues({feature.id() : {prov.fieldNameMap()[attrName]: np.mean(group)}})
            layer.updateFeature(feature)

layer.updateFields()
layer.commitChanges()


texts: Removing Intersection of Buffers in Same Layers

texts:
I have created a series of buffer regions around line segments in QGIS (in my case a road network):

The buffers are in transparent blue with the darker circles being where the buffer segments overlap. Additionally, the red outline shows the shape of a single buffer polygon.

I'm trying to remove the areas in my map where these buffers intersect. Specifically, I would like a layer composed of all of the non-overlapped polygons, and to keep their attributes (in my picture above it would be all of the lighter-blue polygons). Most of the geoprocessing tools need two separate layers. Also, while the Topology Checker can identify the overlapping areas in the same layer, it cannot select them.

Is there a straightforward in QGIS or GRASS that can remove the overlapping areas of buffer polygons in the same layer?

UPDATE

Following Kazuhito's response, I ran the Union command in QGIS 3.2, but it did not produce the desired output:

From from Union changelog, it should be that union "splits off" this overlapped area from the original buffer polygon. Thus the highlighted intersection should only have 2 polygons (the overlaps for each layer), but instead it has 4 because it retains the original polygon as well. Does anyone know why the Union command would not produce a split up polygon as illustrated in the changelog?


texts: Endmember selection when applying SAM, SVM and MLC to PCA

texts:
I intend to apply the Spectral Angle Mapper (SAM), Maximum Likelihood (MLC) and Support Vector Machine (SVM) algorithms on a hyperspectral dataset in ENVI but I am slightly confused as to how to derive the endmembers (different land covers in this case) after applying Principal Component Analysis. I have performed radiometric correction and a forward PCA rotation which resulted in 3 PCA bands. However, I am slightly confused regarding endmember selection procedures when using PCA.

For example, for the MLC algorithm I understand that if I wasnt using PCA I would simply create new Regions of Interest to select the endmembers from the original radiance image.

However, when using PCA as a dimensionality reduction technique, I am not sure where to derive the endmembers from. i.e is it from the original radiance image or from the PCA components?


texts: Making map use set of coordinates in QGIS

texts:
I have made a custom google map with the google map api which has a few functions that work well within it - one specifically is limiting the user on what they can view within the bounds of the map.

The map is a very simple 2 sides the front and back of the human body.

What I want to do is add on specific zones in my body, I've read that I can do with via geojason and have been pointed towards qgis to generate my GeoJSON data.

Now with is where I fall down I can make the zones and even the export just the coordinates don't match up. How do I make the coordinates match?

I know to use EPSG:4269, but the coordinates still don't match.

How do I make it match i.e. set a specific long/lat that my map is using?


texts: Does FME retain a log of the last run process

texts:
The hard drive where my FMW was has crashed and I was wondering if there is somewhere in the user profile or in the install location of FME that retains a log of what processes were run in the FMW?

In FME I can see the recently opened fmw files on the crashed drive but I can't access the log file that's usually present in the same folder as the fmw.


texts: Checking if a specific leaflet.draw control is enabled?

texts:
Where exactly are the leaflet controls stored that are added to the map?
So at anytime in the app, I can check something like

if(controls.polyline.enabled==true)

My guess is its somewhere in the map.layers. But theres a lot of layers there. Any way to identify which one is our control layer? and then the attribute which checks if it is currently enabled or disabled.


texts: Keep "no data" tiles as "no data" w/ proximity (raster distance) operation

texts:
Working in QGIS 2.8.1 Wien.

I have a raster file with three classes:

"No data" for tiles outside of the country of interest's bounds
"0" for tiles that correspond to roads
"1" for all other tiles.

I want to use proximity (raster distance) to create a new raster layer of distance to roads for all tiles... but want to maintain the "no data" tiles as "no data".

However, when I run proximity, it ignores any calls to the "No data = " option and assigns distance values to all "no data" tiles. I can't type any string input into the "No data = " option, and am unsure of how to keep my "no data" cells.

How do I do this?


texts:
I am trying to open a shapefile with Jupyter notebook, but I can't get it to open. The shapefile is located in the same folder as my workbook, and Python can't identify it even if I use the full path.

import shapefile as shp
import matplotlib.pyplot as plt

sf = shp.Reader(r'C:\Users\Public\Documents\1-11-99n.shp')

ShapefileException: Unable to open C:\Users\Public\Documents\1-11-99n.dbf or C:\Users\Public\Documents\1-11-99n.shp.

Opening it with fiona and geopandas does not work either.

import fiona
shape = fiona.open("1-11-99n.shp")

and

import geopandas as gp
shp = gp.GeoDataFrame.from_file('1-11-99n.shp')
print (shp)

both end up with

CPLE_OpenFailedError: b'Unable to open 1-11-99n.shx or 1-11-99n.SHX.Try --config SHAPE_RESTORE_SHX true to restore or create it'

I know this might be a trivial issue, but I am new to using shapefiles and spatial Python libraries so I really don't know how to work around this issue.


texts: Join Attributes By Location inside QGIS modeler still fails to compute sums

texts:
Recently, a question was asked about the spatial join tool, and its failure to calculate the desired data.

I ran into the same problem, but none of the solutions that were offered there solved it.

I ran the tool as a standalone process, and it worked fine
I made sure both my layers (a polygon and a point) are the same CRS, still failed
downgrading is not an option

Running 64bit standalone 2.12.1 Lyon on Windows 7


texts: What resampling method is used in gdal_merge?

texts:
I am using the gdal_merge python script to merge rasters that have slightly different, though similar, input resolutions. The documentation advertises that it can handle different resolutions, and in most cases, this works out reasonably well. However, I have had a few cases where the final merged raster looks significantly degraded, despite a relatively small resolution change (~10%). I first looked through the gdal_merge documentation to see if there were any options to set the interpolation method used to resample rasters whose resolutions change, but did not see anything. I then realized I do not even know what resampling method gdal_merge is using. How does gdal_merge handle resampling?

Also, if anyone has any thoughts on how to improve the results from using gdal_merge on rasters with different resolutions, I would be happy to hear those. I figured that would be a more challenging question though without providing many more details.


texts:
Maybe this could be a duplicated, but I can't find the right answer.

I'm using EPSG:25830 to default projection of my map. For this I've used proj4js :

proj4.defs("EPSG:25830", "+proj=utm +zone=30 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs");

var bounds25830 = [712397.500, 4374797.500, 742202.500, 4394802.500];
var proj25830 = ol.proj.get('EPSG:25830');
proj25830.setExtent(bounds25830);

Now I'm declaring some WMS Layers and adding to my map :

var ignBase = new ol.layer.Tile({
    name: 'IGN Base',
    visible: true,
    source: new ol.source.TileWMS({
        url: 'http://www.ign.es/wms-inspire/ign-base',
        crossOrigin: 'anonymous',
        params: {'FORMAT': 'image/png',
                'VERSION': '1.1.1',
                transparent : false,
                tiled : true,
                LAYERS: 'IGNBaseTodo',
                STYLES: '',
        }
    })
});

var layerVectorVacia = new ol.layer.Vector({
    title:'Vacía',
    name: 'base'
});

var groupCapasBase = new ol.layer.Group({
    name: 'Capas Base',
    layers: [layerVectorVacia, ignBase]
});

var mousePositionControl = new ol.control.MousePosition({
    coordinateFormat : ol.coordinate.createStringXY(2)
});
var map = new ol.Map({
    layers  : [
        groupCapasBase
    ],
    target  : 'map',
    controls: ol.control.defaults({ attribution : false }).extend([mousePositionControl]),
    view    : new ol.View({
        projection: proj25830,
        zoom      : 5,
        center    : [727300.000, 4384800.000]
    })
});

In the map, nothing is rendered, and if I look at the network panel, it appears that images have not been loaded. My first thought was that the request is badly formatted, but if I copy and paste some of the failing request in a browser, it returns the right image.

EDIT

  I've noticed that if I use "FORMAT" : "image/jpeg" intead of
  "image/png" is rendering properly.


texts: How do I automatically generate a spatial data catalogue from shape files and .gdbs?

texts:
I've recently inherited a whole bunch of spatial data (mostly shapefiles but some .gdbs) which relate to a bunch of different cities. I'm looking for a way to take an inventory of what data I've got for each city, preferably automated!

Does anyone know of any (free) tools that would allow me for example, to create a catalogue of my data (and metadata) from the files contained within my directory structure?

I'm using Linux (Ubuntu 14.04) and tend to use QGIS/PostGIS.


texts: QGIS OpenLayers Plugin shows black for OSM?

texts:
the QGIS OpenLayers Plugin version 1.4.3 shows all OSM WMS Layers broken with black- or redscreen "tiles". A colleaque had version 1.3.6 installed and it worked fine. I was already on 1.4.3 with broken layers. With the new version (1.4.2?) OSM layers was split up in two groups (OpenStreetMap, OSM/Thunderforest). After updating to 1.4.3 my colleaques OSM layers were also broken. We both use QGIS 2.18.0 on Windows7-64bit machines.

I searched Google, the Github-Repository and StackExchange but I can't find a bug report on this issue.


texts: Importing multiple shapefiles into one single Postgres table

texts:
I am trying to import in one single table different shapefiles using the command:

 shp2pgsql -D -s 4326 -c my_shapes_table -I *.shp  > my_shapes_table.sql

I have tried different combinations of the parameters but what I got are always multiple tables with the name of each single file. Any idea how to get one single table only that imports all the files?


texts:
In this link xml-link, there is a .xml file from a website (WFS) which contains vector data regarding water bodies for a specific area. I saved it as .xml by just right-clicking -> "Save as" and I tried to load it in QGIS as WFS (Open Data Source Manager (it is the icon below blank page on the top left of the window for QGIS 3.2 version) ->  WFS -> Load -> Selecting the .xml file which I saved before. However, I have the error which you can see in the image.

I tried a workaround and it worked out but I do not know if I missed information. I right-clicked on the WFS from the Table of Contents -> New Connection -> I put a name (it does not matter) and the link:  https://rijnland.enl-mcs.nl/arcgis/services/Legger_Oppervlaktewater_Vigerend/MapServer/WFSServer? from the .xml file in the URL position.
 ).

There are many links in the .xml file but I resorted to it by trial and error. It worked and I can see five layers. Are all these layers that the current WFS has to offer or did I miss any information?


texts: Building error on enum of the QgsMapToolIentfy Class

texts:
I am working on a project in c++/QT4.8.5 based on the api c++ 2.8.2 of QGis
It is a stand alone application

I had to get information of a layer by clicking on it in the canvas

So I decide to use the class QgsMapToolIdentify ( i saw some sample in the book building mapping application with Qgis from Erik Westra )

here is my code :

mIdentify =new QgsMapToolIdentify(mpMapCanvas);
QList <struct IdentifyResult *> lResultList;

lResultList = mIdentify->identify(lCoord.x(),lCoord.y(), TopDownStopAtFirst, VectorLayer);

but when I want to build it i have an error on the enum value TopDownStopAtFirst and VectorLayer

the error message :

  src/IHM/mainwindow.cpp:174:83: error: ‘VectorLayer’ was not declared in this scope
        lResultList = mIdentify->identify(lCoord.x(),lCoord.y(), TopDownStopAtFirst, VectorLayer);

and the same message for the TopDownStopAtFirst

I have included the header qgsmaptoolidentfy.h in main.cpp.


texts: How to align directions of connected lines?

texts: I have a vector line made of different pieces of lines with different directions. I merged the lines to one. But this merged line still have different direction in its "invisible" parts. I want to plot a profile, so I need ONE direction. Any idea?


texts: Performance issues Tilemill with PostgreSQL and large datasets

texts:
I have a Tilemill project the size of europe. I developed my style-sheet by using a smaller dataset (size auf Austria) and everything was perfect. But now I did apply the same styles to my big europe dataset which is in a PostgreSQL database and it's not rendering anything.

I had a look at what takes so long and it seems as if it is PostgreSQL. I already added all kinds of indexes to my tables, but that doesn't seem to help.

As a starting point I used OSM-Bright and that added a layer called landuse-gen0. It pulls in all the polygons of a certain size.

SELECT way, way_area AS area, COALESCE(landuse, leisure, "natural", highway, amenity, tourism) AS type
  FROM planet_osm_polygon
  WHERE way_area > 100000
  ORDER BY way_area DESC;

This query result in 3,4M rows. Each has about 10k... which would result in about 34 GB of data just for this one query. But Tilemill is doing a lot of those queries.... And tilemill seems to also be limiting the Bounding-Box AFTER my query... which is pretty inefficient, isn't it?

Any way I can improve my performance? Cause right now nothing is rendered... even if I wait one hour... not even one tile...

Btw. the computer it's running on has 16 cores and 120GB RAM... so hardware shouldn't be the big problem here...

Thanks for any tip!


texts:
Context:

I've imported a DEM of residential composed in several buildings.
I've identified the max height of buildings by addind points in a layer for every building. Then I've integer an attribute field "max_height" and the value is the absolute altitude of the max Z of the building.

Aim :

For every building (around 150) I have to make a map that shows any object on the site higher than the max Z of the building (max_height).

My need:

I'm very new into GIS and the only way I know to exploit a DEM is to make a colored map.
That's why I would like to create colored maps where the color scale is based on a different height for building and the min value is the max z of the building (max_height), for every building.**

How can I automate this task?


texts: Configuring QGIS 2.18 processing script locations

texts:
I've just installed QGIS 2.18. I used the ltr version from from the advanced options of the OSGeo4W network installer. I went to set the folders where I keep my python scripts but when I click 'ok' I get the following error...

Has anyone had anything similar and know how to resolve?

An error has occurred while executing Python code:

RuntimeError: maximum recursion depth exceeded in cmp
Traceback (most recent call last):
  File "C:/OSGEO4~1/apps/qgis-ltr/./python/plugins\processing\gui\ProcessingToolbox.py", line 163, in updateProvider
    item.refresh()
  File "C:/OSGEO4~1/apps/qgis-ltr/./python/plugins\processing\gui\ProcessingToolbox.py", line 377, in refresh
    Processing.updateAlgsList()
  File "C:/OSGEO4~1/apps/qgis-ltr/./python/plugins\processing\core\Processing.py", line 184, in updateAlgsList
    Processing.reloadProvider(p.getName())
  File "C:/OSGEO4~1/apps/qgis-ltr/./python/plugins\processing\core\Processing.py", line 189, in reloadProvider
    algList.reloadProvider(providerName)
  File "C:/OSGEO4~1/apps/qgis-ltr/./python/plugins\processing\core\alglist.py", line 55, in reloadProvider
    p.loadAlgorithms()
  File "C:/OSGEO4~1/apps/qgis-ltr/./python/plugins\processing\core\AlgorithmProvider.py", line 56, in loadAlgorithms
    self._loadAlgorithms()
  File "C:/OSGEO4~1/apps/qgis-ltr/./python/plugins\processing\modeler\ModelerAlgorithmProvider.py", line 73, in _loadAlgorithms
    folders = ModelerUtils.modelsFolders()
  File "C:/OSGEO4~1/apps/qgis-ltr/./python/plugins\processing\modeler\ModelerUtils.py", line 46, in modelsFolders
    folder = ProcessingConfig.getSetting(ModelerUtils.MODELS_FOLDER)
  File "C:/OSGEO4~1/apps/qgis-ltr/./python/plugins\processing\core\ProcessingConfig.py", line 209, in getSetting
    if name in ProcessingConfig.settings.keys():
RuntimeError: maximum recursion depth exceeded in cmp

Python version: 2.7.5 (default, May 15 2013, 22:44:16) [MSC v.1500 64 bit (AMD64)]
QGIS version: 2.18.19 Las Palmas, 3a17f72ba5


texts: QGIS Clipping raster with vector-getting all touching pixels

texts: Is there a setting I can change in order to get all the pixels of the raster layer which are touched by the vector boundary in the clipped image? Right now, only those pixels which lie completely within the boundary are included, and it's causing issues.


texts: Performing Raster Calculator functions using Python and open source modules?

texts:
I would like to perform mathematical functions similar to those you would use Raster Calculator for in ArcMap, however I'd like to do it in Python without using the ArcpPy module. I have the osgeo module, and it seems like it should be possible to do however I'm having a hard time finding any useful syntax.

An example function would be to multiply every cell value in an input raster by 2.


texts: How to enter text into attribute tables with QGIS graphical model?

texts:
I am using QGIS 2.14 to open NC files for chlorophyll satellite data with the end goal of compiling chlorophyll values for many dates (each its own initial file) in a spreadsheet. I know how to do the work manually, but I want to automate the process in the model due to the bulk of files I am working with.

Currently, my process is to manually open each file in QGIS, then I have been using the model. I run zonal statistics for each file based on the study area, then I add a date attribute in each file. My problem is that I can't figure out how to fill in the date values (using the format 2000_01, 2000_02, etc) using field calculator or advanced python field calculator. I need the date attributes in order to merge the files for export.

I think this would be possible in ArcMap, so can I do it in QGIS with minimal/no python?

I've included the model that I've started with.


texts: SQL: Select aggregated Bounding Box of multiple Features

texts:
I would like to calculate the aggregated bounding box of multiple features using OGC Simple Feature SQL or SQL API of ArcSDE. The following statement will only query the bounding box of each feature:

select ST_Boundary (shape) from polygon_table where country like '%usa%'

What I needed is one bounding box aggregated over all features, resulting to the condition in the where clause.


texts:
My AddJoin operation is failing:

arcpy.CopyRows_management(lup_file,out_dir+'lup')
try:
    # Perform copy rows on merged csv file
    print 'start'
    aa = arcpy.Raster(SEIMF)
    arcpy.AddJoin_management(aa, "VALUE", out_dir+'lup', "VALUE", "KEEP_ALL")
    print 'end'
    #arcpy.AddJoin_management(gis_dir+SIF, "VALUE", out_dir+'lup', "VALUE", "KEEP_ALL")
except:
    print arcpy.GetMessages()

error messages:

  Executing: AddJoin C:\Users\Documents\Projects\AB\Inputs\gis\sif VALUE C:\Users\Documents\Projects\AB\Outputs\C\output\lup VALUE KEEP_ALL
  Start Time: Sat Feb 22 02:27:50 2014
  Failed to execute. Parameters are not valid.
  ERROR 000840: The value is not a Table View.
  ERROR 000825: The value is not a layer or table view
  ERROR 000840: The value is not a Raster Catalog Layer.
  ERROR 000840: The value is not a Mosaic Layer.
  WARNING 000970: The join field VALUE in the join table sif is not indexed. To improve performance, we recommend that an index be created for the join field in the join table.
  Failed to execute (AddJoin).

Does anyone know what is happening? I get this error only in the python code and can execute AddJoin just fine if I use the ArcGIS manually.


texts: FME in task scheduler fails to get ESRI license

texts:
We have FME2016 and an ArcGIS 10.3.1.  I am trying to run a task in task scheduler that writes to either a file geodatabase or an ArcSDE enterprise database. This has been successful until about 2 weeks ago when the ArcGIS license manager was updated to ArGIS 10.3.1. from 10.2.1. and FME started showing the following messages in all of the FME log files during a task scheduler run:
"The Esri Product license used by this machine does not include licenses for the 'Engine Enterprise Geodatabase' licensing level
The Esri Product license used by this machine does not include licenses for the 'ArcInfo' licensing level
The Esri Product license used by this machine does not include licenses for the 'Engine' licensing level
|ERROR |Unable to perform licensing related tasks. Please ensure that ArcGIS is correctly installed."
The ArcGIS Administrator does show that licenses are available and is pointing to the correct license manager. The task scheduler settings are set to 'Run with highest privileges'.
We are thinking this has something to do with permissions settings in task scheduler because the same FME workbenches run successfully when the application is open and manually run.  Attached is a screenshot of the task scheduler settings.
Has anyone else run into a similar issue?


texts:
I want to group and count features with the same code attribute value.

I can classify the data and enable feature count in the layer list but I cannot get counts in the attribute table.


texts:
The Data

Two sets of 1000 points A.geom and B.geom in PostGIS

The Problem

From Set A take the first point and find its nearest neighbour in Set B
Calculate the 3D distance between these two points
Save the result to the dist column of Set A

I've been testing out various python libraries, I thought I had it working in shapely just to discover it ignores the z-dimension.

I think SFCGAL is the way to go, but I am inexperienced, is this a simple query, can anyone help me with this?


texts: Layer with custom UI not working properly

texts:
I have a problem with QGIS 2.4 when introducing layers with a custom UI. More precisely, widgets such as CheckBoxes are not recognised as such (appear as Text Edit in Attributes) and the corresponding Checked value appears as NULL (although checked by default in the UI file).

When I am manually changing the specific field:

instudy_id = layer.fieldNameIndex("instudy")
layer.setEditorWidgetV2(instudy_id, 'CheckBox')

Everything is normal. The corresponding part of the UI file is:

Sorry if my explanation is poor, I am quite new to QGIS. Hope the screenshots help.

To conclude:

some options of my Widgets (defined in my .ui file) are overwritten by QGIS (e.g. "Checked" or "Unchecked" state of a QCheckBox)
others, which are introduced in the corresponding .py file (like QgsEditorWidgetWrapper.fromWidget( in study ).setEnabled(False) ) are not respected by QGIS.

What is the proper way to introduce a custom / customise an .ui file ? Is there anywhere detailed information/examples on the use of newer methods like QgsVectorLayer.setEditorWidgetV2, setEditorWidgetV2Config, QgsEditorWidgetWrapper ?


texts: qgis:deletecolumn in QGIS 3 graphical modeler

texts:
I want to use the algorithm qgis:deletecolumn in a graphical modeler.

How do I have to fill the case Fields to delete :

["field1","field2"]
"field1","field2"
field1,field2
[field1,field2]

And is there another way to say to QGIS : delete all fields except field1, field2, ... ?


texts: OD Cost Matrix doesn't generate all the lines

texts:
I have 177 points of origin and 177 of destination. There's should be 31.329 lines generated by OD cost matrix, but there's only 2077 lines. I've made sure that all the points are on the traversable network, and i already checked the network with topology. I've unchecked the oneway prohibition and hierarchy option, but it still only generated 2077 lines, and when i isolated all the 128 points which don't have lines and ran OD cost matrix, only 33 from 128 points now have lines. What should i do to fix this?

the pink points are all the 128 error points which don't have lines (blue selection mark) shown in picture below

I'm using arcmap 10.2 with network analyst extension.


texts: Shapely strange splits when splitting LineString and Polygon

texts:
I noticed a strange behavior while splitting a LineString with e Polygon: in the below image I have a Polygon (in green) and two lines. If I split the blue Linestring with the polygon the result of the split is a list of 3 LineString (as expected), but if I split the orange Linestring with the polygon the result is a list of 16 LineString. Instead, I expected only 4 linestring; what am I missing?

The polygon is the result of line_blu.buffer(<value>).intersect(line_orange).buffer(<value>) so I am sure that all the orange point you see in the green area are actually in that polygon.
Just as prove of that if I calculate the points outside the polygon with the following code
xy = [xy for xy in line_orange.coords if not polygon.contains(Point(xy))]

The result is

Below you can see two of the 16 line obtained splitting line_orange with the polygon:


texts: GEE Verdet (Temporal Segmentation Algo) Failure

texts:
I'm struggling to get the Verdet Temporal Segmentaion Algorithm to run (https://developers.google.com/earth-engine/apidocs/ee-algorithms-temporalsegmentation-verdet). I get vague "Internal Server Errors" when displaying the output.
As per the documentation I think I'm providing the correct style of input – 1 image per year.   I'm thinking this is a bug in the Verdet code, given that I have noticed other similar posts on Verdet.
Link to code: https://code.earthengine.google.com/fdf8bdd1ecac7551b2dfd89f52d7417c
Started a thread in GEE developer's forum as well: https://groups.google.com/g/google-earth-engine-developers/c/cORUs-lFyCs/m/2UWkTX_DAAAJ
Code used (everything works, except adding the verdet layer to the map):
function mask(img) {
  var clearBitMask = 1 << 1;
  var qa = img.select('pixel_qa');
  var masked_img = img.updateMask(
    qa.bitwiseAnd(clearBitMask)
  ).copyProperties(img, ['system:time_start'])
  return(masked_img)
}

function addNDVI(img) {
    var ndvi = img.normalizedDifference(["nir", "red"]).rename('NDVI')
    img = img.addBands(ndvi)
    return(img)
}

var start_yr = 2015
var end_yr = 2020

var coords = [
    [-74.31961, 43.94047],
    [-74.1584, 43.94047],
    [-74.1584, 44.07021],
    [-74.31961, 44.07021]
]

var aoi = ee.Geometry.Polygon(coords)

var yearlyComposites = ee.ImageCollection([]);

for (var i = start_yr; i < end_yr; i++) {

    var ls8Col = ee.ImageCollection('LANDSAT/LC08/C01/T1_SR').select(
            ['B2', 'B3', 'B4', 'B5', 'B6', 'B7', 'pixel_qa'],
            ['blue', 'green', 'red', 'nir', 'swir1', 'swir2', 'pixel_qa']
        ).filterMetadata(
            'IMAGE_QUALITY_OLI', 'not_less_than', 9
        ).filterMetadata(
            'CLOUD_COVER', 'not_greater_than', 33
        ).filterBounds(
          aoi
        ).filterDate(
          ee.Date(i + '-06-15'),
          ee.Date(i + '-8-30')
        )
    ls8Col = ls8Col.map(addNDVI)
    ls8Col = ls8Col.map(mask)
    ls8Col = ls8Col.map(function(img) {return(img.clip(aoi))});
    yearlyComposites = yearlyComposites.merge(ls8Col.median());
}

yearlyComposites = yearlyComposites
print(yearlyComposites.getInfo())
var centerCoords = aoi.centroid().coordinates()
print(centerCoords.get(0))
print(centerCoords.get(1))
var verdet = ee.Algorithms.TemporalSegmentation.Verdet(yearlyComposites.select('NDVI'));
print(verdet.getInfo())
Map.setCenter( -74.23900500000029, 44.00534468585141, 12)
Map.addLayer(verdet, {}, 'verdet')


texts:
I am currently attempting to put together a web mapping application for a client through the AGOL web app builder.  One of the layers in the application is a set of grid cells covering a large geographic area.  Another layer is a set of points containing a scoring value.  The client would like to be able to summarize the points that fall into each grid cell and symbolize the grid cells based on the aggregated point score for each cell.

I am aware there are tools in AGOL that can perform spatial analysis, but every time the analysis is run new data is created.  I would like a solution that updates and symbolizes the grid cells as the underlying point data changes.  Has anyone attempted something like this or had any luck with something similar?


texts: Speeding Up pgr_dijkstra using bounding box in PostGis2.0

texts:
I'm using PostGIS for the first time to calculate distances between two points, for a driving school booking system I am working on. I have loaded in all of the maps for Australia, but the query I wrote takes a good 5-6 seconds to execute. Most of the journeys are very short so I was thinking of using a bounding box but the relevant function (dijkstra_sp_delta) appears to be deprecated. Does anyone have any insight on how to implement this or speed up these kinds of functions more generally? I would like to get the query time down to 1 second or less.

SELECT SUM(
        (SELECT ST_length(ST_Transform(the_geom,28355))
          FROM ways w
          WHERE w.gid = id2
        )) as distance_in_metres
      FROM pgr_dijkstra(
        'SELECT gid AS id, source, target, ST_length(ST_Transform(the_geom,28355))/1000*60/maxspeed_forward AS cost
        FROM ways', 247201, 298103, false, false);


texts: Publish to AGOL with Editing Enabled (and other Operations and Properties) using ArcPy (ArcGIS Pro)

texts:
I have a script that accesses an APRX and publishes a feature service to AGOL. I cannot figure out how to enable editing in the script or any other of the Operations and Properties in the screenshot below (with exception of the Export Data which is a property of getWebLayerSharingDraft)

Is there a way to Enable editing, choose one of the radio buttons, Enable Sync, and check on or off the Properties such as Preserve editor tracking info using ArcPy? Essentially I want the script to perform enable the below options.

I am using ArcGIS Pro version 2.7 and the standard install of python that comes with it.


texts:
Is there a way to replicate the Split Vector Layer function from QGIS in R?

I have a Shapefile containing around 20 polygons that I'd like to split into 20 separate files.

I guess since it's in QGIS there would be a solution in Python, but I'd like to integrate it into a R workflow if possible.


texts: How to compute areas of influence in QGIS?

texts:
I am trying to create a polygon layer showing polygons that correspond to a closest store.

For a set of 30 store locations, the result should be a polygon layer with 30 features, one for each store. Each feature will represent an area where it's corresponding store is the closest. For example, an address within polygon 12 means that store 12 is the closest.

I have my store locations layer and OSM streets layer with max speed column. The resulting polygon layer should be based on the OSM streets layer and not simply linear areas.

The idea is that given a fixed set of stores and roads, the closest store should remain constant at any given point on a map. Thus, I am hoping that the resulting polygon layer will be gapless.

This video of Maptitude software provides a perfect example of what I'm trying to accomplish (fast forward to 1:55 of the video). Alternatively, see the image below for an example:

For each hospital, there is a corresponding area where anything within that area is closest to the hospital located within that area.

How can I replicate this using QGIS?


texts: Creating Maximum-Value Composite using QGIS?

texts:
I am looking to create a Maximum-Value Composite raster in QGIS

I have three raster images:
Oct 2012
Nov 2012
Dec 2012
All are in TIFF Format, Float32, EPSG:4326 - WGS 84

I want to stack these images, then on a pixel-by-pixel basis,examine each value, and retain only the highest value for that pixel location to create a Maximum-Value Composite (refer to http://en.wikipedia.org/wiki/Maximum-value_composite_procedure)

e.g.(these are made-up numbers for illustrative purposes):

Pixel in Row 205 and Column 106 will have three values:
90
75
100

I want to choose 100, then move to the next pixel and do the same. In the end I should have a final raster with maximum values only, a maximum value composite!

I am have QGIS 1.8.0 & 2.0.1, GRASS 6.4.3RC2 and SAGA 2.0.8 at my disposal.


texts:
I have population data (inhabitants per city and year) of 31 cities of the DR Congo for the years 1960, 1970, 1980 and 1990. I want to interpolate the population of all cities for the year 1990, to estimate the growth of them in order to locate the sites that will minimize the competition for resources.

I have used an IDW and a layer of existing lakes as a barrier. However, it seems to me that the resulting layer does not have sense, and so I would like to ask if IDW is the best interpolation tool for my data or it will be better to use other, such as Kriging.

Thanks!

Update:

I have to create a protected area in the DR Congo, so I want to estimate the possible growth of the cities in order to find the more suitable place to avoid competition for space and resources.


texts: Arcpy - process large shapefile in groups of 100 via cursor

texts:
I am working with a large(ish) shapefile containing 30,000 features. I am trying to run a zonal histogram and have been running into many problems. Running this function with the entire set of features crashes the tool. I've tried rasterizing the feature by FID but this also gets hung up and never finishes. I've written some python code that runs a zonal histogram on each feature but this has not really been finishing and is taking a long time.

I would now like to try running the tool by taking groups of 100 features and running the zonal histogram on each of those groups, outputting a uniquely named DBF file. I can then process all the DBF files in R to create one final table.

Below is my code for running the tool feature by feature. I found writing a temporary shapefile is more efficient than creating a temporary lyr file so I'd like to keep this approach. How could I adapt the code below to run in groups of 100 features?

import arcpy
import os
import os.path
from arcpy import env
from arcpy.sa import *
import gc

rasterLayer = "D:/pathToMyRaster/rasterLayer.tif"
outputDirectory = "D:/myOutputPath/"

newlayer = arcpy.mapping.Layer('D:/pathToShapefile/parcels.shp')

arcpy.env.parallelProcessingFactor = "100%"

with arcpy.da.SearchCursor(newlayer,['OID@','statsZone']) as cursor:
    arcpy.env.addOutputsToMap = False
    for row in cursor:
        print(row[1])
        tempFileName = outputDirectory+row[1]+".dbf"
        sql="""{0} = {1}""".format(arcpy.AddFieldDelimiters(newlayer, arcpy.Describe(
            newlayer).OIDFieldName),row[0])
        arcpy.Select_analysis(in_features=newlayer, out_feature_class=os.path.join(outputDirectory,'TempShapefile.shp'.format(row[0])),
                             where_clause=sql)
        try:
            arcpy.gp.ZonalHistogram_sa(os.path.join(outputDirectory,'TempShapefile.shp'.format(row[0])), 'statsZone', rasterLayer, outputDirectory+row[1]+".dbf", "")
            del row
            gc.collect()
        except:
            print(str(errors))
            del row
            gc.collect()


texts: Copy Multiple (not all) attributes of a polyline to multiple polylines of same feature

texts:
I have multiple fields in a feature. I need to copy multiple (5-6) attributes of a polyline to a different polyline in the same feature. To make things worse, I will have around 30 different polyline to copy the attributes to.

In this answer, I found copying of all attributes of a polyline to a different polyline. The only difference is that, I need to copy specific (multiple) attributes of a polyline to multiple polylines in the same feature.


texts: Globe draped layer updating - ArcObjects vs Desktop

texts:
I have a globe document that includes several layers that are draped on the globe's surface, but no elevation layers. If I open the document in ArcGlobe (desktop) and add an elevation layer, the appearanceof the draped layers changes - they are draped on the elevation. However, if I use this same globe document in a Visual Studio application and add the elevation layer programmatically, the appearance of the draped layers doesn't change. Does anyone know why this is the case?

Here's the code I'm using to add the elevation layer programmatically:

IGlobe theGlobe = axGlobeControl1.Globe;
 IRasterLayer mergedElevLyr = new RasterLayerClass();
 mergedElevLyr.CreateFromFilePath("D:\@ALPACAS\AlpacasCT\experimental\3D\data\MergedElev_111711.img");
 mergedElevLyr.Visible = false;
 theGlobe.AddLayerType((ILayer)mergedElevLyr, esriGlobeLayerType.esriGlobeLayerTypeElevation, true);
 theGlobe.GlobeDisplay.RefreshViewers();

Any thoughts will be appreciated!

Mark


texts: QGIS-bin.exe stops working when trying to run NN-Join

texts:
I am trying to calculate nearest neighbor using the NN-Join add-in. This has worked well until now, but after I updated my data it only partly works.

I have a point layer, and I have three other layers, two point layers and one polygon layer, and I need to calculate the distance to the nearest point in the joint layers, and in the polygon layer to the nearest polygon.

My input layer (point) is in a geographic CRS (EPSG: 4326, WGS 84), and all my join layers are in a projected CRS (EPSG: 32633, WGS 84/UTM zone 33N). This used to give me the distances from the nearest features in the join layers in meters without errors.

Now, after having updated the data, and running the same procedures over again, this works fine with two of the point layers, but on the last point layers, the program aborts ("qgis-bin.exe has stopped working"), and with the polygon join layer the distances are 0 (zero), except for a few records (37 of 3723) where I get the correct distances. In the last case I haven't been able to find any system in the manner how the procedure returns the correct distances in ca 10% of the cases.

I have tried to reinstall QGIS (version 2.18.14). I had version 2.18.x previously as well, but probably not 2.18.14. I have juggled CRS's, I have read point data in over again, but the errors seem to be consistent.

Any help out there?


texts: Black edge (NoData pixels) in Landsat 8 color correction in ERDAS

texts:
My boss asked me to mosaic and color-balance multiple Landsat 8 tiles for presentation. I tried to do that in ERDAS but it returned weird result as below:

Obviously it's the black edges or nodata pixels were counted while balancing the color. Is there a way to eliminate those annoying pixels?

The tiles were exported from the direct 16-bit USGS ".tar.gz" data files into 8-bit images. And the 8-bit images were then brought in to the Mosaic Pro session for color correction. All processes were done by ERDAS. However I also tried to export the 8-bit images through ArcGIS and then process in ERDAS Mosaic Pro but the result was the same.


texts: Seeking Python library (other than ArcPy) for geoprocessing such as buffer?

texts: Excluding ArcPy, is there any python library which can do geoprocessing, such as buffer/intersect, with shapefiles?


texts: I want to create ontology for the INSPIRE data.As the there is a large amount of data I wanted to know what technique should be used for that


texts:
I'm trying to run a map algebra expression in a vb.net add-in, and can't seem to get it to work. It's a roving window that interpolates noData areas in a DEM. The expression works in the Python window, and in a Python stand alone script (after replacing the bound raster [inputDEM], with "inputDEM"). In vb it throws a syntax error. Any thoughts?

Thanks so much,

DR

Try
                Dim demPath As String = demLayer.FilePath.Remove((demLayer.FilePath.Length - demLayer.Name.Length), demLayer.Name.Length)
                Dim inputDem As IRasterDataset = OpenRasterDataset(demPath, demLayer.Name)
                Dim mapAlgebra As IMapAlgebraOp = New RasterMapAlgebraOpClass()
                Dim environment As IRasterAnalysisEnvironment = mapAlgebra
                Dim workspaceFactory As IWorkspaceFactory = New RasterWorkspaceFactoryClass()
                environment.OutWorkspace = workspace
                mapAlgebra.BindRaster(inputDem, "inputDEM")
                Dim outputDem As IRaster = mapAlgebra.Execute("Con(IsNull([inputDEM]),FocalStatistics([inputDEM],NbrRectangle(16,16),'MEAN'),[inputDEM])")
                Dim saveDem As ISaveAs2
                saveDem = outputDem
                saveDem.SaveAs("fillGaps", workspace, "GRID")
            Catch ex As Exception
                MsgBox(ex.Message, MsgBoxStyle.Exclamation, "Error")
            End Try


texts: ArcGIS ModelBuilder: How to control data overwriting from within a model?

texts:
At 10.7.1 I would like to control the Geoprocessing Options > Overwrite the outputs of geoprocessing operations setting from within a model.

Here's why: I've created a model, stored on a shared network drive, that is being used by a constantly-changing variety of individuals.  They may need to re-run the model over time from different computers.  They've been instructed to manually check the ArcMap Geoprocessing Options Overwrite option to "on", but - being human - they often forget, and resulting runs of the model fail due to the inability to overwrite.

All data is stored as FGDB feature classes.  I'm deleting all intermediate data within the model, so that is not a problem.  However, the final output of the model must be saved between model runs.  It is this final output data that I would like to overwrite the next time the model runs, regardless of the Geoprocessing Options setting.

I've searched, but haven't found a way to control this setting from within a model.  I suppose that I could modify the model to delete any existing data as the model's first step, but this seems kinda kludgey.  My model-building skills are weak, so I'm turning to SE for assistance.  How do experts solve this situation?


texts: Changing projection on each page of data driven pages

texts:
I am working with data driven pages and a raster image of the world. I have a separate polygon shapefile of all of the UTM zones (with a polygon for zones 1-60, north and south, for a total of 120 polygons).

I have added a field to the polygon shapefile where I defined the Well-Known ID (WKID) of each feature (each zone, eg, 1N, 2S, etc). These WKID's are based on the WGS 1984 UTM zones, north and south (eg 32731). I set this field to be a "short integer" per this article: http://support.esri.com/EN/knowledgebase/techarticles/detail/40411

In the data-driven pages settings I have set the "Spatial Reference" to this WKID field.

However, when I make the data driven pages, the projections don't change from page to page, and are all still rectangles (because my coordinate system for the data frame is just set to the basic WGS 1984). I'm expecting them to use this WKID to change the projection to the WGS 1984 UTM projections for each zone, so that each page of the data driven pages shows the correct projection for that zone.

Can you think of anything I'm missing to accomplish this?


texts: Why are ArcGIS 10.2 Color Ramps missing?

texts: I'm using ArcGIS 10.2 Trial and found that there is no Color Ramps. Is it a bug or a restriction?


texts: Building polygons and the lines that represent their edges in QGIS3

texts:
I want to create polygons from intersecting lines using QGIS 3, each polygon having a unique ID attribute. In addition to the polygonal areas, I need the lines that represent their edges; each of those lines needs to have attributes that are the ID’s of the polygons that lie to their left and their right. Note that the original lines are often divided into several smaller lines, each with a different pair of ‘left’ and ‘right’ attributes.

See the illustration below to get an idea of what I’m shooting for. Ideally, the next step would be to join consecutive edges that bound the same areas (e.g. the four blue lines that bound areas 0 and 1.

I did this in ARC/INFO 25-30 years ago. I just can’t figure out how to do it in QGIS 3.


texts: Using WhiteBox tools for Raster Calculation?

texts:
I understand the philosophy of WhiteBox Tools (one command for one process) but am wondering how best to use its math tools to create a slightly more complicated function.

E.g., (A/B + C/D)

I'm creating band ratios as given on the Remote Sensing index database which sometimes require functions that are more complex than the functions of Mathematical and Statistical Analysis (e.g., divide two rasters and output a TIF, add two rasters and output a TIF, etc..).

I don't want to simply load by 10 or 15 bands into memory and work with them via numpy; WhiteBox is nice because it just loads the rasters required, when necessary.


texts: Use Python to change datasource from a feature class to a view built on the feature class

texts:
I'm trying to replace the datasource for all of the layers in my map document using Python but I keep running into an error.  Right now, all of the layers point to a feature class called GISDB.DBO.TESTFEATURECLASS and I want them to point to a view that I built using the sdetable -o createview command.  It has both the object id and geometry columns from the original feature class.  It shows up in ArcCatalog with the same icon as a feature class and I can add it as a layer in ArcMap.

  Runtime error <type 'exceptions.ValueError'>: Layer: Unexpected error

mapdoc = arcpy.mapping.MapDocument("CURRENT")
workspace_path = r'Database Connections\TestGIS.sde'
workspace_type = "SDE_WORKSPACE"
dataset_name = r'GISDB.DBO.TESTFEATURECLASS_VIEW'
for layer in arcpy.mapping.ListLayers(mapdoc):
    layer.replaceDataSource(workspace_path, workspace_type, dataset_name)

How can I determine which property is causing the error?


texts: How to recenter extents of my coordinate system in QGIS?

texts:
How do I change the extents of a projection of the world so that longitude 90W is centered?
I'm working in WGS84.
Perhaps I need to add some clarification of my objective.
I have a dataset of the world

Which I would like to reproject to look like this.

One of the packaged CRSs, EPSG:3832 shows something like what I want to do, but I can't figure out how to change the extents at all.

Let me know if further clarification will help answer this question.


texts: How can I use GRASS r.catchment in QGIS?

texts: I'd like to use the r.catchment tool from GRASS within QGIS. It looks like this is a tool that hasn't yet been developed into a processing tool under the GRASS commands. Is there some other way to load GRASS tools within QGIS?


texts: Converting MapXtreme FeatureGeometry to something writable to SQL Server geometry type

texts:
I'm new in MapXtreme and I'm bit lost.

I have polygon in MapXtreme 7.0 (in C#) saved as FeatureGeometry type. I need to save this in SqlServer 2008 as geometry type, but I don't have direct access to DB. So, I need to push data through wcf to be read by some process (in C#) on server with DB access.

What could be order of types conversions which lead from MapXtreme's FeatureGeometry through something serializable for wcf to SqlServer's Geometry?

Of course, polygons will be read from DB and pushed back to be displayed on the map. So conversions in this direction needed as well.


texts: Can ArcGIS Desktop or ArcPad use the Windows Location Provider?

texts: Are there any plans for ArcGIS Desktop or ArcMap to get the GPS Coordinates via the 'Windows Location Provider' in Windows 7/8 instead of reading a virtual com port?


texts:
I'm writing some software that allows users to draw polygons over the raster from a geotiff file.

I can calculate the area of the polygon in pixels; in fact, this calculation is a byproduct of finding the polygon's centroid (which I also report).

However, I would like to convert area in pixels to acres.

Presumably geotiffs have some metadata linking the raster to real world units. I'm vaguely familiar with the concept that internally geotiffs use a "Model (M)" and "Raster (R)" and that there is metadata relating them.

I'm trying to find exactly what metadata there is, and how to query it (using libgeotiff).

For example, if I had a value correlating pixels to square metres, it should be straight forward math to go from pixels -> metres -> acres.

If it is not possible to get an exact correspondence between pixels and real-world units, a method to calculate a reasonable approximation would suffice.


texts: Build arcpy query to select features containing one OR another one string value?

texts:
I have a table. Based of field name "TYP", I would like to select rows which contains string values "a" or "b". I need to specify the query to my arcpy code.

fc = "C:\\Users\\type_jan18.shp"
field = "TYP"
val1 = "a"
val2 = "b"

I understand how to build expression if I have one condition:

# select by one condition
whereClause = '"TYP"' + " = '" + str(val1) + "'"

# by using FieldDelimiters, again one condition
whereClause = """{} = str(val1) """.format(arcpy.AddFieldDelimiters(fc, field))

But how to correctly state it with the OR statement?

This one seems not working...

  whereClause = field + " = '" + val1 + "'"+ "or"  + field + " = '" + val2 + "'"


texts: where are the rasters that make up a layer?

texts:
I used raster calculator to add up 12 rasters made up of 0's and 1's to make a fire frequency map. Then I saved the output as a layer (lyr).

However, I need the raster output.

Did I add it up wrong, save the files wrongly or are they saved some where?


texts: Is there a zonal statistics tool to calculate median and mode raster values for each polygon?

texts:
I'd like to calculate zonal statistics for a raster where the zones are defined by polygons in a vector layer within QGIS. Specifically I'm interested in calculating the median and mode raster values for each feature. QGIS includes a 'zonal statistics' plugin, but it only seems to calculate mean and sum statistics.

Is there another plugin that could do this task for me? I've looked at the source code for the existing statistics plugin and it looks like it would be fairly easy to modify, although my skills in C++ are fairly lacking.

QgsZonalStatistics source:
http://www.qgis.org/api/qgszonalstatistics_8cpp_source.html

How to use QGIS Zonal Stats Plugin from Python Console?

Raster statistics available in ArcGIS:
http://resources.arcgis.com/en/help/main/10.1/index.html#//009z000000w7000000


texts:
I have a plugin that loads layers from SQL Server 2008. It has worked before but when I reset my PC, i wasn't able to backup my plugins. I am left with an older copy. Below is my python code to load the layers:

harvOps = "MSSQL:server=FSDSVR;database=AGDAT_FARMMAP_BIOPOWER;tables=dbo.V_HARV_OPS;trusted_connection=yes"
harvOpsLyr = QgsVectorLayer(harvOps,"Harvest States (Old)", "ogr" )
QgsMapLayerRegistry.instance().addMapLayer(harvOpsLyr)

When i run the plugin in my new QGIS, it gives me this error in the Log Messages Panel under OGR tab:

2018-07-28T14:23:24 1   Data source is invalid (Unable to initialize connection to the server for MSSQL:server=FSDSVR;database=AGDAT_FARMMAP_BIOPOWER;tables=dbo.V_FARM_CUR;trusted_connection=yes,
            [Microsoft][ODBC Driver Manager] Data source name not found and no default driver specified
            Try specifying the driver in the connection string from the list of available drivers:
            SQL Server
            SQL Server Native Client 10.0
            )

Even this old version is already working previously. I have already checked my Data Sources (both 32 and 64-bit) and they all seem to be fine.

When i try to add layers from the QGIS interface, it works fine.


texts: How to work around Field Calculator's character limit?

texts:
I'm trying to use the field calculator in ArcMap, however my code is way too long (1337 lines) and there seems to be a character limit:

I'm only allowed to paste the first 103 lines.

Is there an easy work around? I'm using ArcMap 10.1.


texts: How do I build a custom menu to change preferences for a MapBasic application?

texts:
I would like to be able to input a few different tables for an application.  I want to be able to set the path before the app runs.  How do I create a preference box and what code do I use to call the tables that I set in the box?
I want it to resemble this:


texts: Do Label Expressions Slow Down Drawing Speed?

texts:
Further to my earlier question about using different fonts on the same PC, I am wondering if label expressions would slow down drawing speed for a MXD?

Here is my label expression:

def FindLabel ( [GEONAME] ):
      import os
      if str(os.getenv("COMPUTERNAME")) == 'GIS-HP':
        return ("<FNT name='Adobe Fan Heiti Std B' size='8'>" + [GEONAME] + "</FNT>")
      else:
        return ("<FNT name='Arial' size='8'>" + [GEONAME] + "</FNT>")

This works like a charm and I haven't noticed anything yet, but I feel like adding this to all of my layers will dramatically slow down the drawing speed on my MXD. Drawing speed is one of the more important things for me as I share my MXD's with other users for them to do basic analysis. Usually I save a separate MXD, open it and change all the labels to something easy before sharing. If this doesn't slow down the drawing speed, that would be fantastic and save me some time down the line. Anyone have any ideas?


texts: Joining line features of streets and parking lots?

texts: I have a feature class (line) that include street and parking lots. How would I create it so that wherever the parkinglot touches the street - it creates a merged drawing (so it doesn't look like there is a wall there instead it flows through)? Picture of what it looks like now is below.


texts: Overlaying two raster data sets in ArcGIS for Desktop?

texts:
I am using ArcGIS 10.3

I have 2 raster datasets with different coordinate systems, but they are the same county (just different databases). When I add them to the same map layer they are far away from each other. I want to be able to match up the two datasets so I can analyze them together.

I have tried Project Raster to change the coordinate system of one of the datasets but after 2 hours of trying to execute the function I gave up (it still said it was executing, no errors popped up while I was waiting).

Is that the correct function, should I just keep waiting, or can I use georeferencing to match my counties up even if they are in different systems?

*Another issue I am encounter is that the two images have different cell sizes (98X98 vs. 1X1). Is there a way to change cell size of one of the images?


texts: I used the tool "Extract by Mask" to clip a raster. It resulted in 5 tif images: mask, maskc1, maskc2, maskc3, and maskc4. What are the maskc1-4 files? They're large and I'd like to delete them if they won't effect my normal mask tif. I'm trying to preview them in ArcCatalog but it's taking forever to build pyramids.


texts:
I am new to use Java and Geoserver.

I have a point sampling location (N,P,K Nutirients value for few fields spanning acroos few acres) on which I want to execute IDW or Natural Neighbourhood Interpolation (Vector To Raster Interpolation) in GEOSERVER and show result on WebMap.

I have done the same thing in ArcGIS using ArcGIS spatial Analysis Tool. Now I have to do the same thing in Geoserver

The Geoserver provide Barnes Surface Interpolation but it does not seems to be suitable for my case since smpling locations  are too close.

After searching on Google I understand that Interpolation can be done using GDAL library (GDAL_Grid) by Creating WPS service in geoserver.

I am unable to find proper examples/documentation on How to use GDAL libraries from within GeoServer (I am using Java Platform)

It will be great help for me if anyone can provide pre-requisites of using GDAL libraries with GeoServer and examples of the same.

Also let me know if there is any other way of achieving IDW/ Natural Neighbourhood Interpolation in GeoServer.


texts: Sharing the same file in QGIS and Google Earth

texts:
Is there some layer format that can be edited both in QGIS and Google Earth?
So for example, when I change style of some line in QGIS, it appears also in Google Earth and vice versa?


texts:
My goal is to select all of the features in a layer with certain names.

I have a list of 86 Subdivision names and a I am using a county Subdivision dataset, the names in the dataset include the subdivision names but also some text before and after the name (i,e. REPLAT, and PLAT 1, etc).

I created an expression to return the names from my list but it is displaying the following parser errors in QGIS (my query is displayed in the image as well).

Query text (sample):

SELECT *
FROM 'SUBDIVISION_POLYGON_614b57b6_a27b_4aaa_973c_aa7da31cedd6'
WHERE "PLAT_NAME1"  = 'austinbrooke acres%' OR
'%austintown park%' OR
'%beverly hill%' OR
'%bridget heights%' OR
'%broadview%' OR
'%burgess run heights%' OR
'%candywoods estates%' OR...

What do these errors mean?
Can this work or is another approach better?