texts: Trying to use a python script to manipulate a value from attribute table

texts:
I am using ArcGIS ArcMap 10.2 and I am trying to learn how to write a python script to use a value from an attribute table and perform some checks to it and assign an output value. At this moment all I want to do is multiple the input value from the attribute table by 2 and output it. This is as the below screenshot.

What the error message I get in the result is

Executing: CalculateField taranaki_dairy_farms stuff thickness PYTHON_9.3 "thickness =
!River_ID! * 2"
Start Time: Sun Jun 07 20:34:23 2015
Failed to execute. Parameters are not valid.
ERROR 000989: Python syntax error: Parsing error SyntaxError: invalid syntax (line 1)
Failed to execute (CalculateField).
Failed at Sun Jun 07 20:34:23 2015 (Elapsed Time: 0.00 seconds)


texts:
We are interested at the "Pole of inaccessibility" (https://en.wikipedia.org/wiki/Pole_of_inaccessibility) of our country.

Our approach is to use QGIS, and download coastline ([ngdc.noaa.gov]), and repeatedly create "contour lines of equal distance" until there is a point which is the farthest point away from the coastline.

Do you know how to generate a polygon which is an inner polygon and has equal distance the outer one?

An example of contour lines of equal distance on the Wikipedia page:

After a reminder from @underdark, I finally get the result. Here are my steps:

Get Mainland Coastline Shapefile

$ curl https://www.ngdc.noaa.gov/mgg/shorelines/data/gshhg/latest/gshhg-shp-2.3.5-1.zip -o gshhg-shp.zip

[action] use QGIS to add a field value with value 1 to the polygon of mainland coastline
Rasterize the Polygon of Mainland Coastline

$ gdal_rasterize -a value -tr 0.003 0.003 coastline.shp mainland.tif
Proximity the Rasterized Mainland

$ gdal_proximity.py -values 0 -distunits GEO -of GTiff mainland.tif proximity.tif
Contour Based on the Proximity Result

$ gdal_contour -a distance -i 0.005 -f "ESRI Shapefile" proximity.tif proximity_contour.shp

You now have the contour lines of inaccessibility!


texts: Copying Multiple Selections (by Attribute) to One Output Shapefile using QGIS Python scripting?

texts:
INTRODUCTION

The code below makes a selection and then saves to an 'Output' file (shapefile?) each time through the "FINAL ITERATION LOOP".  Each 'Output' file will contain all selected civic address points (from the "Input" shapefile) on one of the streets assigned to employee Joe Blow (Joe Blow's streets come from the CSV file "FSS_Streets_Test_Sample" and are appended to the "joeblow" list).

QUESTION

But what if I wanted only a single output shapefile (or layer)?  One with ALL civic address points on ALL streets in the 'joeblow' list?

Perhaps the trick is to give "selectbyattribute" multiple street names as the final parameter?  Is that possible?  If not, how do I make multiple selections and essentially concatenate them one after the other until all (z) civic address points on all streets belonging to "Joe Blow" (in the 'joeblow' list) have been highlighted before then saving them in a (single) new shapefile?

Tip - most of the script is just 'laying the ground work' such as building lists (e.g. unique list of all streets ("streets"), and streets assigned to "Joe Blow" ("joeblow").  Perhaps you'll prefer to skip to the end, where the two "qgis" tools are executed.

CODE

from PyQt4.QtCore import *
from qgis.core import *
from qgis.utils import iface
import processing

##Input=vector
##Input_Field=field Input

# Assign .CSV of streets assigned to employees, "FSS_Streets_Test_Sample", to variable 'vl'
vl = QgsMapLayerRegistry.instance().mapLayersByName('FSS_Streets_Test_Sample')[0]

# Set vl as the "Active layer"
iface.setActiveLayer(vl)

# Instantiate 'streets' list
streets = []
# Get List of Unique streets from 'STRNAME' column in 'vl'
for f in vl.getFeatures():
    if f['STRNAME'] not in streets:
            streets.append(f['STRNAME'])

# Instantiate list of streets assigned to "Joe Blow"
joeblow = []
# fill that list with streets belonging to "Joe Blow" as indicated in 'vl'
for f in vl.getFeatures():
    if f['FSS'] == 'Joe Blow':
        joeblow.append(f['STRNAME'])

#Count number of streets belonging to "Joe Blow"
z = len(joeblow)

# FINAL ITERATION LOOP: Iterate through streets belonging to "Joe Blow", each time selecting from the "Input" shapefile, all fields (civic addresses) where the "Input_Field" (STRNAME), equals ("0"), particular street(s) belonging to 'joeblow' (street(s) within the 'joeblow' list)
for k in range(0,z):
    processing.runalg('qgis:selectbyattribute', Input, Input_Field, 0, joeblow[k])
    processing.runalg('qgis:saveselectedfeatures', Input, 'Output')

# End Script

Notes:

Inputs given to GUI when script is run...

The 'Input' file is a shapefile of civic address points (along streets obviously).

The 'Input_Field' is a column called "STRNAME" within the 'Input' shapefile

The CSV file - "FSS_Streets_Test_Sample" is in the Layers Table of Contents of my Map Project.  It basically has 2 columns: employee name ("FSS"), and street name (also called "STRNAME").

Incidentally when the script finishes executing...

a) Civic addresses relating to the last street (but not all streets) in the "joeblow" list are now highlighted on the 'Input' shapefile layer.

b) The GUI disappears without throwing any error messages in the GUI's "Log" Tab (although there are messages in the various tabs of the "Log Messages" window.

c) But... I'm not seeing any 'Output' files (temp layer or shapefile) generated by the "qgis.saveselectedfeatures" tool.


texts: Calculate area of intersecting polygons in QGIS

texts:
I have two collections of shapefiles:

1 - postal codes, administrative and political areas

2 - various land use categories like parks & allotments, and points like trees

I want to produce statistics on the area of those land use categories broken down by postal code, admin & political area, and the number of points by the same means. For example, there are 4.1ha of park in SE15, 879 trees in SE22, and so on.

I am using QGIS. I've processed my raw data to get all the shapefiles ready, but after much fruitless searching around on this and other web sites I can't see a way to do this.


texts: Comparing two tables and appending only the new data

texts:
I need to compare the data in two identical tables in different copies of a geodatabase and select the rows in Table 1 that aren't included Table 2 so that they can be added on to Table 2 (part of the locally stored, main copy of the geodatabase).  I will have to do this task regularly, so I would like to automate it as part of a model (I've already built a model to do this task for the shapefiles in the geodatabase).

So far I've tried Add Join and Compare Table but it doesn't look like either will do what I need.

What are other ways to accomplish this task?


texts: Driving distance between 2 points using commercial roads

texts: How can I get Driving distance between 2 points using only roads available to commercial vehicles?


texts: Use Python in ArcGIS Field Calculator to return value based on numeric value in any of three other fields

texts:
I am trying to calculate a new field in ArcGIS 10.5 Desktop using the field calculator, based on the value in any or all of three other fields in the same table. I have not used Python for a while, and am forgetting how to do this.

In essence, this should calculate if any of three types of a crop are being grown in an area. If at least one has an area greater than 1 hectare, then it should return "1" (if two or three types are present, it should give the same answer). The input fields and criteria:

Dry season rice > 1 hectare
Wet season rice > 1 hectare
Upland rice > 1 hectare
If any of these are true, then return "1" to new field I have already created, otherwise, new field is 0.

I have tried a few options in Python, but none have worked so far, e.g.:

def calc(Dry_area_1,Wet_area_1,Up_rice_1):
    if Dry_Rice > 1:
        Grow_Rice = 1
    elif Wet_Rice > 1:
        Grow_Rice = 1
    elif Up_Rice > 1:
        Grow_Rice = 1
    else
        return 0

Am I overlooking something very simple at a late hour of the night here?


texts: Using NNJoin-Plugin through Python console

texts:
I want to write a Python script which I can execute from within QGIS 2.18 (from the Python console).
The Python script should calculate the nearest line (out of a set of lines) of a given point (out of a set of points) and I have already accomplished this task by manually using the NNJoin-Plugin. Now I want to automize this process within a python script, but I do not know how to run the NNJoin plugin within a Python script.
Is this even possible? If yes, could you give me some advice on how to accomplish this?


texts: I get a transparent/empty image with Mapnik in python from a raster file

texts:
I start by selecting a raster from the website of the Swiss Federal Office.
Therefore, I know, with certainty, the projection in which this layer is located.

Source : https://shop.swisstopo.admin.ch/en/products/height_models/alti3D
-> instructions: click on 'read more', 'sample data' and 'TIF + TFW'.

Then I write the xml file defining the rendering of my layer as below :

<?xml version='1.0' encoding='utf-8'?>
<Map  background-color='#00000000' srs="+init=epsg:4326">
  <Style name="generalMapStyle">
    <Rule>
        <RasterSymbolizer>
            <stop value = "0"/>
            <stop color="red" value = "25"/>
            <stop color="red" value = "500"/>
        </RasterSymbolizer>
    </Rule>
  </Style>
  <Layer name="Any raster" srs="+init=epsg:2056">
    <StyleName>generalMapStyle</StyleName>
    <Datasource>
      <Parameter name="file">rater_file.tif</Parameter>
      <Parameter name="type">raster</Parameter>
    </Datasource>
  </Layer>
</Map>

Finally, I launch the python program below:

import mapnik

m = mapnik.Map(256, 256)
mapnik.load_map(m, "stylesheet.xml")
m.zoom_all()
mapnik.render_to_file(m, "map_xml.png")

I obtain in return a png file containing only a transparent layer.

Can someone please help me to figure out what's going wrong?
Feel free to share others tools for map rendering in Python.


texts: Setting selection color transparent and border color red using PyQGIS

texts:
I have a layer with different informations. I want to change the selected polygon color to "Transparent". I am able to achieve that using:
iface.mapCanvas().setSelectionColor(QColor("Transparent"))

But I want a different border color such as red or any other color. For an example in ArcMap, there is an option to set the border color, border width and many more.
How can I achieve that using Python script in QGIS?


texts: ArcGIS REST Services Directory - this type of layer cannot be turned on or off

texts:
I'm new to arcgis.

I've downloaded the layer with all the geological information from here

map server

I only need to see some of the information of this map, but I get "this type of layer cannot be turned on or off".

How can I go around this problem?


texts:
I found this helpful post that got me started using ogr2ogr in python: Using ogr2ogr to convert GML to shapefile in Python? as well as this one: Converting ArcGIS Server JSON to GeoJSON? I am trying to take JSON from arcgis service query results and convert to geojson. I am querying OBJECTIDs in increments of 1,000 as the query results limit is set to 1,000. I tested this and it works in command line, but when I run the script, I get an error saying unable to open data source. Can someone help?

import os
import ogr2ogr

def main():
    for i in range(0,83000,1000):
        url = ('https://gisrevprxy.seattle.gov/ArcGIS/rest/services/SDOT_EXT/sdot_parking/MapServer/5/query?'
               'text=&geometry=&geometryType=esriGeometryPoint&inSR=&spatialRel=esriSpatialRelIntersects&'
               'relationParam=&objectIds=&where=OBJECTID%3E+'+str(i)+'&time=&returnCountOnly=false&returnIdsOnly='
               'false&returnGeometry=true&maxAllowableOffset=&outSR=&outFields=OBJECTID%2CShape%2CSEGKEY%'
               '2CELMNTKEY%2CDISTANCE%2CWIDTH%2CSIGNTYPE%2CSTANDARDTEXT%2CCATEGORY&f=pjson')

        outfile = os.path.join('/Users/kirsten_lawrence/Documents','recs_'+str(i)+'up.json')
        # note: main is expecting sys.argv, where the first argument is the script name
        # so, the argument indices in the array need to be offset by 1
        print 'creating geojson starting at record',i
        ogr2ogr.main(['','-f', 'GeoJSON', outfile, url, 'OGRGeoJSON'])


texts: QGIS "invents" coordinate decimal places from imported PostGIS geom

texts:
I have a PostGIS database with an geom-accurancy of 9 decimal places for X and 8 decimal places for Y.
Reading out the coordinates of an PostGIS-Polygon in PyQGIS shows:
import psycopg2
conn = psycopg2.connect('spb', 'postgres', 'postgres', 'localhost', 5433)
cur = conn.cursor()
sql = f"""SELECT ST_AsText(geom)
          FROM geom
          WHERE id = 19695"""
cur.execute(sql)
cur.fetchone()

# Output:
('MULTIPOLYGON(((349749.532007553 5575769.14298604,349777.370999882 5575792.32299964,349792.406197559 5575794.8902163,...)',)

But when I import the polygon from PostGIS as a layer in QGIS, and read out its coordinates afterwards, I see that QGIS "invents" new decimal places, 17 for X and Y.
uri = QgsDataSourceUri()
uri.setConnection("localhost", "5433", "spb", "postgres", "postgres")
uri.setDataSource ("public", "geom", "geom", "id = 19695")
vlayer = QgsVectorLayer(uri .uri(False), "geom", "postgres")
QgsProject.instance().addMapLayer(vlayer)

lyr = iface.activeLayer()
features = lyr.getFeatures()
for f in features:
    f.geometry()

# Output
<QgsGeometry: MultiPolygon (((349749.53200755297439173 5575769.14298603963106871, 349777.37099988199770451 5575792.32299964036792517, 349792.40619755897205323 5575794.89021630026400089,...>

So what QGIS is doing here? How can I prevent that QGIS "invents" non existing decimal places?


texts: ESRI ArcGIS for Server AMI for Amazon EC2 reserved instances

texts:
I read several online documents on using ESRI's Cloud Builder to optmimally set and fire up ArcServer AMIs on Amazon EC2. However, everything seems oriented toward selecting and launching on-demand EC2 instances, not reserved instances.

Does anyone know how to get started with ArcServer AMIs on reserved EC2 instances? Does Cloud Builder not allow this type of set up?


texts:
Im looking to create a map that can be interactive.

The best option that I found is leaflet, the thing is I don't find any resource explaining how to create my own map.

Im looking to create an indoor map where user can see all the building

How could I achieve that ?


texts:
All captions are displayed in native language when using Openstreetmap as source and that makes map too hard to read when looking at Asian countries.
Is it possible to use user selected language or English transliteration for all map texts?


texts: How to automatically import field value of polygon to point field based on spatial query?

texts:
I have two layers, points and polygons. Point layer is empty. When I add new point feature I want to update it's field by values of polygon field if new point is within polygon. I need solution in QGIS.

How can I do this?


texts: How to show elevation on mouse hover in ArcMap?

texts:
I need to show elevation of a DEM on the bottom bar (or wherever) on mouse hover in ArcMap 10.2.2. Is it possible?

I know I can access to it by opening "Identity", but that is not practical. It would be better to have information of it without clicking.


texts:
I'm here because I'm trying to create a GTA V map with Leafletjs API. In deed, I've got my tiles that are ready. But the aim is to put markers on the map by using GTA V in-game coordinates. To achieve it, I need to convert the coordinates to LatLong but after many researches, I can't find a solution.

At Stackoverflow, they advise me to ask my question here.

For now, I have a simple function that help me to try to "calibrate" my coordinates :

function GTAVConverter(x,y)
{
    var ex = 0.01255;
    var ey = -0.01260;

    var x = ex * x;
    var y = ey * y;

    var x = x - 486.97;
    var y = y + 408.9;

    return map.containerPointToLatLng(new L.Point(x+1440.45,y+127.45));
}

But I think there is a better way to achieve my aim. This function isn't really accurate, and if the browser window is resized, all the coordinates are moved.

Example that I want to convert :

I know that Leaflet API has a coordinates system but I don't know how it works.


texts: Converting Layout/Bookmark python code to script tool?

texts:
I have a python script that will position the scale bar, north arrow, title, and legend on a map layout in ArcGIS Pro. This script also references a bookmark and then exports the layout to a PDF. I need to turn this script into a script tool so that it can be used in ModelBuilder, but I have no idea where to start.

Here is the python script I need to convert to a script tool in ArcGIS Pro.

# input name of layout

    p = arcpy.mp.ArcGISProject("CURRENT")
    lyt = p.listLayouts("Test")[0]

# Reposition the scale bar

    scaleBar = lyt.listElements("MAPSURROUND_ELEMENT", "Scale Bar")[0]
    mf = scaleBar.mapFrame
    scaleBar.elementPositionX = mf.elementPositionX + 0.0
    scaleBar.elementPositionY = mf.elementPositionY - 0.5

# Reposition the north arrow

    northArrow = lyt.listElements("MAPSURROUND_ELEMENT", "North Arrow")[0]
    mf = northArrow.mapFrame
    northArrow.elementPositionX = mf.elementPositionX + 8.8
    northArrow.elementPositionY = mf.elementPositionY + 0.7

# Align the title with the center of the map frame

    title = lyt.listElements("TEXT_ELEMENT","Name of Map Text")[0]
    mf = lyt.listElements('MAPFRAME_ELEMENT',"Map Frame")[0]
    title.elementPositionX = mf.elementPositionX + (mf.elementWidth / 3.7)
    title.elementPositionY = mf.elementPositionY + (mf.elementHeight / 0.98)

# Reposition the Legend and fix legend title

    legend = lyt.listElements("LEGEND_ELEMENT", "Legend")[0]
    legend.title = "Legend"
    legend.elementPositionX = mf.elementPositionX + 7.7
    legend.elementPositionY = mf.elementPositionY + 7.15

# setting layout to bookmark

    aprx = arcpy.mp.ArcGISProject("Current")

# add name of layout

    lyt = aprx.listLayouts("Test")[0]
    mf = lyt.listElements("MAPFRAME_ELEMENT")[0]

# add name of bookmark and export as PDF

    bkmks = mf.map.listBookmarks("Castro_Py")
    for bkmk in bkmks:
      mf.zoomToBookmark(bkmk)
    lyt.exportToPDF(r"C:\arcGIS_Shared\Exports" + "\\" + bkmk.name + ".pdf")


texts: OpenStreetMap Python API Bounding Box Request

texts:
I am new to OSM.

When I make an API bounding box call like below:

import overpass

# bounding box coordinates
min_long = min(longs)
min_lat = min(lats)
max_long = max(longs)
max_lat = max(lats)

# query overpass within OpenStreetMap
api = overpass.API()
map_query = overpass.MapQuery(min_lat,min_long,max_lat,max_long)
response = api.Get(map_query)

The output is not wrapped.  Here is a sample when I print response:

<member type="way" ref="50345468" role=""/>
<member type="way" ref="119882166" role=""/>
<member type="way" ref="7664653" role=""/>
<member type="way" ref="7661912" role=""/>
<member type="way" ref="7662660" role=""/>
<member type="way" ref="6983960" role=""/>
<member type="way" ref="6975837" role=""/>
<member type="way" ref="6982355" role=""/>
<member type="way" ref="6983849" role=""/>

Anybody have an idea of what is going on here?


texts: Calculate stop times/distances of transport lines - tools and procedure?

texts:
I need to calculate stop times of bus lines (extrapolation)

I have:

Shapefile with all stops
Shapefile with control points where bus time was taken
Shapefile with all bus routes
CSV file with times in first stops on the route, final stops and
control points where time was taken.

All of the Shapefiles have an attribute "route number" which identifies to which transport route they refer.

I need to find a tool (and procedure) that helps me extrapolate times in each of the stop on the route, using the departure time, arrival time and times at control points on the route (which are not stops). I suppose I actually need to calculate distances first and then times based on those distances.

I am unix user so OS tools would be preferable

Any suggestion how to crack this? Attributes described below.

Stops shapefile

Description: Locations of bus stops

Attributes:

Cod_ubic_p -    stop location
Code_variant -  bus line (referenced in *Bus* shp file)
Ordinal -   order corresponding to the stop in the journey of the variant (bus line)

control_points shapefile

Description: Locations of control points

Attributes:

Cod_ubic_c -    control point code
Code_variant -  bus line type (referenced in Bus_lines)
Ordinal  -  order of given control point in the path of the variant

(bus line)
    Desc_locat -    Description of Control Point location

Bus_lines shapefile

Description: Routes of bus lines

Attributes:

Cod_line -      Transportation line code
Desc_line -     Description of the transport line (eg 112, D16, etc)
Ordinal -   Correlative number of the sub-line of the line
Code_subline -  Code of the transportation sub-line
Subline Desc -      Description of the sub-line
Cod_variant -   Code of the route (to link with *Stops* shp)
Desc_variant -      Type (Variant) description (A = Ascending, D = descending)

csv sheet Timetables

Description: Timetables

Attributes:

COD_LINE -  Transportation line code
CODE_DAY_TYPE -     Code for the day of the week
LEAVES_FROM -   Origin of the bus route
DEPARTURE -     Depature time
DESTINATION -   Destination of bus route
ARRIVAL -       Arrival time
DESC_LOCAT -    Description of Control Point (Stop) location
COD_UBIC_C -    control point code (referenced in *control_points* shp file)
TIME    -   Time at control point location
ORDINAL -   Correlative number of the sub-line of the line
DESC_VARIANT -  Type (Variant) description (A = Ascending, D = descending)


texts: Leaflet side by side plugin is not working with imageOverlay

texts:
I'm trying to use the side by side plugin. I already check in the issues and change the getContainter() with getPane() in the .js file.

this._divider.style.left = dividerX + 'px'
    this.fire('dividermove', {x: dividerX})
    var clipLeft = 'rect(' + [nw.y, clipX, se.y, nw.x].join('px,') + 'px)'
    var clipRight = 'rect(' + [nw.y, se.x, se.y, clipX].join('px,') + 'px)'
    if (this._leftLayer) {
      this._leftLayer.getPane().style.clip = clipLeft
    }
    if (this._rightLayer) {
      this._rightLayer.getPane().style.clip = clipRight
    }
  },

I'm using imageOverlay but i don't know where to add the map.createPane('left');map.createPane('right'); that is mentioned in the issue.

var map = L.map('map', {
    center: [23.140, -101.887],
    zoom: 5
});

var catedralVieja = 'http://www.lib.utexas.edu/maps/historical/newark_nj_1922.jpg',
    imageBounds = [[7.9409, -131.1589], [29.2144, -82.6558]];
L.imageOverlay(catedralVieja, imageBounds, {opacity:0.7, pane:"left"}).addTo(map);

var catedralNueva = 'https://img.vixdata.io/pd/jpg-large/es/sites/default/files/a/atlantis-ciudad-perdida.jpg',
    imageBounds = [[7.9409, -131.1589], [29.2144, -82.6558]];
L.imageOverlay(catedralNueva, imageBounds, {opacity: 0.7, pane:"right"}).addTo(map);

var hash = new L.Hash(map);

L.control.sideBySide(catedralVieja, catedralNueva).addTo(map);

This is the jsfiddle that i have and not even the divider appears :(


texts: Loading GeoJSON from local server to show by OpenLayers

texts:
I've written some code to create a vector layer based on a GeoJSON file located on my computer. I've deployed the tomcat server and put my GeoJSON file in some path like:

  webapps/OpenLayersAsset/countries.geojson

Next using the following code to show the file but I get the error "Cross-Origin Request Blocked: The Same Origin Policy disallows reading the remote resource at http://localhost:8080/OpenLayersAssets/countries.geojson". Does it cause the byGeoJSON location or something else?

var countrySource = new ol.source.Vector({
    format: new ol.format.GeoJSON(),
    url: 'http://localhost:8080/OpenLayersAssets/countries.geojson',
    projection: 'EPSG:3857'
});
var countryLayer = new ol.layer.Vector({
    source: countrySource
});

var center = ol.proj.transform([51.34, 35.65], 'EPSG:4326', 'EPSG:3857');
var view = new ol.View({
    center: center,
    zoom: 8
});

var map = new ol.Map({
    target: 'map',
    layers: [countryLayer],
    view: view
});


texts: Finding holes in proposed districts (Polygons/Multipolygons)

texts:
Thinking through something that I seek some SQL on:

One of the ideas in redistricting is that a proposed district boundary should have no holes in it, i.e., a chunk of land with no district_id attached to it, where each chunk is made up of a tabulation_block.   Each row in the proposed_districts table would have a pid, district_id, geom , with the geometry created to represent the UNION of every tabulation_block (separate table) assigned a district_id.

There are two things that I am looking for in a SQL query:

Tabulation blocks whose geometries in the district do not have a district value assigned
Naturally existing interior rings.

I have appended the DESCRIBE TABLE results for the two tables, tab blocks and a district. I suspect that the latter involves something involving ST_DUMPRINGS(), but I cannot quite put it together.

Intended result: Table with each row

primary_key, district_id, geom_of_hole

   Column   |            Type             | Collation | Nullable |                      Default
------------+-----------------------------+-----------+----------+---------------------------------------------------
 gid        | integer                     |           | not null | nextval('rapis_five_districts_gid_seq'::regclass)
 district   | character varying(50)       |           |          |
 population | integer                     |           |          |
 geom       | geometry(MultiPolygon,4269) |           |          |
Indexes:
    "rapis_five_districts_pkey" PRIMARY KEY, btree (gid)
    "rapis_five_districts_geom_idx" gist (geom)```

```                                         Table "public.tl_2019_06_tabblock10"
   Column   |            Type             | Collation | Nullable |                      Default
------------+-----------------------------+-----------+----------+----------------------------------------------------
 gid        | integer                     |           | not null | nextval('tl_2019_06_tabblock10_gid_seq'::regclass)
 statefp10  | character varying(2)        |           |          |
 countyfp10 | character varying(3)        |           |          |
 tractce10  | character varying(6)        |           |          |
 blockce10  | character varying(4)        |           |          |
 geoid10    | character varying(15)       |           |          |
 name10     | character varying(10)       |           |          |
 mtfcc10    | character varying(5)        |           |          |
 ur10       | character varying(1)        |           |          |
 uace10     | character varying(5)        |           |          |
 uatype     | character varying(1)        |           |          |
 funcstat10 | character varying(1)        |           |          |
 aland10    | double precision            |           |          |
 awater10   | double precision            |           |          |
 intptlat10 | character varying(11)       |           |          |
 intptlon10 | character varying(12)       |           |          |
 geom       | geometry(MultiPolygon,4269) |           |          |
Indexes:
    "tl_2019_06_tabblock10_pkey" PRIMARY KEY, btree (gid)
    "tl_2019_06_tabblock10_geom_idx" gist (geom)```

  [1]: https://postgis.net/docs/ST_DumpRings.html


texts: I have a shapefile file of the "cadastre" of Spain, I use the following command in Linux to get the sql file to next insert in database with postgis (shp2pgsql -s 900913 PARCELA.SHP> PARCELA.sql), but next to the insert, the polygons appear in Africa and not in Spain, he Observed that the projection is ETRS89 / UTM zone 30N, but I work with projection 900913, how can I transform the shapefile projection to 900913??


texts: Conditional chart formatting for featureCollection in GEE (e.g. green for positive values and red for negative values)

texts:
I am trying to figure out a way to format a line chart drawn from a feature collection depending on the plotted value. For example, I want to plot the featureCollection below using green for positive values and red for the negative values.
var f1 = ee.Feature(null, {var1: 1})
var f2 = ee.Feature(null, {var1: -2})
var f3 = ee.Feature(null, {var1: 3})
var f4 = ee.Feature(null, {var1: 4})
var f5 = ee.Feature(null, {var1: -2})
var f6 = ee.Feature(null, {var1: -1})

var features = ee.FeatureCollection([f1, f2, f3, f4, f5, f6]);

var chart = ui.Chart.feature
  .byFeature({
    features: features,
  })
.setChartType('LineChart');
print(chart);

The default is the plot below:
I would like it to be plotted in this way (edited on paint):

I tried adding a 'positive' property that is null when negative and the value when positive, and did the same with a 'negative' property but this creates gaps in the plot.
var categorizeFeatures = function(feature) {

  var value = feature.get('var1');

  return ee.Algorithms.If(ee.Number(value).gt(0),
    feature.set({pos: ee.Number(value),neg:null}),
    feature.set({pos: null,neg:ee.Number(value)}));

};

var features2 = features.map(categorizeFeatures);
print(features2);

var chart2 = ui.Chart.feature
  .byFeature({
    features: features2.select(['pos','neg']),
  }).setOptions({
    colors: ['red','green']
  })
.setChartType('LineChart');
print(chart2);

I'm not even sure this is possible, but I figured this forum would confirm or refute whether it is indeed possible.
The link to the code is here.


texts: convert geoJson data to interactive map tiles containing icons,labels and popups

texts:
In my search for a solution for showing up a tremendous amount of points data on a map, I ended up with an idea of generating GeoJSON data, then convert GeoJSON to map tiles that can be served using a tile server.
I saw many sites doing the same but I don't know how.

My questions are :

What tools that can be used on server to generate tiles on fly from GeoJSON data?
The tile must be interactive, clickable and contain popup messages and titles for animated points.
Can any expert here give me an idea of how to do this and what is the right way?


texts: PostgreSql query problem with aggregate functions

texts:
I need help with a PostgreSql query. I need to find all parcels which have bigger area then average. I'm using this query:

SELECT OBJECTID
FROM parcela
WHERE ST_Area(geom)<(SELECT Sum(ST_Area(geom)/count(*)))
GROUP BY OBJECTID;

but it gives me error:

  aggregate function calls cannot be nested

I tried few different ways to solve this (using 'having' instead 'where' etc.) but always some new error occurs. Can anybody please tell me how to write right query to get this information. Thanks


texts:
I have the polygon of my catchment of interest, I call it catchment A (Reyran in the schema). (format *.shp).

I have also the x-y coordinate of a dam that is situated in the catchment. I have also a *.tif file that includes elevation information of the catchment A.

I would like to define and split in GRASS (QGIS) the corresponding upstream (catchment I) and downstream catchment (catchment O) that we could define in catchment A at the upstream and downstream of the dam.

I need to get the shapefile (*.shp) of the upstream catchment. How could I do it? Here is the schema of the catchments.


texts:
I'm trying to hide 3 required fields from a layer:

#Making invisible the fields we're not interested in
desc = arcpy.Describe('fl_target_output_1')
field_info = desc.fieldInfo
# List of fields to hide
fieldsToHide = ['OBJECTID','Shape','Shape_Length']
for i in range(0, field_info.count):
    if field_info.getFieldName(i) in fieldsToHide:
        field_info.setVisible(i, 'HIDDEN')
arcpy.MakeFeatureLayer_management('fl_target_output_1', 'fl_target_output_2', "", "", field_info)
fc_target_output_2 = r'D:\ArcGIS_Pro_Projects\Optimity_Alf\Optimity_Temp.gdb\target_output_2'
arcpy.CopyFeatures_management ('fl_target_output_2', fc_target_output_2)

The above code does not throw any error but it doesn't do anything at all.


texts:
What are the recommended ways of simplifying geometries?  Keeping in mind projections, and simplifying state geometries for example.

I've heard about converting to an 'equi-distant' projection that allows simplification without distortion, and then converting back to your chosen projection.


texts: R: Krovak projection in R - identical in ArcGIS, different in R?

texts:
I have two files: shapefile and raster, both in the same coordinate system in ArcGIS: EPSG 5514.

Check information on projection in ArcGIS for raster:

S-JTSK_Krovak_East_North
WKID: 5514 Authority: EPSG

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

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

Check information on projection in ArcGIS for Shapefile:

S-JTSK_Krovak_East_North
WKID: 5514 Authority: EPSG

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

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

in ArcGIS totally identical coordinate systems.

If I import them in R by

library(sp)
library(rgdal)
library(rgeos)
library(raster)

# import raster
r1<-raster("F:/.../r1_2013.img")

#import shp
d.shp <- readOGR(dsn = "F:/.../dummy_data_arcgis_r",
                     layer = "np_by_manag")

and I want to check data projection, there is difference in projection in my PREVIOUSLY IDENTICAL CRS:

> proj4string(r1)    # RASTER
[1] "+proj=krovak +lat_0=24.83333333333333 +lon_0=49.5 +k=0.9999 +x_0=0 +y_0=0 +ellps=bessel +towgs84=0,0,0,0,0,0,0 +units=m +no_defs"

> proj4string(d.shp)   # SHAPEFILE
[1] "+proj=krovak +lat_0=49.5 +lon_0=24.83333333333333 +k=0.9999 +x_0=0 +y_0=0 +ellps=bessel +units=m +no_defs"
>

which cause problems in any spatial operation I want to perform further.

Why R modify my projection? Is it caused by various nature of my files: shapefile / raster? How can I get rid of it and which one of Krovak CRS is correct?

My aim, after finishing this transformations, is to:

convert raster r1 to polygon p.r1,
intersect polygon p.r1 with d.shp
calculate areas of each polygon

# convert raster to polygon
p.r1<-rasterToPolygons(r1, dissolve = T)

# intersect polygon from raster with another one
int.r <-intersect(p.r1,d.shp)

# calculate areas (in m2) for every polygon in intersected shp
int.r$area <- sapply(int.r@polygons, function(x) x@Polygons[[1]]@area)

My dummy files are attached here: http://ulozto.sk/xKoewAzg/dummy-data-arcgis-r-zip


texts: If you are writing a ModelBuilder model with a lot of geospatial data in it, and if one of the datasets fails or doesn't work, does the model stop and thus not produce any output?


texts:
Is there any way to use analysis from Sextante with PostGIS Geometry Type?

Edit:

I am trying to create a webapp that make some analysis (like interpolation for example) and returns the results. So I don't want to use a desktop solution.

Python could be a solution... But how?


texts:
I am getting the following error:

The log is not providing lot of information.  What is meant by data type? I can't figure out where the problem is situated.

The layers I want to merge are all line-layers produced by the algorithm 'Polygons to lines'.


texts:
Too bad the plugin is 'probably currently in used' because it worked fine all last week and is a great productivity booster when it does work.

Does anybody have an idea why this happens?


texts: Select by attributes within the fiona Python module

texts:
Is it possible to select polygons from a shapefile based upon their attributes using the Fiona Python module? I can't seem to find anything in the docs, but it seems like a strange thing not to be included in the module.

All I can think of at the moment is iterating through all the polygons, doing the calculations I want to do, and then storing the data in a database, or some other data structure, and then running my SELECT query over that.

Alternatively if there is another Python module that will do this easily then that would be great to know about.


texts: How to determine building types and add them to the attribute table

texts:
I have a city in which the buildings are defined as polygons. I have determined the area based on these polygons, however I would also like to determine the building type (i.e. terraced, semi-detached or detached building) and add this to the attribute table (i.e. link the building type to the respective building ID). Can I also do this based on the polygons?

I assume it would be possible based on a spatial analysis/query, however I just can not find a tool or a way to do this. As there are a huge amount of buildings it would take a lot of time to visually check the building type and adding this manually in the attribute table.

Could somebody help me out?


texts: What is the best map provider solution for repetitive routing use?

texts:
[edit]

What is the best solution in terms of cost effectiveness, stability & speed for running a lot of routing requests and comparing them. lots of times over and over?
I understand now that it can be a combination of several tools.

A simple solution is always better...

[original]

I want to take a route I've gotten via a routing service like Google-maps or CloudMade and add several way-points to it. then check how they affect the time and distance of the route. the problem is that i need to do this Many, many times per day. (i have big list of routes and way-points) and this can get expensive. the services are not free..

a) What is the best map / routing provider for such a thing?

b) Alternatively, is it possible to download an entire database of such maps & a routing solution to use it in-house? (only if it's cost effective)

Thanks!


texts: Problem adding add-apt-repository ppa:ubuntugis/ppa

texts:
Since upgrading my Ubuntu to focal (Ubuntu 20.04.1 LTS) I get the following error every time that I update/upgrade my system:
Err:8 http://ppa.launchpad.net/ubuntugis/ppa/ubuntu focal Release
  404  Not Found [IP: 91.189.95.83 80]
Reading package lists... Done
E: The repository 'http://ppa.launchpad.net/ubuntugis/ppa/ubuntu focal Release' does not have a Release file.
N: Updating from such a repository can't be done securely, and is therefore disabled by default.

Is there now a different PPA that I should be using for my GIS updates?  Alternatively, is there another way that I can automate GIS updates?


texts: How is arcpy.Clip_analysis() so fast against non-arcpy implementations?

texts:
I'm using a version of Solution 2 of the selected answer here to build a new shapefile of all points that fall within a given polygon object (i.e. a clip).

When run against 27K points, it takes about 9 minutes. In arcpy it takes < 10 seconds.

I cant find any documentation on the Clip_analysis() code (is it proprietary?) so I am curious as to how it can execute the checks of 27K so fast?

Is it multi-processed? Is it calling some type of map() on the iterable of points?


texts: Grid-like artifact visible in contour layers?

texts:
I am creating contours from an Aster DEM in ArcMap 10.1

I see a grid pattern appearing, particularly in areas of dense contours:

I see the same artifact if I create a Hillshade layer so presumably it is a problem in the DEM.

The grid artifact doesn't correspond to the sample size of the DEM.

I am using the original DEM - no resampling or focal statistics.

The DEM has not been reprojected (displayed in the WGS84 UTM projection it was delivered in).

Any ideas about the cause of the artifact?

If so could anyone suggest a workflow (other than resampling/ focal statistics) for removing it?


texts:
I'm new to Python. I want to create GIS app with Python using QGIS library on Window 7 64 bit.

The way I do:

install QGIS 1.7.0

I set PYTHONPATH=D:\Program Files (x86)\Quantum GIS Wroclaw\apps\qgis\python
I run python from cmd in D:\Program Files (x86)\Quantum GIS Wroclaw\bin
I type import qgis.core, a pop up window appear with text

  The program can't start because qgis_core.dll is missing from your
  computer. Try reinstalling the program to fix this problem

and text in cmd is

  ImportError: DLL load failed: The specified module could not be found.

Can anyone help me?


texts: Dissolve output not dissolved in ArcMap?

texts:
Even though polygons of an input feature have the same attributes/dissolve fields, the output does look like in the picture below.

Does anybody know what the problem is?

This problem occurs only in ModelBuilder. When I use the same data and same default parameters and do not calculate it in ModelBuilder, it works normally.

The original input feature looks like in the picture below. So the "cross division" in the picture one does not come from the original dataset but it is created by the Dissolve tool.

I tried to delete whitespaces - still the same problem.
There are no gaps between common edges.
I tried to run dissolve multiple times in ModelBuilder - still the same problem.
When I try it normally (without ModelBuilder, the same parameters), it works normally.

Any other ideas?

Since I work with many input features, I need to work with ModelBuilder.

This is what the model looks like. (Model works fine till the last part - dissolve. Goal is to generate polygons that represent areas of municipalites with the same closest facilities - in my case central cities with specific types of public institutions, also one central city can have more than one facilities)


texts: Google Reverse Geocoding API doesn't work on mobile internet

texts:
I was able to get Google Geolocations and Geocoder API to work via desktop and on my phone when its connected to wifi, but when I disconnect from the wifi, geocoder doesn't display the results.
just to explain a bit more about what my code does, i use the geolocation API to get lat/long values and then i use that output that is stored in a variable and i use it to get the city information using the reverse geocoding google API.
my code below is as follows:
jQuery(function() {
      var GOOGLE_API_KEY = 'MY_KEY_HERE';
     var outlatlong = null;
      jQuery.ajax({
        url: 'https://www.googleapis.com/geolocation/v1/geolocate?key='+GOOGLE_API_KEY,
        data: JSON.stringify({ "considerIp": "true" }),
        type: 'POST',
        contentType: 'application/json',
        success: function(data) {

    outlatlong = [data.location.lat, data.location.lng].join(', ');

          if(data.location) {
            var out = [data.location.lat, data.location.lng].join(', ');
            jQuery('#gresponse').html(out);

    var urlTwo = "https://maps.googleapis.com/maps/api/geocode/json?latlng="+out+"&key=MY_KEY_HERE";
    jQuery.get(urlTwo).success(function(data) {
       var loc1 = data.results[0];
       var county, city;
         jQuery.each(loc1, function(k1,v1) {
            if (k1 == "address_components") {
               for (var i = 0; i < v1.length; i++) {
                  for (k2 in v1[i]) {
                     if (k2 == "types") {
                        var types = v1[i][k2];
                        if (types[0] =="sublocality_level_1") {
                            county = v1[i].long_name;
                            //alert ("county: " + county);
                        }
                        if (types[0] =="locality") {
                           city = v1[i].long_name;
                           alert ("city: " + city);
                       }
                     }
                  }
               }
            }

         });
         jQuery('#lattie').html(city);

    });

          } else {
            jQuery('#lattie').html('didnt work');
          }
        },
        error: function(e) {
          console.log('error!', JSON.parse(e.responseText));
        }
      });
        return outlatlong
    });


texts: Is there a United States waterbody index number available?

texts:
In searching for this, I see that in some countries they have a waterbody index number, a unique ID for each water body.

Does such a system exist in the US and is it available online or from any government agencies? My online searches have yielded nothing.


texts: Manual range selection for graduated style vector layer in QGIS?

texts: I'm working with QGIS 2.12.1-Lyon in Ubuntu 14.04 LTS. I really like the features of this nice program but I'm badly missing one thing. Quite regular I'm working with differences. For displaying I use a graduated polar style with white for zero and red/blue for the ranges in +/- direction. This works pretty well with raster data, where I can manually set the min/max values. Unfortunately as you can see on the screenshot, setting this range is not possible for vector data (or did I just not find it?). This has the consequence that if data is not equally distributed around 0, the colours can't be interpreted that easily. Is there another way to manually set the colour range for vector layers? (I don't want to manually edit each single class as I'm loading such data quite often and this could get very annoying).


texts: PDAL error: "Couldn't create filter stage of type 'filters.cpd'"

texts:
I would like to use the cpd filter to transform a point cloud, but I am heading to the following error:

  RuntimeError: Couldn't create filter stage of type 'filters.cpd'

which comes from the validation pipeline routine.

I installed pdal via conda on a specific virtual environment. My json pipeline is currently a straight copy/paste of the example available from the documentation:

import pdal, json

pip_cpd_est_trans = json.dumps(
{
"pipeline": [
    "data/Middal_2017_06_14_clean.las",
    "data/Middal_2017_07_12_clean.las",
    {
        "type": "filters.cpd",
        "method": "rigid"
    },
    "output_cpd.las"
]
})
pipeline = pdal.Pipeline(pip_cpd_est_trans)
pipeline.validate()
pipeline.execute()


texts: Query Feature Class from Oracle table with no geometry column

texts:
Using ArcMap 10.4.1 I am trying to create a query layer from an Oracle table (Oracle info at bottom). The table has no geometry column, but has x and y coordinates as double and I am trying to create the geometry on the fly by returning a geometry field.

In the 'Edit Query Layer' wizard I have the following query:

SELECT MSLINK, sdo_geometry(2001,6469, sdo_point_type(CAE_DBA.NETDVD_OPEN_MAINS.FROM_X_COORD, CAE_DBA.NETDVD_OPEN_MAINS.FROM_Y_COORD, NULL), NULL,NULL) as ShapeField FROM CAE_DBA.NETDVD_OPEN_MAINS

MSLINK is my unique identifier (Long Integer, Not nullable)
ArcGIS is able to discover the spatial properties for the layer and recognizes the SQL as valid. When I try to open the Attribute Table I get the following error:

"Error reading OID from Table. Reading rows has been stopped. Check that the datasource is valid. Underlying DBMS error[ORA-13226: interface not supported without a spatial index".

Do I need a spatial index on the aliased Field "ShapeField"? If so, any ideas about how to do this for a query layer?

Oracle version info:
Oracle Database 12c Enterprise Edition Release 12.1.0.2.0 - 64bit Production


texts: Create map(s) to help determine site selection

texts:
I have the task of creating a map in which the objective would be to help determine a site for a conference. I have the following data:

Addresses of possible conference attendees (I have already geocoded them).
6 possible city sites in the Northeast.

I have the Spatial Analyst extension but not network analyst.

I think a major factor in determining a site would be proximity to the city site and keeping driving time to a minimum for as many people as possible, say less than 3 hours or so or 150 miles).

How can I make a heat map as a general overview of clusters? Answered below I used the kernal density analysis.

How can I figure out the best site based on least driving distance without using network analysis? Answered below.


texts: I have been having a lot of issues trying to do a heatmap.  At first the data would disappear when I created it, so I figured out that I had to reproject the CRS.  I reprojected it, ran the heatmap, and there was no output.  The screenshot is how the "Coyote_Heat2" heatmap appears in the menu, but as you can see it hasn't created any heat map.  What else could be going wrong?  There were no error messages this time.


texts: Creating ArcCatalog-friendly xml metadata file from url metadata?

texts:
How do you take metadata information from a url like this and create an xml file that ArcCatalog's "import xml" metadata tool will accept?

I have tried a few variations on pasting the url information into a text editor and saving as xml but haven't found the right method.

I am using ArcGIS Desktop 10.


texts: Choosing observer points for Viewshed tool in ArcMap to identify visibility of landform

texts:
I am using ArcMap 10.8.
I want to create a polygon that shows all the areas within a region from which an observer could see a specific mountain.
I created a polygon around the mountain, and I'm going to create random points within the polygon to use for the observer points in the viewshed tool. The polygon is about 690 square kilometers, and I want to minimize the number of points I use for the viewshed so it doesn't crash my computer.
Is there a way to estimate how many observer points I should use based on the size of the landform polygon to get an accurate representation of visibility of the whole mountain?
Edit: I have to do this for a few features, including a single peak (one point on top should work), a basalt spire with a flat top (one point on top is insufficient), a ridgeline, and a small mountain range. For the ridgeline and small mountain range, maybe a point on each peak or high point would be enough?


texts: Printing selected Region of map with several layers in QGIS?

texts:
I creaetd a map with 6 vector layers (each is a .shp-file), the geometry type varies: 4 are polygons, 2  are lines.

I would like to print a region - to be specific approx. 20 individual regions -  of the map.

Until now, I didn't find a proper way to do it. This is my workaround so far:
One vector layer (polygon) has the outlines of the 20 regions I want to print. I color the 19 regions I don't need white, and leave the one I need transparent. Hence, I kind of have the region I want and can print it with the print composer.
The drawback: I would like to print Streetnames, which are included in one vector layer (geometry type: lines; "roads" from OSM). If I turn on labels for the street names, they are not covered up by the 19 "white" regions, hence I have street names in the "white" regions that should be invisble for printing.

Do you have an idea how I can get rid of the names or - even better - a suggestion how to reach my desired result in a better way?


texts: Loss of quality of aerial photo when removing black boarder of georeferenced tiff

texts:
Can someone assist? When georeferenceing a black and white aerial tiff, I can remove the black background, but it makes the aerial slightly transparent. How do I make the aerial less trasnparent again, and maintain the quality of the original tiff?

Qgis


texts: Finding the lat/long bounds of a 1km square (for all of Europe)

texts:
I can code this, given figures/a formula; but I just flat out don't have the G.I.S knowledge for those figures/formula.

I have to come up with some Python which parses the SRTM elevstion data at 30 arc second resolution, and divide it into (approximately) 1km squares.

I then have to examine each elevation in each square and find the highest point in that square.

I am vague aware that I can't just multiply by a constant to convert - but that's my limit.

How can I divide Europe (the political entity) into roughly km squares using Python?


texts: Add a .lyr file to multiple mxd's within a specified folder using python

texts:
I need to create a python script that will look in a certain folder or workspace, and add a certain .lyr file to each mxd in the folder.  I am able create a python script that adds the layer by specifying the path to each .mxd, but ideally, I don't want the name of each .mxd to have to be hard coded in the script.  I've tried using the arcpy.ListFiles("*.mxd") but I can't get my script to read the list when I need to specify the mxd's to point to when identifying the dataframe in the arcpy.mapping.ListDataFrames function.
This is what I have so far...Not sure what I'm missing.

import arcpy
import arcpy.mapping
import os
import sys
from arcpy import env

env.workspace = workspace = "M:\projects\mcag\Project_Site_Packet"
mxdlist = arcpy.ListFiles ("*.mxd")
for mxd in mxdlist:
    mapdoc = arcpy.mapping.MapDocument(mxd)
    df = arcpy.mapping.ListDataFrames(mapdoc, "Layers")
    addlayer = arcpy.mapping.Layer(r"M:\projects\mcag\Project_Site_Packet\Subject Property.lyr")
    arcpy.mapping.AddLayer(df,addlayer,"TOP")
    mapdoc.save()
    del mapdoc, addlayer


texts: climatic data from MODIS or other sources

texts: I'm looking for climatic information in the form of raster like MODIS daily LST (Land Surface Temperature) data. I'm looking for other parameters like precipitation, humidity, air pressure and so on. Does anybody have information where could I find MODIS (or other sources) data for these parameters and what is the best choice for each?


texts: OpenLayers 3: Refreshing OSM layer upon resolution change

texts: Upon changing the resolution of my Map's view with View.setResolution(API doc), the map zooms as expected, but the OSM layer which is included in the Map does not refresh the tiles according to the new resolution level.  Manually (e.g. with the mouse wheel) zoomin in and out once shows the OSM layer as expected.  How can I tell the map (or the OSM layer) to refresh?  I already tried to dispatchChangeEvent() on the OSM layer but it did not change anything.


texts: Convert X, Y into Latitude,Longitude with unknown system

texts:
I have a set of X, Y values that indicate some locations in Hong Kong. But I am not able to convert the X/Y into latitude/longitude since the coordinate system is unknown. Is there any way I can fit the X/Y value into many Geographic information systems to find out which system the data provider is using?

The data I have:

X       Y
840787  816064
839740  819091
839291  819498
840383  816178
835702  818880
833502  816629
837048  815708
838910  823019


texts: I created one polygon shape file to sub-water areas on the river then I created 5 small area layover the first area I created, I need to cut each area to be separate area because when I move the cursor away from them they merge and become only one area that I created first.


texts:
I would like to disaggregate a list of raster files in R. I write the following code but receive an error.

# Import libraries
require(raster)

# Setup filelist and naming convention based off infiles
infiles <- list.files("infilepath here", full.names=TRUE)
names <- extension(basename(infiles), '.tif')

# Set output file naming convention
outpath <- "outfilepath here"
outfiles <- file.path(outpath, names)

# Loop function to disaggregate and write files
for (i in 1:length(infiles)) {
  r <- disaggregate(infiles[i], fact = 2, method='')
  writeRaster(r, paste(outfiles[i]), format ='GTiff', overwrite = T)
}

The error I receive is this:

  Error in (function (classes, fdef, mtable)  :    unable to find an
  inherited method for function ‘disaggregate’ for signature
  ‘"character"’


texts: Zooming Google Maps to fit the extents of a Fusion Table layer

texts:
How can I zoom a Google Maps control to fit the contents of a Fusion Table layer?
The layer only contains pushpin / marker point data.

One method is to perform a REST query on the Fusion Table and to calculate the bounds, and zoom to these (i.e. as illustrated here: http://gmaps-samples.googlecode.com/svn/trunk/fusiontables/fitbounds.html ). However this seems very inefficient - especially for larger datasets, because all of the table data has to be sent to the client and compared.

This question appears to ask the same query, but the answer does not apply. The answer simply queries the data for the current map extents. Nice and efficient, but the opposite of what I want. I want to zoom to display all available data.


texts: How to add coordinate system to shapefile source layer via Web API

texts:
I have a Web API that creates a shapefile on demand for the user to download over the browser. The geometries are all points and in WGS84 projection. I'm running into an issue where the ogr2ogr is throwing this error when I apply a transform to it:

  Can't transform coordinates, source layer has no coordinate system. Use -s_srs to set one.

My question is how do I correctly create the shapefile so ogr2ogr can at correctly guess which coordinate system I'm using. Because I'm not saving the files to disk, I'm actually building the byte array of the .shp file following the shapefile technical description ESRI Shapefile Technical Description

If I use -s_srs or -a_srs to manually assign it with SRID 4326, then ogr2ogr will work.

EDIT: in response to iant's comment, I do have the .prj file in the same directory as the .shp file and they share the same filename.


texts: How to digitize a feature to include elevation data and interpolate slope in QGIS

texts:
I would like to take this archaeological drawing and use QGIS to create a something similar to a DEM or TIN, so that depth/elevation values can be quickly and clearly understood.

I do not know what would be an appropriate workflow to accomplish this. Although lines do represent different depths, they are not regularly spaced like a topographic map. Additionally, there are some areas of the drawing which indicate different depths within the same feature (in the southern wall trench depths of 34, 29, 20, and 23cm are shown together).

At this point, georeferencing is not required, but the ability to do so in the future is one of the reasons I would prefer to use QGIS over other software alternatives. What recommendations do you have for converting this drawing into a DEM or TIN?

Here's an example of the type of drawing I'm working from.


texts: API for Java and JavaScript app integration

texts:
Are there any APIs to make calls or integrate JavaScript functions to a Java written application?  I have come across Rhino's liveconnect.js, but this is uncharted territory for me, and would like a bit more information.  My scenario is using a Java application parallel to a JavaScript application wherein the data that is recorded from the Java app is made available in the JavaScript application as well.  This is all CRM/SOA stuff.

The customer would enter a request in the java application, but the same request is passed to the JavaScript app.

I am using an Esri ArcGIS JS API Web Application with ArcGIS Server and SQL Server 2008.

Node.js could be used, however now, think of it this way;
The ArcGIS JS API web application is the orgin of all spatial components for a call center application which is public facing and internal, there's a Java application that is running parallel to it, we would like to pass any requests originated in the Java application the JavaScript application's database.


texts:
I'm working on hatching lines at 500 foot increments, but I've encountered a confusing error.

I have a start number for the stationing and a number that equals the difference between the start number and an even (effectively 0) point. This ensures the numbers that appear on the hatch labels are at the same 500 foot increments as-builts I am referencing. These numbers in the attribute tables all have a decimal digit value of 2, so they're rounded off as appears on the attribute table.

One of my routes is hatching just fine. The other route hatching seems to have a strange trailing digit error where it isn't going to the even 500 increment numbers. I have checked the attribute properties and decimal digits are 2. Why does it have more significant digits than what I've set it to?


texts: gdal2tiles.py generates wrong tiles / MapTiler the right ones

texts:
I'm trying to generate tiles with gdal2tiles.py and the resulting tiles are stored as 18/132639/164219. If I do the same thing with MapTiler (the GUI variant of gdal2tiles) the tile result is 18/132639/97923. I want the results from MapTiler (97923), but just with gdal2tiles, as I want to do this in a batch process and not with a GUI.

I tried a lot of things already, warping, translating, whatever, and it's getting very annoying, especially as it says that MapTiler is just a graphical interface for gdal2tiles and does basically the same.

gdal2tiles.py -z 17-19 -s EPSG:900913 -p raster overlay-level1.png
gdal2tiles.py -z 17-20 -s EPSG:3785 overlay-level1.png
gdal2tiles.py -p geodetic -s EPSG:4326 -z 18-21 -w none overlay-level1.png

Most of those result in an ERROR 5 Access window out of range in RasterIO() or ERROR 5: Illegal values for buffer size.

This one works, but results in the wrong Y values:

gdal2tiles.py -s EPSG:4326 -z 16-19 -w none overlay-level1.png

Any ideas what's wrong there?


texts: Sequential IDs with field calculator: Pad a prefixed field to specific length

texts:
I have an ID field that is a text datatype. The IDs in the field have a prefix, and are padded with zeros to always be a specific length. Examples:

INSP000098
INSP000099
INSP000100
INSP000101

I would like to use the Field Calculator to generate sequential IDs for new records (records are batch loaded from various sources by the thousands).

I have a python script that almost does this:

Modified from: Create sequential numbers in a field using Python in the Field Calculator

prefix = "INSP0000"
lastnumber=97

rec=0
def autoIncrement():
 global rec
 pStart = 1
 pInterval = 1
 if (rec == 0):
  rec = 1 + lastnumber
 else:
  rec += pInterval
 return prefix+str(rec)

autoIncrement()

Steps:

Manually enter the prefix (with zeros for padding)
Manually enter the last number in the existing IDs
Run the script

Result:

Unfortunately, the tool isn't smart enough to dynamically adjust for the difference in the number of digits between numbers (99 to 100, etc.):

INSP000098
INSP000099
INSP0000100  << Problem: too many zeros/digits
INSP0000101

How can I use the field calculator to generate sequential, prefixed IDs that are padded with zeros to a specific length?

I'm not married to Python; VBScript would be an acceptable alternative.


texts: Cannot set map extent in ArcGIS JavaScript API 3.17

texts:
I have my map set to a custom projection SVY21/3414
This map service only supports Tiled and not Dynamic.
I manually convert the projection coordinates using proj4 and it works fine for my graphics.

I can get an extent object using esri/GraphicsUtils of the graphics on map which seems to be valid. However if I try to set map extent to this object I get error

undefined ""

in dojo.js

If I leave out the projection change and stick to 4236 it works fine.
I'm using version 3.17

Here's some sample code

 require(["esri/graphicsUtils"], function(graphicsUtils) {
    map.setExtent(graphicsUtils.graphicsExtent(arr));
  });

  Aforementioned code only works if I leave the projection change out.
  Following code is how the graphics array is formed

var symbol = new PictureMarkerSymbol();
symbol.setUrl(thisicon);
var latlng = proj4(projection,[longitude,latitude]);
var point = new Point({x:latlng[0],y:latlng[1]});
// "projection" is the wkt string
var graphic = new Graphic(point, symbol, attr, template);
map.graphics.add(graphic);


texts:
I am trying to connect a WMS layer from GeoServer with OpenLayers 3 like below, however, my JavaScript file does not show my WMS layer.

My Chrome says it failed to load response data (WMS). My source is redirected to below URL. This URL works on my server, but not through the JavaScript source.

http://localhost:8080/geoserver/Mars/wms?SERVICE=WMS&VERSION=1.3.0&REQUEST=GetMap&FORMAT=image%2Fpng&TRANSPARENT=true&layers=Mars%3AMercury_basemap&tiled=true&format=application%2Fopenlayers&srs=EPSG%3A4326&transparent=true&WIDTH=512&HEIGHT=512&CRS=EPSG%3A3857&STYLES=&FORMAT_OPTIONS=dpi%3A180&BBOX=-15028131.257091932%2C5009377.085697312%2C-12523442.714243276%2C7514065.628545968

Can someone please help?

It seems like the problem is caused from security settings in GeoServer and I tried to create roles to getread() as well. (http://docs.geoserver.org/stable/en/user/security/service.html)

var layers = [
 new ol.layer.Tile({
 extent: [-13884991, 2870341, -7455066, 6338219],
 source: new ol.source.TileWMS({
  url: "http://localhost:8080/geoserver/Mars/wms",
  params: {'layers': 'Mars:Mercury_basemap', 'tiled': true, 'format': 'application/openlayers', 'srs': 'EPSG:4326', 'transparent': true},
  serverType: 'geoserver'
 }))
 })
];

var map = new ol.Map({
 layers: layers,
 target: 'map',
 view: new ol.View({
   center: [-10997148, 4569099],
   zoom: 4
  })
});


texts: Extract characters from an existing string field to create a new field

texts: I have a string field named CHESS2016_N94S_FL1_C__20160407230029.562GMT_COLOR-8-BIT.JPG and I would like to select characters 16-18 to populate a new field. I've been trying to use the field calculator index function in ArcMap but have only been getting errors.


texts: QGIS 2.18 - How to reference a svg symbol which is in the qml file's directory?

texts:
Relative path for svg symbols in .qml files works as expected if the svg file is in another directory than the .qml file's.
Example :

<prop k="name" v="../icons/symbol.svg"/>

However, if the svg file is in the same directory as the .qml file, it doesn't.
I tried with this syntax :

<prop k="name" v="symbol.svg"/>

And this one :

<prop k="name" v="./symbol.svg"/>

With no success...
Is it a bug or do I miss something ?


texts: Is it possible to redefine a projection in QGIS?

texts:
Is it possible to redefine a projection in QGis?

It seems like the shape files offered by Agiv (a Flemish agency) differ slightly from the definition that QGis has of EPSG:31370.

I'm NOT looking to create a custom projection, only to replace the existing one.

Update: the definition Agiv (official Flemish agency) uses, becomes this in QGis if you add a shapefile. It is considered a user defined coordinate system, generated CRS, USER:100000.

Generated CRS (+proj=lcc +lat_1=49.8333339 +lat_2=51.16666723333333 +lat_0=90 +lon_0=4.367486666666666 +x_0=150000.01256 +y_0=5400088.4378 +ellps=intl +towgs84=-106.869,52.2978,-103.724,0.3366,-0.457,1.8422,-1.2747 +units=m +no_defs)

The definition of the EPSG:31370 in QGis for EPSG:31370 is:
+proj=lcc +lat_1=51.16666723333333 +lat_2=49.8333339 +lat_0=90 +lon_0=4.367486666666666 +x_0=150000.013 +y_0=5400088.438 +ellps=intl +towgs84=-106.869,52.2978,-103.724,0.3366,-0.457,1.8422,-1.2747 +units=m +no_defs

So there's only a very very minor difference, but it seems QGis won't recognize it by default. That's what I'd like to change. (I know I could define it manually). That's how the files are provided, and how they are recognized in ArcGIS, which shows the same False_Easting and False_Northing.


texts:
I am new to remote sensing. I downloaded a landsat image from earth explorer. It was in tar.gz format. I used 7-zip to extract the files. I got 7 tiff images and 1 jpg image.

I added those 7 images in arc map 10.4.1 and made a composite layer. I made it permanent as the one generated first is always a temp file. But i found that （both temp as well as ）the permanent file has only band 1 mentioned in front of red green and blue.

When I try to change the band numbers I get 7 bands all of which are numbered as 1. Although I can still make it an fcc image by clicking on the fourth, third and the second band in the layers window corresponding to red green and blue bands respectively.

But why the bands are not numbered from 1 to 7?


texts: D3 - Loading data from JSON doesn't work

texts:
I'm starting to learn D3 and it seems like a powerful framework. I'm trying load the data from json. I can see data but nothing happen after loading data into D3 library. There are no error messages in the JavaScript console.

Here's my code:

<html>
  <head>
        <meta charset="utf-8">
        <title>Map</title>
  </head>
<style>

</style>
<body>
<script src="http://d3js.org/d3.v3.min.js"></script>
<script>
    var width  = 1000;
    var height = 1000;

    var svg = d3.select("body").append("svg")
        .attr("width", width)
        .attr("height", height)
        .append("g")
        .attr("transform", "translate(0,0)");

    var color = d3.scale.category20();

    d3.json("examples/us.json", function(error, root) {
        console.log('read');
        if (error)
            return console.error(error);
        console.log(root);
         var group = svg.selectAll("g")
            .data(root)
            .enter()
            .append("g")

          var projection = d3.geo.mercator(); //.scale(7300).translate([0,1980]);
          var path = d3.geo.path().projection(projection);

          var areas = group.append("path")
            .attr("d", path)
            .attr("class", "area")
            .attr("fill", "steelblue")

    });

</script>

</body>
</html>

Plunker : http://plnkr.co/edit/wHIsaB?p=preview


texts: How to export vector layer containing polylines into polygons in Google Earth

texts: I'm using QGIS 3.16.5 Hannover and I already knew how to convert vector layer (from AutoCAD containing polylines) into KML format. When I export it as KML and open it up to Google Earth Pro, the resulting KML object is a path (not a polygon as I would want it to be). There is an option to set geometry type as polygon (in the export vector window) but it doesn't make any difference. How do I set it as polygon and not a path?


texts: Recalculate a Field with data from another field, ignoring NULL values

texts:
I have a Feature Class in a Geodatabase with an field for pipe diameter. I also have a shapefile with a field for diameter. The shapefile is from field verification work that verified a variety of pipe characteristics. I have joined the shapefile to the Feature Class in ArcMap, see the image below. I would like to recalculate using Field Calculator and the Python Parser the first diameter field with that information in the second diameter field, any null values would keep the original data in the first diameter field. Through some research I have hodgepodged the following script that fails upon execution when using the Filed Calculator. ArcMap spits out an Error:

  Python syntax error: Parsing error SyntaxError: invalid syntax (line 1)

def IgnoreNull (GravitySewer.Diameter, assets.Diameter)
 if assets.Diameter is None:
      return GravitySewer.Diameter
 else:
      return assets.Diameter


texts:
My request in postgis :

create table "schema1"."table1" as select p.field1 ,p.field2, l.field3, l.id from
"schema1"."table2" p inner join "schema1"."table3" l on (p."GEOM" && l."GEOM" and st_intersects(st_buffer(p."GEOM",0),st_buffer(l."GEOM",0)));

I want to execute this request on Oracle so st_intersects is SDO_ANYINTERACT but I don't find for && operator.

What is the equivalence of && operator in oracle sdo ?


texts: Recomputing footprints on old satellite imagery

texts: I am trying to create an image service out of some 1976 satellite imagery. But first I need to get all of the no data cells out of my mosaic dataset. I have tried to recompute the footprints on my mosaic dataset in ArcMap and only some of the footprints are changing on the tiles. I don't understand why this is happening, and I am not sure how to fix it. Any suggestions would be helpful.


texts: Is there a way to embed a Leaflet map into a Wordpress theme website?

texts:
I recently decided to make my own website, but chose to use a Wordpress theme to save me time in designing. Unfortunately it primary uses php which I am not familiar with and end up breaking my page after attempts to edit. Is there a relatively simple way to place my webmap from QGIS2Web in wordpress? Potentially in my portfolio or another page?

This question here is fairly similar, but I'm not working in html.


texts: Estimate area from raster different resolution

texts:
I am trying to calculate the area that each pixel in the raster has, in order to know the total area per country. This area refers to those pixels equal to 1 which refers to where protected areas in the world should be. According to the authors, they constructed this raster with a resolution of 20km at the equator decreasing it towards the poles. The raster is in GCS_WGS_1984.

I have tried raster statistics with a shapefile of all the countries in the world. I have now the count of pixels equal to 1 and 0 but I do not have a resolution to multiply with.

It can be found here: https://drive.google.com/drive/folders/1jmqbUj95iYCH_jTfoKfi17-J-rJAgq1f?usp=sharing

Protected areas = 1 (blue), not protected = 0 (black)

Protected areas = 1 (blue), not protected = 0 (black)


texts:
I am writing a python plugin for QGIS 1.7 (started it on QGIS 1.6 but switched to 1.7 recently).
Since I switched to 1.7 my unload function doesn't seem to be working anymore. Whenever I implemented changes to the code I use the "Plugin Reloader" plugin to load the newest version of my code.

In the MenuBar this causes an additional, new instanec of the plugin to load while it doesn't unload the old ones. After doing this a couple of time I end up having a list of endless outdated outdated references to my plugin. For the menu bar icon the same does NOT happen...

The code worked fine in 1.6 and I seem to be correctly following all the online examples.
It shouldn't be a problem, that I'm creating two actions, at least the PluginReloader plugin seems to do the same.

def initGui(self):
    [...]
    # GemEB Tool
    self.action_tool = QAction(QIcon(iconpath), u"GemEBTool", self.iface.mainWindow())
    self.action_tool.setWhatsThis(u"GemEBTool")
    self.iface.addPluginToMenu("&GemeindeEnergieBeratung 1.5", self.action_tool)
    self.iface.addToolBarIcon(self.action_tool)
    QObject.connect(self.action_tool, SIGNAL("triggered()"), self.run_tool)

    # Solarrechner
    self.action_solarrechner = QAction(QIcon(iconpath), u"GemEB Solarrechner", self.iface.mainWindow())
    self.action_solarrechner.setWhatsThis(u"GemEB Solarrechner")
    self.iface.addToolBarIcon(self.action_solarrechner)
    self.iface.addPluginToMenu("&GemeindeEnergieBeratung 1.5", self.action_solarrechner)
    QObject.connect(self.action_solarrechner, SIGNAL("triggered()"), self.run_solarrechner)

def unload(self):
    # GemEB Tool
    self.iface.removePluginMenu("&GemEBTool",self.action_tool)
    self.iface.removeToolBarIcon(self.action_tool)

    # Solarrechner
    self.iface.removePluginMenu("&GemEB Solarrechner",self.action_solarrechner)
    self.iface.removeToolBarIcon(self.action_solarrechner)


texts: I am new in GIS, but I need extract addresses (city, street, housenumber, lat, lon) from my osm file. Can you give me advice what is the best way to do that? I found project openaddresses, but a lot of addresses of my area is missing.


texts:
After running a lot of scripts in QGIS' python console, I have more than 10 dictionaries and lists, and I need these for further analysis. Is there a way to save python session, with all these variables?

In general, my question is how can I export my variables from python console to something, so that they can be loaded in again at a future time?

Thanks


texts: Large scale geoprocessing through ESRI REST interface?

texts:
A typical limitation of web based processing is that urls are limited to 2 thousand characters for a GET request.  From my understanding of how the ESRI Javascript API works, JSONP is utilized to interact with ArcGIS Server Geoprocessing services via the REST API, requiring the use of a GET request?

This seems to pose a problem if I'm attempting to send large quantities of user defined runtime data (> 100 coordinate pairs and several other user parameter strings) to be worked upon in a geoprocessing service via the standard ESRI Javascript API.

What's the best way to tackle this problem?

I'm thinking about writing a proxy service that will receive GET requests, but will issue POST requests to the REST API.


texts: Slow performance when selecting features

texts:
I am trying to use the following code to select the features that have the attribute values that I pass in. The performance is really slow when it comes to a lot of values e.g >16k values to match features.

Is there something I'm missing or a better way of doing this?

// Assuming values.Any() is true:

var featureSelection = (layer as IFeatureSelection);

var whereClause = new StringBuilder();
var count = 0;
while (true)
{
    var values = valuesToSelect.Skip(count*1000).Take(1000).ToList();
    if (!values.Any()) break;
    var thousandsClause = new StringBuilder();
    foreach (var value in values)
    {
        thousandsClause.Append(",");
        thousandsClause.Append(value);
    }
    thousandsClause.Remove(0, 1);
    thousandsClause.Insert(0, fieldName + " IN (");
    thousandsClause.Append(")");
    whereClause.Append(" OR " + thousandsClause);

    count++;
}
whereClause.Remove(0, 4);
IQueryFilter queryFilter = new QueryFilter();
queryFilter.WhereClause = whereClause.ToString();

featureSelection.SelectFeatures(queryFilter, esriSelectionResultEnum.esriSelectionResultNew, false);
featureSelection.SelectionChanged();

Sample WHERE IN clauses:

761 values to select

122,072 values to select

Edit--

So far with recommendations from some of the comments, I've improved the code and the results are faster but still not desirable. Here is the code I've begun trying using the same WHERE IN clause:

IQueryFilter queryFilter = new QueryFilter();
queryFilter.WhereClause = whereClause.ToString();
queryFilter.SubFields = "OBJECTID";

var featureLayer = layer as IFeatureLayer2;
if (featureLayer == null) continue;

IFeatureClass featureClass = featureLayer.FeatureClass;
var dataset = (IDataset) featureClass;
IWorkspace workSpace = dataset.Workspace;
ISelectionSet selectionSet = featureClass
.Select(queryFilter, esriSelectionType.esriSelectionTypeHybrid,esriSelectionOption.esriSelectionOptionNormal,workSpace);

var featureSelection = layer as IFeatureSelection;
if (featureSelection != null)
{
    featureSelection.SelectionSet = selectionSet;
    featureSelection.SelectionChanged();
}


texts: QGIS attribute table displaying wrong features?

texts:
In QGIS 2.18.1 I set up Attribute table to 'Show features Visible On Map' but in the form view in the list are definitely more features than visible in the map.

In the screenshot below two features are visible in the map (Layer Road_Erlangen_2016_10_10):

The orange one on the top right corner and
The (selected) yellow one in the middle

... but the screenshot reveals 5 (when I scroll down ~15) feature in the left hand list:

In another case I have too few features in the list (4 in the map, one in the list):

My first thought was that the slow reacting WMS basemap aerial photo layer makes such problems, but turning it off changes nothing.

Data source is a shapefile, I created a spatial index under Layer properties and I'm in editing mode (where leaving editing mode changes nothing). Even 'Reload the table' button has no effect.

Am I missing some settings or what may be my problem?


texts: Geometry errors when importing shapefile with shp2pgsql PostGIS

texts:
I'm trying to import the Census State Boundaries from the census data site although I get errors in the geom data every time I import in the same rows.
I've tried shp2pgsql and ogr2ogr
shp2pgsql -I -d -e /Users/cj/Projects/geospacial_data/cb_2019_us_state_500k/cb_2019_us_state_500k.shp public.state | psql -h localhost -d my_db -U cj

and this
ogr2ogr -f "ESRI Shapefile" /Users/cj/Projects/geospacial_data/cb_2019_us_state_500k/cb_2019_us_state_500k.shp PG:"host=localhost user=cj dbname=mydb" "state"

The errors show up here in the geom column (in the screenshot) showing the loaded size. Are these errors just part of the file? or am I doing something wrong? This seems to crop up when I do Counties and zip codes as well.

I'm running

PostGIS 3.0.1 (homebrew)
PostgreSQL 12.3 (homebrew)
GDAL 2.4.4 (homebrew)
MacOS 10.14.6


texts: Missing coordinates in the attributes list after importing osm-Files

texts:
I imported an osm-File in QGIS 3.12.1 (osm File was created with JOSM).
Then I wanted to export the XY coordinates (Long.-Lat. EPSG4326 => XY UTM). But it seems that the coordinate information of the points is missing.... I exported the layer to csv & dxf but there are no points.
Currently I used the plugin MMQGIS to export the coordinates to csv, import the csv to a new layer and then exported the layer with UTM and Field "as_XY".
I guess there is a way to get those informations directly out of the osm-Data?!
Thanks.


texts: Error when creating data with special characters and numbers

texts:
My code will not project and store the feature layer in the scratch.gdb.
"ERROR 999999 ... The table name is invalid. The table was not found. [17_5330_51961Lake]"

The interpreter doesn't like "-" characters or beginning or ending in a number. Does anyone know of a work around?

import arcpy
import os

mxd = arcpy.mapping.MapDocument("Current")
mxdpath = mxd.filePath
(filepath, filename) = os.path.split(mxdpath)
MCP_Path = str(filepath)

srl = os.path.join(MCP_Path, r"CoordinateSystems\NAD 1983 UTM Zone 17N.prj")

WBYLID = "17-5330-51961"

WBYLIDlabel = WBYLID.replace("-","_")

LakeOut = os.path.join(MCP_Path + "\OutputData\Scratch.gdb", WBYLIDlabel + "Lake")
arcpy.Project_management(fl, LakeOut, srl)


texts: Converting 3dr file into EsriASCII Raster for ArcMap

texts:
I am taking a GIS course and have to convert a 3dr file that is a hardened image of combined fuzzy soil data. I then have to convert it to asc and open it in ArcMap.
My created 3dr image looks like the example, but when I convert it, I lose pixels containing "Basco" soil.
I have done all of this through the class-provided SoLIM program.
This is the before and after of the conversion:

Is there another method of conversion? Or is there something I am missing? I have done it twice now from scratch and got the same results.


texts: Is it true that WAAS provides vertical and horizontal error correction but RAIM only provides horizontal correction?

texts:
Is it true that WAAS (or other SBAS) provides both vertical and horizontal error correction but RAIM only provides horizontal error correction but absolutely no vertical error correction beyond the acquisition of 4 GNSS satellites?

If this is true, then is the grade of horizontal correction from WAAS comparable to RAIM or is WAAS also much more reliable and/or accurate horizontally than RAIM?


texts: Set default projection of geometry objects in arcpy

texts:
I have a script that creates a polyline geometry from scratch using some math to generate coordinate pairs. At the same time, there is a user-defined polygon shapefile that is used to "erase" the polylines. Doing this using the low-level geometry methods (e.g. http://bit.ly/1vBMSjP under the "Methods" section) is much, much faster than using feature classes (e.g. arcpy.Erase_analysis). However I am trying to figure out the fastest way to assign a projection to the un-projected polyline geometry coordinates. Is there a way to set a default projection for all processing in arcpy?

I want to use this geometry operator to erase the polyline geometry:

lineGeometry.difference(polygonGeometry)

However this fails if both geometry objects do not have the same projection. So my question is, is it possible to set a default projection for all geometry and arcpy operations, so that a new geometry created from scratch will be assumed to have that projection?

For reference, below are the three strategies that I have used so far to assign a projection to the geometry, since I have not been able to figure out how to assign a default projection.

use the "Define Projection" function

arcpy.DefineProjection_management(lineGeometry, SpatialReference)

where lineGeometry is a Polyline geometry object created from calculated coordinates, and SpatialReference is read from the polygon shapefile. This is very slow (possibly because the geometry is converted to a feature class internally, then back to a geometry).

copy the line geometry into a pre-existing polyline feature class with a defined projection, then copy it back to a geometry list object. Oddly enough this is much faster than #1 above.

g = arcpy.Geometry()
arcpy.CopyFeatures_management(lineGeometry,linesTemp)
lineProjected = arcpy.CopyFeatures_management(linesTemp, g)
currentLine = linesTemp[0]

Use a polyline feature class with a defined projection, with only one polyline feature in it. Then, copy the line geometry directly into the SHAPE field of the existing polyline feature, overwriting the existing geometry. Then copy it back to a geometry list object. This is the fastest method I have found so far.

g = arcpy.Geometry()
with arcpy.da.UpdateCursor(linesTemp,["SHAPE@"]) as cursor:
    for row in cursor:
        row[0] = lineGeometry
        #Update the current cursor with the changes
        cursor.updateRow(row)

#Copy the projected line feature back out into a geometry list
lineProjected = arcpy.CopyFeatures_management(linesTemp, g)
#Get the actual line geometry from the geometry list
lineGeometry= lineProjected[0]


texts: Missing Conversion | Rasterise option in raster menu

texts:
I'm using QGIS.
Would like to use Raster > Conversion > Rasterise from drop down menu.
Raster drop down menu is smaller today than it was yesterday though... why?

How can I access the missing options?

Note, I'm aware rasterisation can be done through processing toolbox, but I don't want to use that for various reasons.


texts: I need a attribute field to be auto generated incrementally. What to do

texts: I have created a FGDB with many layers. Each layer I need a attribute FeatureID. It should be auto generated with increment number. How can I do this?


texts: Add attributes to polygons that are adjacent to a polygon from a different layer?

texts:
I have a shapefile with polygons which represent a raster which represents all countries on Earth (~70,000 cells) and one shapefile layer from NaturalEarthData which represents the worlds oceans.

How can I add an attribute to each raster cell that is adjacent to the outer border of my open shapefile or, in other words, how can I add a specific attribute to those polygons that entails the coastline?


texts: Find intercept WGS84 coordinates between two segments defined by (Lat1-Long1 and Lat2-Long2) and (Lat3-Long3 and Lat4-Long4)

texts:
I need to determine the WGS84 coordinates for a point that is defined by the interception of two segments represented by two sets of WGS84 coordinates each (lat1-long1, lat2-long2) and (lat3-long3, lat4-long4).
Is there a free application that I might be able to use or even an excel formula?


texts: Is it possible to create standalone installers for QGIS plugins?

texts:
Plugin in QGIS are managed by the internal plugin manager, and this requires active (and fast) internet connection. Another way to do it offline is answered here How to install a QGIS plugin when offline? which requires technical knowledge.

Is there not any way we can get standalone installers for individual plugins just like the qgis package?


texts: Setting Embed Interop Types in ArcObjects to True or False?

texts:
I took over a new project with some legacy ArcObjects code and the original developer set the ArcObjects references "Embed Interop Types" property in Visual Studio to a mixture of "True" and "False".

From my experience, the differences mostly seem to be some syntax changes with things such as creating new objects.

Is there a reason to choose  "True" rather than "False"?


texts: How to get back discrete raster legend symbols in QGIS 3.16

texts:
In the recent QGIS version, layer tree symbols for pseudocolor continuous raster palette have changed, so that it shows as bar with continuous color range with only minimum and maximum value of the raster for the defined color range displayed.
But for raster showing for example concentrations of chlorophyll-a in inland water as in image below, it is important to be able to distinguish if the concentration value is 0 or 15 or 30 or 60 mg/m3 or much higher. Therefore, I created a color legend with continuous color change in the map, which shows well how the concentrations change, but showing the discrete values of 0, 15, 30, 60, and 600 in the legend, each with different distinct color and linear transition between them.
Is there a way to switch the layer tree style back to the old rendering (i.e. similar to what is visible in the Layer properties/Style dialog)? Or at least show more distinct values on the scale in the bar than just minimum and maximum.
I see many users wanted the change the way it is now, since I am basically wanting the opposite of what was requested for older versions for example here.


texts: How to restrict only a Layer Group instead of giving the whole service when loading a WMS on qgis?

texts:
I created a workspace, a store and some layers on GeoServer. Then I created a layer group and added a few layers.
When I try to load a specific layer on QGIS I use the following URL:

http://localhost:8080/geoserver/workspacename/layername/wms

and everything works fine, I do the connection and is listed only the layer I want and not the whole service layers.
Now I want to do the same for a layer group, I want that only a specific layer group be listed when I click on connect button.
However, if I do the same as above:

http://localhost:8080/geoserver/workspacename/layergroupname/wms

it gives an error.

In OpenLayers, the layer group is correctly loaded when I do the layer preview in GeoServer, and the resulting URL is the following:

http://localhost:8080/geoserver/workspace/wms?service=WMS&version=1.1.0&request=GetMap&layers=workspace:layergroupname&styles=&bbox=-25.1745333333,30.02271667,-15.85316667,37.2715666667&width=768&height=597&srs=EPSG:4326&format=application/openlayers

but using this URL on QGIS when loading the service, the whole service layers are listed and not just the group layer I want.

How can I restrict to a layer group when loading a WMS service into QGIS?


texts: ogr2ogr gml file loaded but not in the database

texts:
I am trying to import gml files to postgres by following the steps on AstunTechnology/Loader

loader.config file is updated as;

ogr_cmd=ogr2ogr --config GML_EXPOSE_FID NO -append -skipfailures -f PostgreSQL PG:'dbname=postgres active_schema=public host=data01 user=postgres password=123456789' $file_path

The python loader.py loader.config file execute successfully as,

Config: {'tmp_dir': 'U:/Data/AB/data/Temp1', 'out_dir': 'U:/Data/AB/data/Temp', 'gfs_file': '../gfs/osmm_topo_postgres.gfs', 'ogr_cmd': "ogr2ogr --config GML_EXPOSE_FID NO -append -skipfailures -f PostgreSQL PG:'dbname=postgres active_schema=public host=data01 user=postgres password=123456789' $file_path", 'post_cmd': '', 'debug': 'True', 'src_dir': 'U:/Data/AB/data', 'prep_cmd': 'python prepgml4ogr.py $file_path prep_osgml.prep_osmm_topo'}
Processing: U:/Data/AB/data\AddressBasePremium_FULL_2019-09-12_001.gml
Prepared file: U:/Data/AB/data/Temp1\tmp8aypa_\AddressBasePremium_FULL_2019-09-12_001
Prep command: python prepgml4ogr.py U:/Data/AB/data\AddressBasePremium_FULL_2019-09-12_001.gml prep_osgml.prep_osmm_topo
Loading: U:/Data/AB/data\AddressBasePremium_FULL_2019-09-12_001.gml
OGR command: ogr2ogr --config GML_EXPOSE_FID NO -append -skipfailures -f PostgreSQL PG:dbname=postgres active_schema=public host=data01 user=postgres password=123456789 U:/Data/AB/data/Temp1\tmp8aypa_\AddressBasePremium_FULL_2019-09-12_001
Loaded 1 file

However I cannot see a table created in the public schema. Am I missing something here? I believe the table name would be similar to the file name.


texts: Splitting vector layer into shapefiles (not geopackage) using QGIS?

texts:
Using QGIS 3.6.3, I have a layer with ten polygons. I want to use split vector layer to split the file into ten shapefiles, But the result output file is a Geopackage.

How can I get ten shapefiles by using split vector layer?


texts: Calculating sequential numbers into sorted table using ArcGIS Desktop?

texts:
Is there a way to calculate a sorted field with sequential numbers?  I have seen Sorting feature class to calculate sequential ID field using ArcGIS Field Calculator? that outlines how to calculate sequential numbers, but this is always calculated on FID order, not on sorted order.

#Pre-logic Script Code:
rec=0
def autoIncrement():
    global rec
    pStart = 1
    pInterval = 1
    if (rec == 0):
        rec = pStart
    else:
        rec += pInterval
    return rec

#Expression:
autoIncrement()

An example of what I'm trying to do.  I've used an advanced sort to sort by year, month, day, and now want to have sequential numbers in the Seq field.  You'll see that my OBJECTID field is not in order, so the above code won't work.

Can this be done either in the Field Calculator or using an Update Cursor in arcpy?


texts: (R)SAGA: saga_cmd 2.1.0 error? inconsistency?

texts:
frequently I am using (R)SAGA, but after SAGA became >2.1.0 the R package is no longer working. The last update for RSAGA was 07/2013.
Thus I downgraded SAGA to 2.1.0 (first try was skipping the ppa of johanvdw, second try was from source) in order to be able to use RSAGA again.
Actually this was successful I thought because saga_gui ran well.
Only the more important saga_cmd brings up this:

In bash:

[jejune:~]$ saga_cmd
Error: no arguments for saga call
Error: library

62 loaded module libraries (567 modules):
- garden_webservices
- ihacres

...

- shapes_grid
- table_tools
- ta_compound

    type -h or --help for further information

    Inconsistency detected by ld.so: dl-close.c: 762: _dl_close: Assertion `map->l_init_called' failed!

In R I can load the package but all functions come up with this error message:

    rsaga.get.modules('io_gdal')
    Error: module
    Inconsistency detected by ld.so: dl-close.c: 762: _dl_close: Assertion  `map->l_init_called' failed!
    Error in system(command, intern = intern, ...) : error in running command
    >

So either there is a workaround to get RSAGA runnning with SAGA >2.1.0 or maybe someone is able to tell me what the error messages mean (and hopefully how to get rid of them)?

system in use:

Ubuntu 14.04
R version 3.0.2
SAGA 2.1.0
RSAGA 0.93-6


texts: Batch processing not available in Graphical Modeller QGIS

texts:
I'm currently running QGIS 3.6. I'm trying to set up a model that uses the GDAL 'Clip Raster by Mask' tool to extract individual rasters from a DEM. I'm using geopackage polygons as the mask layer. I have no issues carrying out a batch processing when outside of the graphical modeller.

I can't see any option to select 'batch processing' in the graphical modeller tool interface.

Is there a way to carry out batch processing within the graphical modeller?


texts: Querying PostGIS and OSM to get only major stations and subway of UK?

texts:
I am working with Great Britain OSM data. I have read the classification of UK railway stations. There are a lot of classes in there such as subway entrance, halt, station and many others. I want to get only the major stations and subway of UK. I am using the following Query but it is not returning all station points. I have compared the result with OSM maps service.

SELECT way, name, railway from planet_osm_point WHERE railway IN ('station', 'halt', 'entrance') AND name is not null

What more types should I include?


texts: Increasing/reducing features proportionally using QGIS?

texts:
Is it possible with QGIS to increase or reduce features (polygons or lines already drawn) proportionally? Preventing me to redraw several times new polygons :)

See picture above: Making it possible to resize for example the green polygon to the size of the red polygon (the polygons are build by the plugin "Rectangles ovals digitizing").

In ArcGIS, this can be accomplished by using the Scale tool to resize an individual feature.  I haven't found a solution in the QGIS-Guide.


texts:
A have a number of shapefiles and I want to export their attribute table to .dbf and .mdb (Microsoft Access file). How can I do this with ArcGis or with any other tools?

Notice that I want to keep their actual field names.

E.x. When I use "table to table" tool the "FID" field change to "OID". I do not want this to happen.


texts: GDAL/BUILD can't get result files Geojson to mbtiles

texts:
i try to 2.55gb of GeoJson (maybe 200 over layers in one file) to vector Mbtiles, but fail to create .mbtiles.

so i try convert each layers(.shp) to geojson(over 200 files),
then create each .mbtiles. but works it.

so i try to edit gdal source for console debug message,
end the compile, try to geosjon to mbtiles(Previously successful file),
but now... cannot error, but cannot create mbtiles file!

summary.

2.55gb of Geojson to mbtiles fail(using gdal 2.3.0ver Binaries)

what is problem? many layer? or json's size?

custom gdal(just add console print) can't create mbtiles file.

What should I identify?


texts: ArcGIS - Spatial ETL Tool - Checking fields in input shapefile

texts:
How to check if the input file fields have the expected name and type?

Which transformer will be useful to validate names and types of attributes of reader.

edit:
I have a shapefile that transforms using ETL to gdb.
But at first I want to check if the attributes of this shapefile is required by me names and types.
For example: I want to check if the file road.shp, which is my "reader" has attributes:
'NAME' [type: char, width: 265],
'NUMBER' [type: integer].
How to control input shapefile when I will have a next file to be transformed.

I'm using ArcGIS 10.2.2 with FME 2014