texts: Converting tuple into records in blank feature class using ArcPy?

texts:
Here is the code block that I am using with ArcGIS 10.2 for Desktop, it all works great except for the bottom part. I am not sure how to turn the tuple into records on a blank feature class. The tuple looks like this: ('BRONTE DR', 'BRONTE DR', 'RICHARDSON DR', 'RICHARDSON DR'). I want each part of the tuple to be a new record in the blank feature class.

import re
import arcpy
from arcpy import env
arcpy.env.overwriteOutput = True
if arcpy.Exists(r"J:\STW\Divisions_&_Branches\MSMD\Branches_&_Sections\GIS\users\Anault\intersections_software_development\OUTPUT.gdb\street_assets2"):
    arcpy.Delete_management(r"J:\STW\Divisions_&_Branches\MSMD\Branches_&_Sections\GIS\users\Anault\intersections_software_development\OUTPUT.gdb\street_assets2")
arcpy.env.workspace = r'J:\STW\Divisions_&_Branches\MSMD\Branches_&_Sections\GIS\users\Anault\intersections_software_development\OUTPUT.gdb'
table1 = r"street_assets"
field1 = r"STREET_OR_INTERSECTION"
def unique_values(table, field):
    with arcpy.da.SearchCursor(table, [field]) as cursor:
        return sorted({row[0] for row in cursor})
newtuple = unique_values(table1,field1)
newtuple = str(newtuple)
newtuple = newtuple.split("}{")
newtuple = ([s.replace('\'', '') for s in newtuple])
newtuple = ([s.replace('"', '') for s in newtuple])
newtuple = ([s.replace('u', '') for s in newtuple])
newtuple = ([s.replace('[', '') for s in newtuple])
newtuple = ([s.replace(']', '') for s in newtuple])
newtuple = ([s.replace('{', '') for s in newtuple])
newtuple = ([s.replace('}', '') for s in newtuple])
out_path = r"J:\STW\Divisions_&_Branches\MSMD\Branches_&_Sections\GIS\users\Anault\intersections_software_development\OUTPUT.gdb"
out_name = "street_assets2"
geometry_type = "POINT"
template = r"J:\STW\Divisions_&_Branches\MSMD\Branches_&_Sections\GIS\users\Anault\intersections_software_development\OUTPUT.gdb\street_assets"
has_m = "DISABLED"
has_z = "DISABLED"
spatial_reference = arcpy.Describe(r"J:\STW\Divisions_&_Branches\MSMD\Branches_&_Sections\GIS\users\Anault\intersections_software_development\OUTPUT.gdb\street_assets").spatialReference
arcpy.CreateFeatureclass_management(out_path, out_name, geometry_type, template, has_m, has_z, spatial_reference)
newtuple = tuple(newtuple)
print newtuple

## Problem Here

cursor = arcpy.da.InsertCursor(r"J:\STW\Divisions_&_Branches\MSMD\Branches_&_Sections\GIS\users\Anault\intersections_software_development\OUTPUT.gdb\street_assets2", ['STREET_OR_INTERSECTION'])
i = 0
for i in newtuple:
    cursor.insertRow([i])
    i = i + 1


texts:
I'd like to symbolise a point dataset with a square (or circle) using two colours. This is useful for displaying values collected on both sides of a road (for example).

I created two symbols, a left side and right side of the square as shown below:

SYMBOL
    NAME 'SQUARE_L'
    TYPE vector
    FILLED true
    POINTS
        0 0
        0 1
        0.5 1
        0.5 0
        0 0
    END
END
SYMBOL
    NAME 'SQUARE_R'
    TYPE vector
    FILLED true
    POINTS
        0.5 0
        1 0
        1 1
        0.5 1
        0.5 0
    END
END

Then in my LAYER I can create a CLASS using the two different STYLES:

CLASS
    STYLE
        SYMBOL "SQUARE_L"
        COLOR [LEFT_CLR]
        SIZE 10
    END
    STYLE
        SYMBOL "SQUARE_R"
        COLOR [RIGHT_CLR]
        SIZE 10
    END
 END

However the features do not display correctly - they are not aligned correctly and are not of equal size (due to an overlap?). I thought I had resolved it by setting a "bounding box" in the symbol to make sure they were equal size:

SYMBOL
    NAME 'SQUARE_L'
    TYPE vector
    FILLED true
    POINTS
        0 0
        0 1
        0.5 1
        0.5 0
        0 0
        -99 -99 # increase bounding box so shape is correct size
        0 0
        1 0
        1 1
        0 1
        0 0
    END
END

However I then noticed that the sides of the square were displayed with the wrong values. The MapServer symbol collection at http://github.com/mapserver/mapserver/wiki/SymbologyExchangeVector has some nice examples but none of which describe how to do this. Any ideas on how to achieve this in MapServer?


texts: Is there a OSM tile server that allows you to style like google api?

texts:
The google maps API allows you to style the map you get back. Is there an OSM equivalent?  This question asked how google did it How are Google serving up their Styled Maps?
has there been any progress since then in the OSM community?


texts: Calculating overlapping of lines using QGIS?

texts:
I have a GeoPackage that includes series of lines, each line corresponds to the itinerary of a person from point A to point B. I would like to know how to calculate the overlapping of lines (the superimposition has been made by the sna) and how to get to know which one are the most used itineraries.

I've already tried to use QChainage to get a heat map but I wondered if there was another way.


texts: Invalid SQL-Statement error when querying file-geodatabase table date field with arcobjects

texts:
I want to query a DateTime Field in a table in a file-geodatabase and I keep on getting the error "invalid sql statement" (Error Code: -2147220985) while that sql statement does not give me an error on an empty table with the same configuration on an sde (even though there is no data in that table I think if the sql statement was wrong it should produce the same error there).

I use ArcObjects .Net SDK 10. I construct my whereclause like this, passing in two DateTime Objects as arguments. I've tried various ways of formating them.

string whereClause = string.Format("dat >= '{0}' AND dat <= '{1}'", fromDate.ToShortDateString(), toDate.ToShortDateString());

I then create an IQueryDef on a file-geodatabase workspace and set the previous whereclause. I then call the Evaluate() method of the IQueryDef and that is where it crashes. I was able to successfully launch other where-clauses on that same table in that same workspace but without taking the Date-Field into account. At the moment I don't see where the problem is. Is there a difference between accessing Date columns on a file-geodatabase vs a SDE with SQL out of ArcObjects? Or is something wrong with the way I use SQL here?


texts:
I would like to vectorize a PNG containing lines (DN=1 for line, DN=0 for transparent background). See below. As the data is georeferenced I'm asking it here. I would like to automate the vectorization and I'm currently looking into Python/GDAL/QGIS solutions with a vector result in the form of a shapefile. I found gdal_polygonize but this tools creates a polygon shapefile while I would like to create polylines. Does anyone have an idea how to make this work?

UPDATE: Searching for an open source solution, so without using ArcGIS like this solution.


texts:
I am getting the error-000539: Error running expression: rcexec(). I've seen that other people have this same issue but I am not sure how to fix mine. I just want to average the values of 10 rasters, which are .tiff files. This code worked in my older version of ArcGIS but I had to update to the recent one and now it fails.

Here is the full error message:

  ERROR 000539: Error running expression: rcexec()

  Traceback(most recent call last)P

  File "", line1, in

  File "", line 8, in rcexec

Here is a picture of my raster calculator code.

This calculation is one step in a model since I need to do this many times. Here is a picture of that too if it helps.


texts: Labeling runs out of map extent and onto borders when exporting to PDF from ArcGIS for Desktop?

texts:
When I have tried to export several of the maps I have made but I have a problem with labeling that runs out of the map extent and onto the borders.

I am exporting to a PDF and opening the file in Adobe.

Is there a way to stop this from happening?

I am using ArcMap 10.1.


texts: Save vector layer and update .qgs project file (PyQGIS)

texts:
Need help to solve the following use case:

QGIS started with new project (no layers)
New layer is created programmatically and data is populated:

vl = QgsVectorLayer("Point", "temporary_points", "memory")
pr = vl.dataProvider()
vl.startEditing()
<... put data into vl>
vl.commitChanges()
QgsMapLayerRegistry.instance().addMapLayer(vl)

New layer is saved as shapefile:

QgsVectorFileWriter.writeAsVectorFormat(vl, "c:\\temp\\my_shapes.shp", "CP1250", None, "ESRI Shapefile")

Then I assume I need  to update my layer to change data type from "memory" to "ogr" with the path to shapefile from step #3.
Then update project file with this info; so the next time this project is loaded - the layer will read data from my_shapes.shp.

I "kind of" know how to make these #4 and #5 steps using QGIS UI:

save layer as shapefile plus set "Add saved layer to map" check box.
delete original layer
save project

Shall I implement these steps programmatically or there is better way?


texts: Raster map & vector boundaryline data for Istanbul, Turkey

texts: I'm looking where I can download background raster map for Istanbul (street map?), and then a boundary line file of the the 39 Districts in Istanbul.  Does anyone have any recommendations?


texts: Creating a Greyscale DEM from Multiband (RGB) GeoTIFF in QGIS

texts:
I'm trying to merge some river bathymetry with LiDAR from the surrounding riverbed using QGIS 2.14.0 Essen. The river bathymetric data is a GeoTIFF, and it comes with a separate scale associating the colors with elevations in feet using the vertical datum NGVD 1929.

The RGB color scheme loads into QGIS without any associated elevations, just values that seem to represent the hues used in each of the 3 bands of data.

If the values in the Band Rendering window corresponded to the elevations associated with each color on the scale included with the bathymetry data (see below) instead of hues, I would simply use the Raster Calculator to add up the 3 bands and divide the sum by 3. The problem is that I have no idea how to associate the color values on the scale with their corresponding elevations.

I am seeking guidance about associating color values with the corresponding elevations on the scale. My ultimate goal is to get this bathymetry data in single-band greyscale like the DEM I intend to merge it with. Once the two sets are merged, I should have a DEM that provides elevation data above and below the river surface.


texts: Is there a service that would return me all streets in a city (ZIP)?

texts: I am looking for a service that you know of, that would return me all streets in a city (ZIP) for the purpose of building a drop-down (or auto-complete) address entry field.


texts: Is it possible to bundle SVG symbology in Geopackage or any QGIS single-file vector format?

texts:
I read that the Geopackage format should support attaching a custom SVG symbology in the file.

Does it already work?

If yes, how can I do this?

No matter if GUI or CLI if it works.

Default point symbols are not enough and advising all the users to install additional ones complicates the situation.


texts: Removing black edges while mosaicking using ArcGIS Desktop or QGIS?

texts:
I want to mosaic three images. Each image has a black back ground. The background is still there even after mosaicking the images, so it's like a gap between the image.

How can I remove this black color using ArcGIS or QGIS.


texts: Google Earth Engine MODIS Pixel-sized Spaced Points Position Adjustment

texts:
I have managed to export a time series of MODIS daily surface temperature pixel data from the Google earth Engine Code Editor using a combination of code scripts provided by Stack Exchange, here is the link https://code.earthengine.google.com/0949dc53da45580dfaf2131250187847. However the pixel-sized spaced points are not all inside the area of interest, I can`t figure out why and I feel maybe my exported data set includes data from outside the area of interest.

I have also managed to create another code script with the pixel-sized spaced points inside the area of interest, here is the link; https://code.earthengine.google.com/f73a5afed5e813f8b83a528303ab6506.
But I can`t manage to combine both codes correctly to export a time series of daily surface temperature pixel data which has all the pixel-sized spaced points inside the area of interest.

Here is the feature collection/function code from the first link which does not have all the pixel-sized spaced points inside the area of interest;

function spacedPoints(fc, proj) {
  var latlon = ee.Image.pixelLonLat().reproject(proj);
  var coords = latlon.select(['longitude', 'latitude'])
                 .reduceRegion({reducer: ee.Reducer.toList(),
                                geometry: fc,
                                scale: proj.nominalScale().toInt()
  });
  var point_list = ee.List(coords.get('longitude')).zip(ee.List(coords.get('latitude')));
  var list = ee.List([0]);
  var feats = ee.FeatureCollection(point_list.map(function(point){
    var ind = point_list.indexOf(point);
    var feat = ee.Feature(ee.Geometry.Point(point_list.get(ind)), {'ID': ind});
    return list.add(feat);
  }).flatten().removeAll([0]));
  return feats;
}

How can I adjust the above code so that all the pixel-sized spaced points are inside the area of interest?


texts: Can a tab masquerade as a space in a unicode string?

texts:
I have some code that is intended to remove spaces from some data I'm ingesting from a group of shapefiles into a number of feature classes.  The malformed data consistently looks like this:

ABC12345678 90

and I want to make it look like:

ABC1234567890

Simple, right?  So the code block goes:

with arcpy.da.UpdateCursor(fc, 'FIELD') as cursor:
    for rows in cursor:
        if rows:
            firstValue = rows
            secondValue = firstValue.replace(' ', '')
            rows = secondValue[:13] #this is because some rows have more than 13 characters
            cursor.updateRow(rows)

I've inserted print statements and it is getting into the loop, but the .replace method is not doing anything.

I've also tried swapping out .replace() with

secondValue = (firstValue[:11] + firstValue[12:])

and no dice there either.

The print statements show:

[u'ABC12345678 90']

through the entire loop - neither replace() nor concatenation is having any effect.

Is it possible that this is not actually a space character?  Could this be a tab value or something which is causing this code to fail?


texts: Saving COG with NaN values using Python and GDAL

texts:
I have a numpy array (data) with some NaN values which I would like to save as a COG using GDAL. I am saving it first to GDAL memory layer and then to COG and I tired using both CreateCopyand Translate as shown below:
# Save to memory
memDriver = gdal.GetDriverByName("MEM")
shape = data.shape
ds = memDriver.Create("MEM", shape[1], shape[0], 1, gdal.GDT_Float32)
ds.SetProjection(proj)
ds.SetGeoTransform(geotransform)
band = ds.GetRasterBand(1)
band.SetNoDataValue(np.nan)
band.WriteArray(data)

# Save to COG
fileFormat = "COG"
driverOpt = []
#driver = gdal.GetDriverByName(fileFormat)
#out_ds = driver.CreateCopy(outPath, ds, True, driverOpt)
out_ds = gdal.Translate(outPath, ds, format=fileFormat, creationOptions=driverOpt,
                        noData=noDataValue, stats=True, outputType=gdal.GDT_Float32)
out_ds = None

When I open the resulting file in QGIS I can see that the NaN value is not ignored when creating pyramids because the gaps expand as I zoom out until the whole image disappears (all pixels have no data value). See the images below. This is the same regardless of whether I use CreateCopy or Translate.
Is the method incorrect or is there some limitation/bug in GDAL COG driver?
I am using Python 3.8.8 and GDAL 3.1.4


texts: Cannot open QGIS form from PostGIS Layer with bigInt datatype for my primary keys?

texts:
Using QGIS 2.10.1 Pisa / PostgreSQL 9.4 with PostGIS.

In PostgreSQL I have column oid with datatype is bigINT as primary key. These are 16-digit numbers.

In QGIS I can read the table as a layer and view it. If i use the Identify tool,
it shows all attributes in the Identify Results Window. In this Window I have set the option auto open form. Normaly if I now use the Identify tool a form opens and show the attributes as a form.

With this table the form does not open. If I change the datatype of the primary key to int, with a 8-digit number, QGIS shows me the form as expected.

I need to use bigInt datatype for my primary keys and I need to show the forms.

Any   suggestions?


texts: Labeling issue after converting database

texts:
I am using ArcMap 10.4.1 and I converted from a .mdb to a file geodatabase. Now the label class has stopped working. The following example is what it used to be, but now I am not able to get it to verify. I have tried to remove the brackets from Label_Class and PIN, but that didn't work. Any suggestions on how to correct?

[Label_Class] = 572 AND NOT (mid( [PIN],5,1) = '6')


texts: Preventing ArcPy from changing MXD extent?

texts:
I want to process mxd datasources without changing the extent.

When I open the MXD in ArcMap and click full extent, I get 1:214,890, after I save and re-open in ArcMap I get what I saved with.

When I open the MXD with arcpy the scale changes to 1:570,396, so after processing and saving via arcpy it's zoomed further out when reopened in ArcMap.

df = arcpy.mapping.ListDataFrames(mxd)[0]
print df.extent


texts: Could not access file "$libdir/postgis-2.3" on Travis installation

texts:
I was searching for a long time to test on Travis my symfony application using postgis extension.

I read the Travis documentation to enabled Postgis and my .travis.yml was this one:

# Project language (I use php)
language: php

# To start postgresql
services:
  - postgresql

php:
  # aliased to a recent 7.0.x version
  - 7.1
  # aliased to a nightly version
  - nightly

# Matrix to test in every php version
matrix:
  # Fast finish allows to set the build as "finished" even if the "allow_failures" matrix elements are not finished yet.
  fast_finish: true
  allow_failures:
    - php: nightly

# Define an environment variable
env:
  - SYMFONY_VERSION="3.3.*" DB=postgresql

# Update composer
before-install:
  - composer self-update

before_script:
  # Testing to reinstall it
  - sudo apt-get install postgresql-9.6-postgis-2.3 postgis -y -q
  # Parameters are copied
  - cp app/config/parameters.yml.travis app/config/parameters.yml
  # Database creation
  - psql -c 'create database symfony;' -U postgres
  # Postgis extension creation
  - psql -U postgres -c "create extension postgis" -d symfony
  # - psql -c 'CREATE EXTENSION postgis_topology;' -U postgres -d symfony
  # excution de composer
  - composer install
  # Database migration to current version
  - php bin/console doctrine:migrations:migrate  --no-interaction
  # loading data for phpunit and codecept tests
  - php bin/console doctrine:fixtures:load --fixtures ./src/AppBundle/DataFixtures -n --env=test

script:
  # Tests de versions
  - psql --version
  - psql -d symfony -c 'SELECT PostGIS_version();' -U postgres
  # Phpunit
  - if [ "$TRAVIS_PHP_VERSION" == "7.1" ]; then travis_wait 40 php vendor/phpunit/phpunit/phpunit --coverage-clover build/logs/clover.xml; fi
  - if [ "$TRAVIS_PHP_VERSION" != "7.1" ]; then php vendor/phpunit/phpunit/phpunit; fi

But an error occured when creating postgis extension. Error occured on this line:

$ psql -U postgres -c "create extension postgis" -d symfony

  ERROR:  could not access file "$libdir/postgis-2.3": No such file or directory


texts: How to solve issue with log messages panel in QGIS: "Not logging more than 100 request errors."?

texts:
I am using the following code in pyqgis to catch errors/warning from a WMS-layer, in order to trigger a repaint as soon as an error/warning is detected (based on previous question: How to catch WMS error message from log messages panel in QGIS with python?)

But obviously the "WMS" provider seems to have a restriction of not sending more than 100 error requests to the messages log, meaning after the 100th error/warning I´m no longer able to catch any signal, even if the WMS-layer is still not responding correctly.
Nevertheless, if I send own messages to the log panel there does not seem to be any restriction (see code below).

Is there a possibility to catch the error/warning directly from the instance that is responsible here (I guess it is the WMS-provider), instead of using the messages log panel? Or maybe just clear/reset the log messages panel in a running process or remove the limitation?

I am using QGIS 2.18.2 on Windows 10.

Here is the python code:

# coding=utf-8

from qgis.core import *

wmsLayer_name="wms-dtk50_wgs"
url_with_params ='url=http://sg.geodatenzentrum.de/wms_dtk50?&crs=EPSG:25832&featureCount=10&format=image/png&layers=DTK50&styles='

wmsLayer = QgsRasterLayer(url_with_params, wmsLayer_name,'wms')
QgsMapLayerRegistry.instance().addMapLayer(wmsLayer)

def errorCatcher( msg, tag, level ):
    if tag == 'WMS' and level != 0:     #Warnings or Errors (0: Info, 1:Warning, 2:Error)
        print "WMS error detected!"
        myWMSLayer = QgsMapLayerRegistry.instance().mapLayersByName("wms-dtk50_wgs")[0]
        myWMSLayer.triggerRepaint()

# connect with messageReceived SIGNAL from QgsMessageLog to an errorCatcher custom function
# instantly reacts if error/warning occurs
QgsMessageLog.instance().messageReceived.connect( errorCatcher )

#after 100 times triggering a "wmsLayer.triggerRepaint()",
# I get following warning in log messages panel "WMS":
# "2017-01-17T07:17:52  1   Not logging more than 100 request errors."

#this does not raise any issues and prints all 500 test messages in the log panel:
for i in range(500):
    QgsMessageLog.instance().logMessage("Message #{}".format(i),"Test",2)

UPDATE: I submitted a feature request (see: https://hub.qgis.org/issues/16168)


texts:
I downloaded gridded soils data from gSSURGO and used ArcGIS to export the raster to a GeoTiff of keys for use on Linux.  On my Ubuntu machine, I have imported the GeoTiff to QGIS as a raster layer, and have imported the geodatabase table as a vector, using the instructions here File Geodatabase (*.gdb) support in QGIS?, but they are still separate from one another.

How do I connect the GeoTiff with the geodatabase so that I can visualize values from the geodatabase?


texts: postgis under win 7 "on demand" not as a service

texts: how do I install postgis2.0  under win 7 to run just on demand. I dont want to run it permanently as a service. Similary as the (nonspatial) mysql can be run on demand via the xampp-tools (http://www.apachefriends.org/de/xampp.html)


texts: How to render a transparent raster layer on top of one or more vector layers with PyQgis?

texts:
I been trying to render a map with a 3 vector layers and one raster layer. The idea is that the raster layer should be on top of the other 3 vector layers. The problem is that even when I add the transparency parameters to the Raster Layer I can't visualize the vector layers. In other words, the Raster Layer completely covers the area.

#!/usr/bin/python
# -*- coding: utf-8 -*-
from qgis.core import *
from qgis.gui import *
from qgis.utils import *

from PyQt4.QtCore import *
from PyQt4.QtGui import *
from PyQt4.QtCore import QFileInfo, QSettings

QgsApplication.setPrefixPath("/usr", True)

QgsApplication.initQgis()

# raster layer
filename = "points.png"
file_info = QFileInfo(filename)
base_name = file_info.baseName()

l_raster = QgsRasterLayer(filename, base_name)
l_raster.setCrs( QgsCoordinateReferenceSystem(4326, QgsCoordinateReferenceSystem.EpsgCrsId) )
l_raster.setTransparency(20)

l_a = QgsVectorLayer("a.MIF", "Layer A", "ogr")
l_b = QgsVectorLayer("b.MIF", "Layer B", "ogr")
l_c  = QgsVectorLayer("c.MIF", "Layer C", "ogr")

if not l_raster.isValid():
  print "Raster Layer failed to load!"

if not l_a.isValid():
  print "A Layer failed to load!"

if not l_b.isValid():
  print "B Layer failed to load!"

if not l_c.isValid():
  print "C Layer failed to load!"

QgsMapLayerRegistry.instance().addMapLayer(l_raster)
QgsMapLayerRegistry.instance().addMapLayer(l_a)
QgsMapLayerRegistry.instance().addMapLayer(l_b)
QgsMapLayerRegistry.instance().addMapLayer(l_c)

# create image
img = QImage(QSize(800,600), QImage.Format_ARGB32_Premultiplied)

# set image's background color
color = QColor(255,255,255)
img.fill(color.rgb())

# create painter
p = QPainter()
p.begin(img)
p.setRenderHint(QPainter.Antialiasing)

render = QgsMapRenderer()

# set layer set
# add ID of every layer
lst = [l_raster.getLayerID(), l_a.getLayerID(),  l_b.getLayerID(), l_c.getLayerID()]
render.setLayerSet(lst)

# set extent
rect = QgsRectangle(render.fullExtent())
rect.scale(1)
render.setExtent(rect)

# set output size
render.setOutputSize(img.size(), img.logicalDpiX())

# do the rendering
render.render(p)

p.end()

# save image
img.save("render.png","png")

QgsApplication.exitQgis()


texts: ArcGIS - Custom, rule-based borders for polygons?

texts:
I've been working with a map that looks something like this:

Wherever two polygons of different types (shown above with different colors) meet there needs to be a thick dashed line between them. The solution so far has been to keep a separate feature class for the lines, which of course needs to be updated every time the polygons change. This has gotten to be pretty tedious and is creating a lot of redundancy in our data. Is there any type of map layer that could render the lines on the fly based on spatial relationships and attribute data? I'm guessing a geometric network would be involved - does anyone have any insights?

EDIT: Sorry everybody, I made a mistake with the screenshot and accidentally had another layer on. The polys are already dissolved to begin with. Please see above for a new version.


texts: How to upgrade to 10.5.1 using Azure Cloud Builder

texts: I installed ArcGIS Enterprise 10.5 in Azure using the Azure Cloud Builder, now I need to upgrade my portal to 10.5.1. How can I upgrade to 10.5.1 using Cloud Builder or do I have to manually upgrade?


texts: EPSG Geodetic Parameter Registry - SW supporting transformation EPSGs

texts:
The page www.epsg-registry.org contains EPSG codes for different transformations (i.e. EPSG:1173, EPSG:1244 etc.).

Do you know about some GIS SW (online or desktop), which supports converting coordinates using these transformations?


texts: GEE: mapping over a wind data (ERA5) image collection

texts:
I want to extract average daily wind speeds over the Netherlands (NL below) at 1pm from the ERA 5 dataset on GEE.
I am trying to map over an image collection, and for some reason it is telling me that "u.pow is not a function". The formula is right (I tested it on a single image).
It's the calcwind variable that doesn't run as a result. If I understand correctly, this code should add a windspeed attribute to the dataset which I will then use to identify days with low wind speed.
I will compare this wind data to S5P NO2 data, which should contain u and v wind vectors (see Panoply screenshot below), but unfortunately this seems to be excluded in the GEE L3 NO2 dataset...
//Setup

var countries = ee.FeatureCollection("FAO/GAUL/2015/level0");
var NL = countries.filterMetadata('ADM0_NAME', 'equals', 'Netherlands');

var ERA5 = ee.ImageCollection("ECMWF/ERA5_LAND/HOURLY")
     .filter(ee.Filter.date('2020-07-01', '2020-07-02'))
     .filterMetadata('hour', 'equals', 13);

var ERA5NL=ERA5.map(function(im){
   return im.clip(NL)
})

// Wind calculation
var u = ERA5NL.select('u_component_of_wind_10m');
var v = ERA5NL.select('v_component_of_wind_10m');

var calcwind = function(feature) {
    return feature.set({windspeed:((u.pow(2)).add(v.pow(2)).sqrt()) })
};

var ERA5x = ERA5NL.map(calcwind);

//Mapping
var vis = {
  min: 0,
  max: 5,
  palette: [
    '#FFFFFF', '#FFFF71', '#DEFF00', '#9EFF00', '#77B038', '#007E55', '#005F51',
    '#004B51', '#013A7B', '#023AAD'
  ]
};

Map.setCenter(5.5, 52, 7.5);
Map.addLayer(ERA5x.select('windspeed'), vis, "Wind speed (m/s) 10m above surface");


texts: OpenLayers: Properly use a WMS layer to display a TMS layer from a different server

texts:
I have two servers:
the first server is 177 and that has a html file with the TMS layer.
We will call this tms_layer.html and it has your init function with this in there:
var tmsoverlay = new OpenLayers.Layer.TMS( "TMS Overlay", "",
  {
    serviceVersion: '.', layername: '.', alpha: true, type: 'JPG', getURL:
    overlay_getTileURL, isBaseLayer: true
  });

When i go to 177/tms_layer.html the tiles are being correctly displayed.
So that is working fine.
Now, I want to serve that tms_layer.html file on a different server (server 190), and I have been trying to serve it as a WMS layer. Not sure if WMS is the best approach but I think it is. Correct me if I'm wrong.
We'll call this tiles_revision1.html and within the init function I currently have:
layer = new OpenLayers.Layer.WMS( "My TMS layer in WMS",
  "http://177/dev/tiles/tiles_1a.html", {layers: ['tmsoverlay'], transparent: "true", format: "image/jpg"} );
map.addLayer(layer);

For the layers: [] do I even need to reference this since it's just a bunch of jpgs making up the tile(s)/map? Leaving it blank has not done anything.
With the above code I currently get a blank OL map, see here (When I pan around the map I get a Uncaught TypeError: cannot read property 'wrapDateLine' of null. I don't think that is the root of my issue though):

Edit: the most important thing here is to be able to access those 177 tiles from anywhere as a baselayer. There is quite a few associated jpgs, and the tms_layer.html file is pretty beefy on its own.
What is the correct way to do this?


texts:
How can I import a shapefile in python?

from PyQt4.QtCore import QSettings, QTranslator, qVersion, QCoreApplication
from PyQt4.QtGui import QAction, QIcon, QFileDialog
# Initialize Qt resources from file resources.py
import resources_rc
# Import the code for the dialog
from Import_tool_dialog import ImporttoolDialog
from functools import partial as pto
from Tkinter import *
from tkMessageBox import showinfo, showwarning, showerror
from xlrd import open_workbook
import os
import shutil
import time
import psycopg2
import shapefile as shp

Here I am getting error :

ImportError: No module named shapefile.


texts: Writing Python toolbox that gives two options for one parameter?

texts:
I'm working on a Python toolbox that buffers features using slope distance based on an input DEM rather than horizontal distance. The tool works great with buffer widths derived from a field in the input feature.

I'd like to add the option to buffer by one given width like the ArcGIS proximity tool, Buffer, resulting in a radio button choice:


texts:
I try to create an application for taxis. The application will be GMap based.

On the map will be marked the region in which the taxi driver is working in.

The region will be a polygon designed by the application owner.

Also the client charge will be based on kilometer.

The end user can define the departure and destination points on the map. In case that both locations are inside the region all are fine.

In case the destination point is outside the region, the client must be charged in double rate for the kilometers outside the region.

What I like to ask is the following:

Is there a way to know what is the distance from the point inside the region until the region border?
And what is the distance from the point outside the region until the region border?

To be more specific, I attach an image here:

What I need is to get the distance from departure location to the border of the region (the yellow circle) and from yellow circle to the destination location.

I know that this is very complex, but I don't even know if is possible.

Can somebody to help me with that problem ?


texts: Converting styles from ArcGIS Pro (.stylx) to ArcMap (.style)?

texts:
I want to open .stylx file from ArcGIS Pro in ArcMap.

Is there any tool in ArcGIS Pro to convert .stylx files to .style format?


texts: QGIS 3.1 Producing a proportional circles legend - No option for "Data defined size legend"

texts:
I am making a global map in qGIS 3.1 using proportional symbology with the following:

Data = Number of hospitalizations - classified as 'Real' type
Symbology = Single Symbol: Fill -> Centroid fill -> Marker -> Simple Marker

For Size under 'Simple Marker' I have checked the assistant tool and selected my hospitalizations variable as the source.
Unfortunately I can't figure out how to produce a legend for this. All solutions I'm seeing direct people back to 'Marker' and then 'Advanced' for 'Data defined size legend' but I do not have this option, only says Symbol Levels.


texts: Hillshade tool in ArcMap returns terraced, pixelated image

texts:
I am applying Hillshade tool (ArcMap 10.1) to SRTM in Caucasus Mountains (border between Russia and Georgia). It gives me weird image with 'terraces' (see attachment). I reprojected SRTM from WGS-1984 to UTM 37 zone using cubic convolution, but it didn't help. Also, I tried to apply Z-factor 0.00001171 for 40 degr. like it says here : http://blogs.esri.com/esri/arcgis/2007/06/12/setting-the-z-factor-parameter-correctly/, but it says 'The numeric value is outside of the valid range.'
Any suggestions on how to fix it would be very much appreciated.


texts: mod-tile: How to set up two independent tileservers?

texts:
Reading mod-tile and rendered documentation, especially the commented lines in the renderd.conf file, it seems to be possible, with the same renderd and mod-tile, to define two different style sheets that query two different databases and write tiles to two different folders.

But I couldn't find exactly how this worked. Does someone have an idea or link explaining this?


texts:
I've been trying to merge some raster tiles at work and I've been getting an error message saying I don't have the correct permissions. This works fine on my home pc.

Any obvious reasons for this?

QGIS at work runs on an application virtualisation client. Also 3rd party repos won't work. Set proxy settings to allow standard QGIS repo and that works fine.


texts:
whenever i try to add a buffer to a certain layer it does that in degrees , i realized that i need to change the CRS but there are a ton of them
what do i need to know to choose the correct CRS for that layer to change it into meters?

The area that I am working on is Baghdad in Iraq.


texts:
I currently rely on a mixture of mailing lists, twitter and a couple of blogs.

What do you use for your GIS media intake, and why?


texts: Closing polygon with right click on Mac OS using QGIS?

texts:
I frequently use QGIS to digitize satellite imagery or georeferenced maps. Recently, however, the Add Polygon Features and Add Line Features have stopped working - I can Add a Polygon Feature and trace the polygon, but QGIS will not respond to my right click which should finalize the trace and initialize the polygon feature.

Right-clicking (using control + , as this is a Macbook OS 10.13.6) works fine for other operations in QGIS, such as right clicking on a Layer to export, etc.

I have tried QGIS versions 3.4, 3.8, 3.10. And I have used every conceivable keyboard shortcut to see if it would register as a right click (command click, shift click, . click, etc). I have tried to find a way to change how QGIS interprets the keyboard, but I haven't found anything.

What's more, after trying long enough to 'close' the polygon, QGIS stops responding to the Esc key to exit polygon creation mode.

If anyone has encountered this issue, how did you resolve it?


texts: What is the width of the beach of a little island polygon?

texts:
This is a "classic" problem, that have no practical utility (see note-1 for applications), but illustrate another GIS problems...

We can do some approximations to modeling the situation,

the real center (coconut) is not at the ST_Centroid point, but can be approximated to;
the real shape (a island) is not a circle, neither a square, they are optional references.

So, translating the question to an equation, and using OGC standard functions... "What is the value of w?" at the equation

ST_Area(ST_Buffer(ST_Centroid(littleIsland), w, quad_segs)) = ST_Area(littleIsland)

When quad_segs=1, the reference is a square, when quad_segs>=8 the reference is near to a circle (see examples of generated polygons).

NOTE-1: For understand this problem into a context of applications,
see "Is there a st_buffer inverse function, that returns a width estimation?" and
see "How to fit a polygon to a “best buffer” about a reference geometry in it?".

NOTE-2: the "inverse function" is so general and complex, because can operates with point, line and polygon inputs. Here we can discuss the simplest case, when the input is like a point, with no area or length.


texts: How to export vector features to grid using proxy and WFS protocol?

texts:
I am trying to export vector features (current view) as a grid, something like: http://workshops.opengeo.org/geoext/wfs/grid.html

so far layer features are not being retrieved.. I got empty grid in the window. would anyone take a look and let me know if I miss something ?

gridPanel = new Ext.grid.GridPanel({
        title: "Feature Grid",
        store: new GeoExt.data.FeatureStore({
        fields: [
            {name: "id", type: "int"},
            {name: "name", type: "string"}
        ],
            proxy: new GeoExt.data.proxy.Protocol({
            protocol: new OpenLayers.Protocol.WFS({
                url: "http://IP/geoserver/wfs",
                featureType: "myLayer",
                featureNS: "http://mapmap.org",
                srsName: "EPSG:3857",
                version: "1.1.0"
            })
        }),
                autoLoad: true
                }),
                sm: new GeoExt.selection.FeatureModel(),
        width: 320,
    columns: [
        {header: "id", dataIndex: "id"},
        {header: "name", dataIndex: "name"}
    ]
    });

        //window to display grid
        var gridWindow = new Ext.Window({
              items: [
                     gridPanel
                       ]
            }).show();


texts: How to convert WMS layers to raster layers QGIS 1.8.0

texts: I need to export a series of maps for a report. The layers imported into QGIS are all WMS files provided by the GeoBC website. As I have now discovered, I cannot export these maps using the print composer because of the limitations of WMS. I also cannot add a scale bar to these maps as the WMS CRS is WGS 84 not UTM. If I understand the problem correctly, the solution would be to convert the WMS layers in to raster layers initially, and then to vector layers? Does this sound correct? If so, is there a tool to convert WMS to raster I can use as the first step?


texts: Simplest way to share map and data with ArcGIS 10.x and 9.x users

texts:
What is the simplest way to share a map and data with both ArcGIS 10.0 and 9.3 users, portably?

This question could be rephrased as "Instead of waiting for Esri to (maybe) implement backwards compatible map packages, what can we do?"

Assume that one is starting with project resembling the following, all
in v10:

D:\workspace:
    MyCoolMap.mxd:
        default.gdb:
            Transportation\roads
            Transportation\trails
            Hydro\Lakes
            Hydro\Rivers
            some_random_points
            a_linked_table
      a_shapefile.shp
      a_linked_table_for_fun.csv
      {some_coded_value_domains}

I've tried saving the mxd as 9.3, and then create a map package from that, but get an "error compressing map package".

Using the Feature Class To Geodatabase (Multiple) tool is unsatisfactory because it can't handle Feature Datasets.

Another failed approach was to export the v10.gdb to binary xml and import it into a v93.gdb, ("Import xml data failed. The domain field type does not match that of the field it is being assigned to [feature class name]")

Got any other approaches that might bear fruit?


texts: Updating Description of MXD Layers after metadata updated on source using ArcCatalog?

texts:
Is there a way to force the update of the Description content in the Properties/General/Description of each feature layer in an MXD?

I frequently update metadata for each feature class in ArcCatalog, but these changes are not coming through to the feature layer in MXDs.

Using ArcGIS 10.2.1


texts: "add feature" tool QGIS 2.0 - polygon fill blocks view of what I want to digitize

texts:
How can I change the settings on the vector "add feature" tool so the tool fill doesn't block my view of what I want to digitize?  This problem popped up after I upgraded to QGIS 2.0

QGIS 2.0
Macbook Pro  OS 10.9


texts:
How do I invert the color ramp for an raster data in QGIS 3.0?

I tried to find the invert menu in single pseudocolor render type in QGIS 3.0 but I don't see it there.

I was using QGIS 3.0.3


texts: FieldPyculator Plugin - How to Classify Fields?

texts:
Could someone tell me how to classify/populate fields with QGIS FieldPyculator Plugin?

In the Attribute Table, I have a field named POPULACAO with population data:

I want to apply a weight to classify these fields. I've installed QGIS FieldPyculator Plugin via repository. Below the Python code to complete this task in ArcGIS 10:

def reclass(POP):

if (POP <= 10000):
return 1
elif (POP > 10001 and POP <= 30000):
return 2
elif (POP > 30001 and POP <= 50000):
return 3
elif (POP > 50001 and POP <= 100000):
return 4
elif (POP > 100001 and POP <= 200000):
return 5
elif (POP > 200001 and POP <= 500000):
return 6
elif (POP > 500001 and POP <= 1000000):
return 7
else:
return 8

In the box below:

reclass( !POPULACAO! )

The ArcGIS 10 Python code works like a charm.
Data shapefile: http://bit.ly/xBW8Pz | In FieldPyculator, I did wrong, I know:

Please help me. Jorge Santos


texts: Clipping line feature to MANY polygons for measurements

texts:
I am looking to measure the length of stream segments(line) on numerous tax lot properties (Polygon).

Is there a tool that will allow me to "Cut" or "Clip" streams to tax lots to estimate not only length, but to create a buffer for each property and calculate area?


texts: Issue with exporting GRASS visibility graph to shapefile

texts:
I am rerunning code in GRASS that has worked in the past (about six months ago) and am now having problems.

In the GUI, I have created a new mapset and imported the 'SimplifedCoastline' file to the location (and set the default region based of the coastline's extent) before entering the code below into the Python Shell within GRASS. The objective is to create a visibility graph of shortest paths/lines between the points (centroids.txt) while going around the coastline as an obstacle.

The code:

import os
import grass.script as gscript

directory = r"E:\CollabProject\VisibilityAnalysis\VisibilityInputs\IntegerInputs"
os.listdir(directory)
centroids = open(os.path.join(directory, "Centroids.txt"), 'r')
data = centroids.read()
vn = "Centroids_visibility"

gscript.run_command("v.net.visibility",input="SimplifiedCoastline", output=vn, coordinate=data)

outname = "E:\CollabProject\VisibilityAnalysis\VisibilityOutputs"
gscript.run_command("v.out.ogr", input=vn, type="line", dsn=outname, output="visibility_output.shp")

The v.out.ogr command results in:

Traceback (most recent call last):
File "<input>", line 1, in <module>
File "C:\OSGEO4~1\apps\grass\grass-7.2.2\etc\python\grass\script\core.py", line 410, in run_command
return handle_errors(returncode, returncode, args, kwargs)
File "C:\OSGEO4~1\apps\grass\grass-7.2.2\etc\python\grass\script\core.py", line 329, in handle_errors
returncode=returncode)
CalledModuleError: Module run None ['v.out.ogr', 'dsn=E:\\CollabProject\\VisibilityAnalysis\\VisibilityOutputs', 'input=Centroids_visibility', 'type=line', 'output=visibility_output.shp'] ended with error
Process ended with non-zero return code 1. See errors in the (error) output.

I'm not sure if it is failing because of an issue with v.out.ogr or perhaps because v.net.visibility isn't running correctly and creating an empty output, so therefore v.out.ogr can't run properly?

I have tried:

running this code on several other computers, as some forum posts suggested this error could result from GDAL conflicts from downloading through OSGeo
downloading GRASS alone onto a computer and running the code; downloading GRASS and QGIS onto a computer and running the code
importing the data from a shapefile rather than text file
running the code on a test dataset with only 27 points so it shouldn't be maxing out data limits
going through the documentation for v.net.visibility and v.out.ogr and searching on forums

I seek any suggestions on how to correct this, as well as view the output of v.net.visibility to confirm it is being created correctly.

I can view the visibility graph that is being created in GRASS' Map Display so I know it exists, but I can't export it using v.out.ogr. The error I receive when trying to export it through the Python Shell is above; the error I receive when trying to export it through the point and click method is:

v.out.ogr input=inputdatatestttt@PERMANENT output=E:\CollabProject\VisibilityAnalysis\VisibilityOutputs\visgraphout format=ESRI_Shapefile
WARNING: No attribute table found -> using only category numbers as attributes
Exporting 2470954 features...
WARNING: 2470954 features without category were skipped. Features without category are written only when -c flag is given.
WARNING: Output layer is empty, no features written
v.out.ogr complete. 0 features (Line String type) written to <inputdatatestttt> (ESRI_Shapefile format).

I'm not sure why the visibility graph is being created without attributes or how to fix this.


texts:
Once in the good old times of QGIS 2.14.17, there was an 'extract by location' algorithm:

In QGIS 2.99.0 neither searching for 'extract' nor for 'location' delivers such result:

(Oops, NO hit when searching for 'location' in a geospatial processing environment?!)

Question: How to 'Extract by location' in QGIS 2.99.0?


texts:
Is there a way to recognize a Tiled Internet Layer using arcobjects in .NET?
I've tryed several ways but I always end up in a dead road...

ILayerFile lyrf = new LayerFile();
lyrf.Open(@"G:\mypath.lyr");
ILayer layer = lyrf.Layer;

MapServerLayer msl = layer as MapServerLayer;

ITiledMapServer tms = msl as ITiledMapServer;

I'm stuck obtaining the MapServerLayer: layer doesn't answer to ITiledMapServer nor to IMapServerLayer.
Is there another way to obtain it?


texts:
I have a basic question regarding Geo data: is the normal assumption mandatory?
Given a dataset I identified outliers and inliers (using Moran local and plot test) so I excluded them from the original data. I checked normality in this new dataset, however it does not fit a normal distribution.
What is the effect of this non-normality in the variogram and krigging?


texts: GeoPandas legend not displaying correctly the colors

texts:
I need a legend where for every bin that I defined there is a color from a colorscale.
I wrote this:
ax= merge.plot(column=dates,
               figsize=[15,15],
               cmap='OrRd',
               scheme='user_defined',
               classification_kwds={'bins':[10,20,50,100,200,500,1000]},
               legend=True
               )

and I keep getting this:

As you can see the first three colours are white and the last three are dark red.
How can I modify my code to get a shade of colour for a single bin?
Here's an example of what I would like it to be:


texts: Layer styling with feature 1 but atlas pages generated by feature 2

texts:
How can I tell the atlas-composer to display on one page only all geometries with feature two and all the other geometries should be hidden?

I have four geometries with two feature-attributes: (1) "Kategorie" and (2) "Merkmal" in the layer with the name "Test".
The layer-style in the map is made with categorised by feature one ("Kategorie").
The atlas-composer should display on one page only all the geometries which compare the same value in feature two ("Merkmale") and all the other geometries should be hidden.

I think I can write some SQL in the atlas-parameters but I'm not that experienced. When I use e.g., dissolve() which generates a new layer, then I am not able to hide the other geometries. Any further help is welcome.


texts: How to convert csv file into geo database file?

texts:
I have a csv file that i want to export to a file geodatabase (*.gdb) in C#.

I am currently looking into GDAL ogr2ogr commands.

They have a way to convert gdb file into csv file but not the opposite.

Is there a way using GDAL or ogr2ogr to export csv to gdb file?

What i really want to do is that my database table's data is already exported into separate csv files and i want to write that data from csv files to FileGDB and export it (the entire directory as a zip). GDAL ogr2ogr has commands to export FileGDB into csv or several formats but not the opposite. So if i have a ogr2ogr command or a way to do that all i will do is run the cmd command from c# code.

I know how to run the command from a code, all i need is a command to export csv to FileGDB. (it doesn't have to be a command, it can be a external driver or tool or configuration or a method to do that)


texts: ArcGIS Javascript API - carrying out multiple query tasks

texts:
I have carried out a QueryTask on a map service that returns a polygon when I input a postcode. Within the same service there is an additional layer that allows the user to query a local authority name on.  I would like to carry out a querytask on this layer at the same time so that the map returns a polygon relating to the postcode input and also tells me what local authority it is contained within. I will ultimately show the local authority in a separate widget.
Can anyone help please?


texts:
I need to select from a layer(s) and I am wondering what way I should serve the layer.

The requirements are that I need to be able to select multiple features at a time by drawing a box or polygon. They should become highlighted and I should display the data for all the features. I don't need to do any editing of them.

I'm using geoserver and openlayers


texts: Values calculated by pyproj and the ones by a website differ greatly?

texts:
By this website http://www.movable-type.co.uk/scripts/latlong.html,

I would like to calculate distances between latitude 37.919 and 37.465, longitude 15.358 and 14.601, the results are 50.48 km and 66.81 km, respectively.

e.g. Point1: 37.919, 14.601 and Point2: 37.465, 14.601, the website gives 50.48 km,

Point1: 37.465, 15.358 and Point2: 37.465, 14.601, I got the distance among them, 66.81 km.

However, I use the pyproj.Proj to calculate the distances, which differs greatly from the results above.

Below is my python code:

import pyproj.Proj
p=Proj(proj='eqc',lat_0=37.692, lon_0=14.980)
p(14.601,37.465) #return (-42190.08701065082, -25269.524410073795)
p(15.358,37.919) #return (42078.76751985713, 25269.52441007167)

The distance between 37.919 and 37.465 is 25269.52441007167-(-25269.524410073795), roughly 50539 m = 50.539 km, which matches approximately the value by the website: 50.44 km

But, the distance between 15.358 and 14.601 is roughly 84269 m, 84.269 km, which differs greatly from the one by the website: 66.81 km.

I don't know why?


texts: Adding styled OpenLayers Google map in QGIS

texts:
I came across a great tutorial about how to add a styled Google map layer to the QGIS OpenLayers plugin (http://gistncase.blogspot.com/2014/02/how-to-add-google-maps-styled-layer-in.html).

However, it appears that some changes have been implemented on the QGIS/OL side since this was posted and it no long works.

I have sorted out how to have openlayers_plugin.py find and list the new custom layer, but I am unable to get the layer html itself to function properly. (Likely due to my complete lack of experience on the html/java side).

QGIS: 2.10.1-Pisa
OpenLayers: 1.3.6

Tutorial example (that doesn't seem to be working):

    <title>OpenLayers Google Streets Layer</title>
    <link rel="stylesheet" href="qgis.css" type="text/css">
    <link rel="stylesheet" href="google.css" type="text/css">


texts: Converting points to Percentage Volume Contours using ArcGIS 10.1 for Desktop?

texts:
I need to convert about 100 different sets of points into percentage volume contours. These are currently in vector format.

What would be the most efficient way to do this? I'd rather do it in ArcGIS but free tools / software is also an option, provided they are reasonably easy to use (I'm not too familiar with coding / prefer user-friendly interfaces).

I'm hoping there's something out there I've missed - so far my only option is to do this in GME after converting the points to raster within ARC. Even with batch conversion it will take me a while to insert each of the files in, one by one (plus I don't fully understand what some of the conversion options actually involve).

As a side problem, I also have in-group overlapping points (points in the same location) which I'd prefer to include in the PVC calculation as they are.

Any help will be much appreciated -

Thanks

EDIT: Do see your point about explaining a little more. Basically I'm looking at naming conventions for different places. For each name, there are a bunch of points which correspond to where people perceive this place to be. Ideally: assuming this follows a normal distribution in space, I'd like to do kernel density estimates of these (gaussian) and then do contour lines which reflected the % distribution (so 95%, 50% and a couple of others maybe). It's exploratory, I can't code, but I'm looking to see if I can do some contours of a bunch of different places and where people think these places are - a sort of vernacular gazetteer. I have at least 100 places though, & each of these should have their own set of contours. Hawth's KDE tool seemed perfect for this as it did both gaussian KDE and PVC's at the same time in batches. Unfortunately it won't work in ArcGIS 10.1 and I can't afford the days it would take me to figure out GME at the moment (would like to though). So I'm a little stuck & need an efficient way to do this.


texts: Adding alternatives to least cloudy image to a collection

texts:
I've made a collection of least cloudy images based on the EE tutorial guide, using Landsat 5 data. However, for the area I'm interested in, some of the returned images are not fully covering the area. I'm using the code below:

var imageFin1984 = ee.Image(
  l5.filterBounds(ROI)
    .filterDate('1984-07-01', '1984-08-31')
    .sort('CLOUD_COVER')
    .first()
    ).clip(ROI);

var visParams = {bands: ['B3', 'B2', 'B1'], max: 0.3};
Map.addLayer(imageFin1984, visParams, 'True Colour 1984');

How would I view the second least cloudy image, or generally improve the coverage of the image for my area?


texts: Get Lon Lat values from a GEOGRAPHY data type in PostGIS

texts:
After storing lon lat values as a GEOGRAPHY data type, how do I retrieve the individual lat lon values?

Failed attempt:

SELECT id, geog, ST_X(geog), ST_Y(geog) FROM locations;

Error:

No function matches the given name and argument types. You might need to add explicit type casts.


texts: Delimiting regions of influence for ports

texts:
I have datasets containing origin-destin (ports) and volumes for trade data.

From these datasets, I can extract points of origins that are related to only a particular port. Interpolating such data, I can generate the images below, which show the "area of influence" of a particular port.

Now, I want to merge all of these rasters, in order to create some sort of what I believe would look like a series of "watersheds" delimiting the area of influence for the different ports

How can I "merge" these different areas in a single map, with all the multiple regions? Is it possible?


texts: I've found there's integration of DKAN into cartodb http://www.nucivic.com/publish-beautiful-maps-cartodb-dkan/ but does it really allow to setup a separate webpage with a cartodb map?


texts: CopyFeatures management creates long name

texts:
I'm trying to join a shapefile to a csv file. The join works but the attribute table has the long name of the shapefiles which makes the fields not readable.

How can i fix this problem?

Calling the scrip from an os window as

C:/Python27/ArcGIS10.4/python csv_2_dbf_tdm.py "D:/01_Inputs/y2035/cube_se2035.csv" "no_bld_tdm_draft1"

where cube_se2035 is sys.argv[1] and the no_b;dm_tdm_draf1 is a folder name

import arcpy
import sys
import os

basename =  os.path.basename(sys.argv[1]).replace(".csv", ".dbf")
dbf_path = sys.argv[1].replace(".csv", ".dbf")
se_shp = os.path.join(os.getcwd(), "data\se_folder_shapefile\se919.shp")
sys.argv[1] = se_path = "D:/cube_se{year}.csv".format(year=2035)

shapefileName = "temp"
arcpy.MakeFeatureLayer_management(se_shp,shapefileName)

###Join the feature layer to a table
arcpy.AddJoin_management(shapefileName, "DASZ2010",dbf_path, "dasz_id")

### Copy the layer to a new permanent feature class

arcpy.CopyFeatures_management(dd, os.path.join(cube_folder, basename.strip(".dbf")))

Here is the table after the join. All the fields have the source name and I just need the field name without the source.


texts: PYQGIS error in file format working in standalone script using PCharm

texts:
i want to run this PYQGIS python script in standalone out of QGIS and i use PyCharm Idle.
but my code have error on the file format.

the script

from qgis.core import *
import sys
import qgis.utils
from osgeo import ogr
from PyQt4.QtCore import *
from PyQt4.QtGui import *
from qgis.gui import *

# supply path to qgis install location
QgsApplication.setPrefixPath("/usr", True)#for linux

# create a reference to the QgsApplication, setting the
# second argument to False disables the GUI
qgs = QgsApplication([], False)
# load providers
qgs.initQgis()

# Write your code here to load some layers, use processing algorithms, etc.
sys.path.append('/usr/share/qgis/python/plugins')
from processing.core.Processing import Processing
import processing
Processing.initialize()
rasterlyr=QgsRasterLayer('/home/qgis2/PycharmProjects/test/aster.tif')
rasterlyr.isValid()
#slope='/home/qgis2/PycharmProjects/test/slope.tif'
#aspect='/home/qgis2/PycharmProjects/test/aspect.tif'
processing.runalg('saga:slopeaspectcurvature', rasterlyr,6,1,1,'/home/qgis2/PycharmProjects/test/slope.tif','/home/qgis2/PycharmProjects/test/aspect.tif',None,None,None,None,None,None,None,None,None,None)
#processing.alghelp("saga:slopeaspectcurvature")
# When your script is complete, call exitQgis() to remove the provider and
# layer registries from memory
qgs.exitQgis()

error message :

ERROR 4: `/tmp/processing969c81d9be194fa49b8dd2c0daebf68b/f9a3240b743f4554b2602e7940bf8fd2/SLOPE.sdat' not recognised as a supported file format.

*** Error in `/usr/bin/python2.7': corrupted double-linked list: 0x0af99d30 ***
======= Backtrace: =========

backtrace have many linux paths.


texts: Python Script for Select by Attributes taking user input

texts: I am trying to write a python script and create a script tool in ArcGIS that will let a user choose the value or values that will be used to extract a shapefile. These selected value(s) will then need to pass through the select tool. I'm sure this involves writing a WHERE clause but I'm new to python and am a little lost on this one.


texts:
I have a series of UTM coordinates in an Excel file that I need to convert to Lat/Long.  They are listed in the spreadsheet as "LONG(XUTM)" and "LAT(YUTM)".  I know the coordinates are located in Bulgaria.  Trying to figure out exactly how I can convert these - so far I'm using tools online and I'm getting lat/longs that put me in the ocean or in sudan or something...

Here are a selection of the coords:

LONG(XUTM)  LAT(YUTM)
5400054.968 2262690.795
5410764.956 2279927.441
5406759.683 2282320.619
5402952.243 2277883.785
5412539.46  2258225.544
5425763.452 2287138.979
5437312.03  2250401.825

I don't know the zone but from looking at a UTM map I think it would be 35T.  Still not getting expected lat/longs when I do the conversion though.


texts: Creating QgsVectorLayer from Pandas dataframe

texts:
I am writing a python script for QGIS.

I have a pandas dataframe that I would like to use in order to create a QgsVectorLayer.

At the moment I am saving the data frame into a .csv and then I load it:

name_csv = 'virtuali.csv'
df.to_csv(name_csv,index=False)
path = "file:///" + name_csv + "?encoding=%s&delimiter=%s&xField=%s&yField=%s&crs=%s" % ("UTF-8",",", "Longitudine", "Latitudine","epsg:4326")
mylayer = QgsVectorLayer(path, "pandas_layer", "delimitedtext")

Is there a way to create the Pandas dataframe directly without saving and reloading it?


texts: QGIS 3.4 merge raster output only 1 px wide

texts:
I have four .xyz files with terrain data. I can convert them to TIFs in QGIS 3.4, but in either format, merging them fails. The merge output always has the X dimension of the merged rasters, but a Y dimension of only one pixel.
This happens in QGIS 2.18, 3.2 and 3.4 and whatever settings I apply. Any idea how to solve this problem?

The input rasters look like this:

The output looks like this:


texts: Pivot table - Get counts of values in field by key field

texts:
I'm trying to create a pivot table which displays counts of unique values per a key field:

I've tried a combo of PivotTable and Statistics but the results are not correct. The PivotTable doesn't produce what I expected. The pivot table produces rows with empty key field values, which shouldn't exist as all input rows have a key value. I'm confused as to what inputs are what for the tool. Statistics does work as expected.

#this code doesn't produce what I was hoping - adding field with key and new fields for each summary value
pivTab = arcpy.PivotTable_management (inputTable, keyFld, summaryFld, summaryFld, pivotTab) [0]
##note the tool produces empty strings for values where they should be null
##I fix this here before stats
sumStatsFc = arcpy.Statistics_analysis (pivotTab, sumStatsTab, "valX COUNT; valY COUNT; valZ COUNT", keyFld) [0]

I don't need the python code for the answer; just which tools to use and what their inputs are.

Note that using cursors, dictionaries, etc isn't ideal as data sets are large and memory issues arise.


texts: Group Layers Do Not Display in GeoExt LegendPanel

texts:
I have a GeoExt.LegendPanel which is working just fine for layers which are not grouped.  However, the grouped layers do not show any icon in the LegendPanel.  I have tried to create the group layer various ways, using the exact name of the layer group in Geoserver and aldo by specifying a layer string in code.  I've been looking around for a while now, but I don't see what I am missing..

Here is relevant code:

//Layer Group
var geography2 = new OpenLayers.Layer.WMS(
                    "Geography2", "http://localhost:8080/geoserver/wms", {layers: 'cite:geog_multipol_vw,cite:geog_multiline_vw', transparent: true}, {isBaseLayer: false, displayInLayerSwitcher: true, visibility: true, minScale: 27734017.045258757}
               );

//LegendPanel config
var legend = new GeoExt.LegendPanel({
                region: "east",
                title: "Legend",
                width: 200,
                autoScroll: true,
                layerStore: mapPanel.layers
            });

//Explicitly point LegendPanel to GetLegendGraphic on WMS Geoserver for GeoWebCached Layers
            var LayerRec;
            for(var l = 4; l < map.layers.length; l++)
            {
                LayerRec = mapPanel.layers.getAt(l);
                LayerRec.set("legendURL","http://localhost:8080/geoserver/wms?TRANSPARENT=true&SERVICE=WMS&VERSION=1.1.1&REQUEST=GetLegendGraphic&EXCEPTIONS=application%2Fvnd.ogc.se_xml&FORMAT=image%2Fgif&LAYER=" + map.layers[l].params["LAYERS"]);
        }

Thanks for your time in trying to help me :),

elshae


texts: Export a large image from Google earth engine platform

texts:
var aoi = ee.FeatureCollection('users/vguyir/NewBegoro');
Map.addLayer(aoi)

var Begoro = aoi.filter(ee.Filter.eq('FID', 0));
var image = ee.ImageCollection('COPERNICUS/S2')
         .filterBounds(geometry)
         .filterDate('2017-01-01', '2017-12-31')
         .sort('CLOUDY_PIXEL_PERCENTAGE')
        ;
//Map.addLayer(image, {bands: ['B4', 'B3','B2'], min:0 ,max: 3000}, 'L1');

var Sort = image.toList(20)
Map.addLayer(ee.Image(Sort.get(1)), {band: ['B1', 'B2','B3'], min:0, max: 0.3},'Landsat 1')

//Map.addLayer(geometry, {}, 'Geometry');

var clip = ee.Image(Sort.get(1)).clip(table2);
Map.addLayer(clip, {band: ['B1', 'B2','B3'], min:0, max:3000}, 'clip');

var image2 = ee.ImageCollection('COPERNICUS/S2')
         .filterBounds(geometry2)
         .filterDate('2017-01-01', '2017-12-31')
         .sort('CLOUDY_PIXEL_PERCENTAGE')
        ;
//Map.addLayer(image, {bands: ['B4', 'B3','B2'], min:0 ,max: 3000}, 'L1');

var Sort2 = image2.toList(20)
Map.addLayer(ee.Image(Sort2.get(1)), {band: ['B1', 'B2','B3'], min:0, max: 0.3},'Landsat 1')

var clip = ee.Image(Sort2.get(3)).clip(table2);
Map.addLayer(clip, {band: ['B1', 'B2','B3'], min:0, max:3000}, 'clip');

var mosaic = ee.ImageCollection([ee.Image(Sort.get(1)),ee.Image(Sort2.get(1))]).mosaic();
Map.addLayer(mosaic, {bands: ['B4','B3', 'B2'], min:0 ,max:3000}, 'mosaic');

Export.image.toDrive({
  image: mosaic.clip(aoi).select('B8', 'B4', 'B3', 'B2'),
  description: 'imageToDriveExample',
  scale: 10,
  region: mosaic.clip(aoi),
  folder : 'Tank',
  maxPixels: 3784216672400,
  crs : 'EPSG:32630'
});

I keep getting this error message:

  Export too large: specified 37842166724310 pixels (max:
  3784216672400). Specify higher maxPixels value if you intend to export
  a large area


texts: Using Python code from QGIS 3.x Graphical Modeler in console?

texts:
I want to use "Join Attribute by Location (summary)" processing algorithm in Python console of QGIS 3.10. After design my algorithm process (join sum of field name "abc" form layer 2 (point layer) to layer 2 (polygon layer), I exported it as script algorithm as follow:

from qgis.core import QgsProcessing
from qgis.core import QgsProcessingAlgorithm
from qgis.core import QgsProcessingMultiStepFeedback
from qgis.core import QgsProcessingParameterVectorLayer
from qgis.core import QgsProcessingParameterFeatureSink
import processing

class Model(QgsProcessingAlgorithm):

    def initAlgorithm(self, config=None):
        self.addParameter(QgsProcessingParameterVectorLayer('layer1', 'Layer 1', types=[QgsProcessing.TypeVectorPolygon], defaultValue=None))
        self.addParameter(QgsProcessingParameterVectorLayer('layer2', 'Layer 2', types=[QgsProcessing.TypeVectorPoint], defaultValue=None))
        self.addParameter(QgsProcessingParameterFeatureSink('Output', 'Output', type=QgsProcessing.TypeVectorAnyGeometry, createByDefault=True, defaultValue=None))

    def processAlgorithm(self, parameters, context, model_feedback):
        # Use a multi-step feedback, so that individual child algorithm progress reports are adjusted for the
        # overall progress through the model
        feedback = QgsProcessingMultiStepFeedback(1, model_feedback)
        results = {}
        outputs = {}

        # Join attributes by location (summary)
        alg_params = {
            'DISCARD_NONMATCHING': False,
            'INPUT': parameters['layer1'],
            'JOIN': parameters['layer2'],
            'JOIN_FIELDS': '[abc]',
            'PREDICATE': [1],
            'SUMMARIES': [5],
            'OUTPUT': parameters['Output']
        }
        outputs['JoinAttributesByLocationSummary'] = processing.run('qgis:joinbylocationsummary', alg_params, context=context, feedback=feedback, is_child_algorithm=True)
        results['Output'] = outputs['JoinAttributesByLocationSummary']['OUTPUT']
        return results

    def name(self):
        return 'model'

    def displayName(self):
        return 'model'

    def group(self):
        return ''

    def groupId(self):
        return ''

    def createInstance(self):
        return Model()

Now I want to copy this code in Python console of QGIS and use it, but with run this code in console, I don't receive any result. How can I use this code?


texts: How to reproject vector datasets in ArcMap?

texts: I want to change one Personal Geo-database projection to another. Kindly give me instructions to re-project vector data.


texts: MapServer (MS4W) test page from shapefile shows only blank white square?

texts:
I am trying to create a test page with mapserver and shapefile. At the moment I just want to get the map shown on my test page.
I have gotten the EXTENT coordinates by opening shapefiles with QGIS - I am sure those are correct.
I have followed this tutorial as an example.

I am accessing my test page via this link:

http://localhost:1025/cgi-bin/mapserv.exe?map=C:\MS4W\ms4w\Apache\htdocs\mapfile.map

Problem I have - map is not loaded - I only see a blank white square.

When I look at contents of ms4w\Apache\htdocs\output , I see that images are actually getting generated but all of them are white (empty), for example:

lat_pmts_region13941906093724.png
lat_pmts_region13941910095144.png
lat_pmts_region13941913405500.png
...

Shapefiles are located at this location: C:\MS4W\ms4w\Apache\htdocs\shapefiles.

An example of what is logged to error file:

[Fri Apr 11 11:19:32 2014].883000 CGI Request 1 on process 7888
msWCSDispatch20(): wrong service (none)
msWCSParseRequest(): request is KVP.
msWCSDispatch(): SERVICE and REQUEST not included
msDrawMap(): rendering using outputformat named png (AGG/PNG).
msDrawMap(): WMS/WFS set-up and query, 0.000s
msDrawMap(): Layer 0 (lat_pmts_region), 0.001s
msDrawMap(): Drawing Label Cache, 0.000s
msDrawMap() total time: 0.002s
msSaveImage(C:/MS4W/ms4w/Apache/htdocs/output/lat_pmts_region13972043727888.png) total time: 0.138s
[Fri Apr 11 11:19:33 2014].24000 mapserv request processing time (loadmap not incl.): 0.141s
[Fri Apr 11 11:19:33 2014].24000 msFreeMap(): freeing map at 0246EFD8.

ms4w\Apache\htdocs\template.html:

<title>Mapserver Template HTML Document</title>

Template

It looks like this:

C:\MS4W\ms4w\Apache\htdocs\mapfile.map:

MAP
    CONFIG "MS_ERRORFILE" "C:/MS4W/ms4w/tmp/ms_error.txt"
    DEBUG 5
    NAME lat_pmts_region
    STATUS ON
    SIZE 800 600

    EXTENT 286878.53 142087.25 799084.10 481123.75
    #   lat_pmts_region

    UNITS METERS
    SHAPEPATH "shapefiles"
    IMAGECOLOR 155 155 155

    WEB
    TEMPLATE "template.html"
    IMAGEPATH "C:/MS4W/ms4w/Apache/htdocs/output/"
    IMAGEURL "/output/"
    END

    LAYER
    NAME "lat_pmts_region"
    DATA lat_pmts_region
    STATUS DEFAULT
    TYPE POLYGON
    END
END

I will gladly provide more information if it's asked.


texts: How can I, in a general way, rescale a raster in ArcGIS Model Builder

texts:
What I want to do is rescale a few rasters from whatever they are now (i.e., values between -99 and +1000 OR -256 - +256 OR 0 - 3.23) to a normalised value (i.e. values between 0 and 1).

I would prefer to do it in model builder with standard tools, but I am not averse to writing a new one in Python.

Any info appreciated.


texts:
I'm not able to run the following PyQGIS script in Windows (it worked perfectly in OS X):

import sys
import os
from qgis.core import *
import qgis.utils
from qgis.analysis import *

app = QgsApplication([],True,"")

# supply path to where is your qgis installed
QgsApplication.setPrefixPath(os.getenv("QGISPATH"), True)

# load providers
QgsApplication.initQgis()

sys.path.append(os.getenv("QGISPATH")+"\\apps\\qgis\\python\\plugins")
import processing

# Load layers with precincts and zip codes
print os.path.exists("..\\external\\precincts\\USA_precincts.shp")

layer1 = QgsVectorLayer("..\\external\\precincts\\USA_precincts.shp", "precincts", "ogr")

if layer1.isValid():
    print "Loaded layer 1"
else:
    print "Layer 1 failed to load!"

layer2 = QgsVectorLayer("..\\external\\zipcodes\\Bnd_2008_q3_region.shp", "zipcodes", "ogr")

if layer2.isValid():
    print "Loaded layer 2"
else:
    print "Layer 2 failed to load!"

# Spatial intersection between layers, save it as a new shapefile
overlayAnalyzer = QgsOverlayAnalyzer()
overlayAnalyzer.intersection(layer1, layer2, "..\\temp\\precinctsZipcodes.shp")

QgsApplication.exitQgis()

The initial lines run fine, but then the code is not able to load the shapefiles. The output from the line print os.path.exists("..\\external\\precincts\\USA_precincts.shp") is True, but then it prints Layer 1 failed to load! when checking whether the layer was loaded correctly.

Does anyone know why it does not run in Windows? I was even able to run it in the python console (only running the lines between layer1 = QgsVectorLayer ... and the second to last line).

Note: I originally had slashes instead of backslashes. I got two answers telling me that the slashes were the problem, but it wasn't: I changed them and the result is the same.


texts: Georeference a pdf without losing image quality

texts:
My recent project is to impelent a custom and georeferenced image to a leaflet map. I finisehd my custom map using inkscape now having  a .svg and a .pdf vector data. As far as i know my next step should be the georeferencing of my image via QGIS, and thats where i need some advice. I already red some threads and tutorials but i couldnt figure out how to fix my recent problem: How to georeference a pdf with QGIS? I tried the georeference plugin for QGIS but this turn my pdf to a raster. The result is a big lack of image quality when zooming.

And furthermore i need to create a tile map with this georeferenced vector data. Is this even possible?


texts: How can I join all data of two raster datasets in ArcMap 10.2?

texts:
I have two raster datasets in ArcMap10.2. One is displaying rails and one that is displaying a 100m Buffer to both sides of it. My cellsize is 100x100m.
What I want is to combine the rails AND the Buffer around it in one raster. All tools that I used only display the overlapping rasters in the new raster dataset, but I want to have ALL data in the new raster.

I tried Multiply and Plus in the raster calculator. I tried 'Mosaic To New Raster'. I tried merge... Actually I tried all possible union, combine, merge, join etc functions that I know about in ArcGIS. Either they just don´t work or I don´t get the result that I want.

Does anyone know how I can join the two rasters to have all data in the new layer?


texts: Identifying missing pixels/ pixels with no data value in .tifs

texts:
I am looking for a way to identify missing pixels in .tifs. For example, in the attached image there is a row and a column of missing pixels across the image.

The no data/ missing pixels after some investigation appear to have the value of 0 and appear across all three bands.

Can anyone recommend a way, via a piece of software/script, that could do this in bulk, say for 400-500 .tifs at a time, rather than having to do it manually/visually for each .tif?

Ideally, I'd like something that could give me a printed output in a .txt which lists all the .tifs that have missing pixels, which I could then investigate visually.

The software I have access to is; FME, ArcMap, QGIS, Python, and excel, however I may be able to get hold of something not on that list.

I have attempted a solution using Python with numpy and gdal (shown below), however I am new to python so my script is probably trash. I've played around with using numpy array and GDAL band statistics and I've got the expected results (both methods identify 0 values when expected), however I'm struggling to find a way to loop through a folder of 400/500 images and have some sort of summary statistics for all the images rather than 400/500 individual results.

import numpy as np
from osgeo import gdal

ds = gdal.Open("single.tif)

print "[ RASTER BAND COUNT ] : ", ds.RasterCount
for band in range( ds.RasterCount ) :
    band += 1
    print "[ BAND ] : ", band
    srcband = ds.GetRasterBand(band)
    if srcband is None:
         continue
    stats = srcband.GetStatistics ( True, True )
    if stats is None:
         continue

    print "[ STATS ] =  Minimum=%.3f, Maximum=%.3f" % (
                stats[0], stats[1])
myArray = np.array(ds.GetRasterBand(3).ReadAsArray())
print myArray

np.savetxt("Output.csv")

if 0 in myArray:
    print "yes"
else:
    print "no"


texts: Export composer map into set page of an existing Word document?

texts:
I'm a new user to QGIS and I'm seeking a specific method for an exporting task.
Currently, we have a 20 page Word document report with a blank page 5 for a location map to be inserted into. Once the document has been manually checked for errors, we convert it to PDF and combine it with the location map PDF (the map is exported from QGIS composer separately). We do 100s of these reports a day. Each location map is of a different area, and each report is different too.

Now, what I want to do is export the location map directly from QGIS into PAGE 5 of the existing Word document report. It will save us valuable time this way.

What's the best way of accomplishing this? I can't find any solutions online or by myself. Is there a QGIS plugin available? Or a python method someone knows of?


texts: Shapefiles - ability to toggle editing and filtering

texts:
I have a shapefile from a county HER which I have successfully loaded.   I have catagorised the point so depending on type they show a different symbol.  I can filter those points in order to show them by period etc.  However, the toggle editing option is not there and I need this to manual move and place labels.

So I save the file as ESRI shapefile. Categorise and the toggle editing is now showing for me to manual move labels.  However, the ability to filter is not there it is greyed out.

Is there any way around this?


texts: Converting GPS points with timestamp to polylines in PostgreSQL?

texts:
I have some GPS coordinates points based on timestamp in a table. My table consists of and id, time_stamp and GPS points as the geometrical data.

I want to create a line by joining all the points from this table. I assume I would have to convert the point table into a polyline table. Should I find out all the points in a particular area using St_Collect() function?

I also think that as table has time stamp related to each coordinate, I should follow the timestamp to get the next point in a line?

Putting it in other words, I have created a table with columns as id(int) , time(int) and geom(GPS coordinates as POINT). I have a set of rows added to this table. Now in order to create a polyline I assume that I should take the time in increasing order and join all the points associated with. Is my view correct?

If yes, I would like to convert these GPS points into a polyline. How should I go about it? Is there any geometric function in PostgreSQL that could support this idea?


texts: Saving outputs from arcpy.da.Walk to the same folder

texts:
I have a number of sub-directories within my "BP Tiles" directory. Each of these sub-directories represents a different year (2016, 2017, 2018, etc.) and within each of these sub-directories is a list of raster layers. I am attempting to use the Walk module to mosaic each list of Raster layers within the individual year directories.
See Code Source

import arcpy, os
from arcpy import env

# SET GEO_PROCESSING ENVIRONMENTS
arcpy.env.workspace = r"C:\Users\BP_Tiles"

walk = arcpy.da.Walk(env.workspace, topdown=True, datatype="RasterDataset")

for dirpath, dirnames, filenames in walk:
    print
    dirpath
    rasterList = []
    for file in filenames:
        raster = os.path.join(dirpath, file)
        rasterList.append(raster)
    try:
        arcpy.MosaicToNewRaster_management(rasterList, dirpath, + "Mosaic.tif", "", "", "", 1)
    except:
        pass

The code is set up to save the output "Mosaic.tif" within each sub-directory. I'd like to manipulate this code such that all of the Mosaics are saved to the same output folder. Id also like to append the sub-directory name (Ex. 2016) to each output such that I know how to differentiate the outputs.

Would anyone happen to know how to do this?


texts: How to install an old QGIS Master (1ab5cdc)?

texts: I need to revert my current version of the QGIS Desktop Master back to 1ab5cdc (about 5 months ago).  Anyone know how I can make this happen?  I'm running windows XP


texts: Making calls the PostGIS functions with Java

texts: We have Java stored procedure in Oracle that make calls to functions offered by Oracle Spatial to perform spatial analysis and return results. I am looking to the same in PostGIS, is this possible?


texts: Adding QGIS Server to an existing Windows web server

texts:
What is the best way to add the QGIS Server capabilities to an existing web server running under Windows?

Specifically, I have WAMPServer 2.2 installed and running on Windows Server 2012, and would now like to add the capability to serve WFS layers. On a different machine (running Windows 7) I have successfully set up QGIS Desktop and Server via the OSGeo4W installer, but now I want to achieve the same result without installing a second version of Apache.

The tutorials I have found are either for linux/unix-type systems, are several years old and use pre-2.0 versions of QGIS, or use the OSGeo4W installer to create a new Apache installation.

I haven't been able to find an individual installer for a Windows version of QGIS Server, and it doesn't appear to be part of the QGIS Standalone Installer.


texts: Creating new layers using ArcPy to define Definition Query for each field value?

texts: I am trying to set multiple definition queries for one layer based on one field using python.  ie.  I have a City layer with a field that holds the city name.  I want a copy of the layer in an MXD with the Definition query set for each value in the City field.  So essentially, I have one source layer, with a copy of the layer as many times as there are values in the City field.  The query would look like this: CITY = 'Atlanta'


texts: How to find coordinates for a lot of locations based on the town or village name?

texts:
I have got an Excel spreadsheet with roundabout 15,000 locations, which are only identified by the towns or villages name and the country. I am intending to present these on a world map.

Therefore, I need to know, whether there is a tool, which can find the coordinates for me based on the location's name.

Pasting every name into Google maps is no solution due to the high number of locations.


texts:
I am using Leaflet with the MiniMap plugin:
https://github.com/Norkart/Leaflet-MiniMap

Often, the zoom level difference between main map and mini map is so big that the bounding box shown in the minimap cannot be seen.

Can I set a marker on the minimap in this case? Or are there other solutions e.g. a minimum size for the bounding box on the mini map?


texts: QGIS 2.14 Composer incorrectly renders svg text position

texts:
I just noticed that SVG rendering in composer since upgrade to QGIS 2.14 is incorrectly handling some text placement.  I have an SVG compass rose, which was rendering fine before upgrading to 2.14 (Essen).

Before upgrade -

Now as displayed in Composer, 2.14 -

Note that some text is correct, the problem seems to be with the numerals arrayed around the circles, which have been transformed by rotation.

This svg was created in Inkscape, exported as a standard svg.  It is rendered correctly in any browsers and other viewing software I have tried.  It also rendered correctly in QGIS version 2.12.  I am using both Inkscape & QGIS on the same computer, running Ubuntu 14.04, so referenced font is available.

As a workaround, I have exported the svg from Inkscape to a png bitmap, which displays fine in Composer.  Another workaround is to convert the problem text to paths in inkscape, which then render as expected, but that seems a little clumsy.

Has anyone else experienced this, or have any ideas on the problem?  I really like the ability to create detailed svgs in Inkscape to greatly enhance the graphic capabilities of Print Composer, but not so handy if I can't rely on the text rendering.

Here is a copy of the problem svg, in case some curious soul wants to check for themselves - NauticalRoseSJ_ink.svg


texts: Workflow to create global mapbox tileset

texts:
I am in process of creating a hiking map rendered from OSM data using Mapbox. It works quite fine already, however some data are missing in the Mapbox Streets V7 datasets, namely:

bus stops
power lines
marked hiking/cycling routes

I am able to get all those data using Overpass Turbo for small regions and upload them as tilesets, however the largest area I am able to process this way is about 200x200 km for highly detailed areas like Europe.

Given the quite generous 50 GB limit for tilesets, I think I should be able to achieve a global (or at least continental) coverage with custom tilesets for those features. However I am not sure how to proceed. I have found an article which describes a workflow like this, however it does not show any implementation details:

Building a Custom Trail Map Using OSM and Mapbox

I think it could be possible to build some tool which would download data using Overpass API part by part. I understand the tools needs to avoid being too heavy on Overpass servers and I am fine with that, I can develop my queries using small region only and I am fine with downloading global data for a few weeks if needed. However I am bit concerned with issues about duplicate / missing data on the part boundaries. Is there any tool which would be able do handle the splitting and merging of the Overpass requests and responses (or at least merging the responses)? If not, is there any advice how to build such tool?


texts: For loop with SearchCursor stops after first iteration

texts:
I'm trying to iteratively select and export a layer using da.SearchCursor, but the loop stops after the first iteration. Here is the code:

arcpy.env.workspace = ScratchGDB #
infc= "PA"
field = 'UniqueID'
arcpy.MakeFeatureLayer_management(infc, "infcLyr")

for row in arcpy.da.SearchCursor(infc, field):
    iterateValue = row[0]
    exp = field + " = '" + iterateValue + "'"
    outfc = ScratchGDB + iterateValue
    if arcpy.Exists(outfc):
        print ("skip - output exists") # skip row if it has been exported already
    else:
        selectRow = arcpy.SelectLayerByAttribute_management("infcLyr", "NEW_SELECTION", exp)
        arcpy.CopyFeatures_management("infcLyr", outfc)

However, when I use a similarly structured code, but instead of exporting, I just printed the output path, the loop worked fine:

arcpy.env.workspace = ScratchGDB #
infc= "PA"
field = 'UniqueID'
arcpy.MakeFeatureLayer_management(infc, "infcLyr")

for row in arcpy.da.SearchCursor(infc, field):
    iterateValue = row[0]
    exp = field + " = '" + iterateValue + "'"
    outfc = ScratchGDB + iterateValue
    if arcpy.Exists(outfc):
        print ("skip - output exists") # skip row if it has been exported already
    else:
        print(outfc)

Update:
I tried to set the arcpy.env.overwriteOutput to FALSE following the suggestion by @Emil Brundage. The loop still stops after both SelectLayerByAttribute_management and CopyFeatures_Management  run for the first time.
Note that the first part of the if statement loops with no problem; i.e., if there are multiple existing outputs, the loop will skip all of them. It only stops after the select and copy tools run through their first input row.


texts: Join Attribute Table with multiple owners with only one Parcel

texts:
I have a problem in QGIS. I have an attribute table (no geographical information), which contains owners of parcels and an ID for every row.

Then I have a shapefile, which contains polygons and an attribute table which contains an ID (Fits with the owner information).

If I want to join these two parts, QGIS joins only one owner per parcel, even if there are 5.

How can I manage it, that finally there is a row in the attribute table of the shapefile, containing all owners and not only one?


texts: Merge raster tiles with different resolution and keep original resolution using Python

texts:
Problem

I have two raster images with different size and resolution, which are overlapping and which I'd like to merge. The smaller image has a higher resolution which I want to merge into the larger image with the lower resolution. However, I don't want to resample the low resolution raster tile, but keep the origin resolution in the area where no overlapping exists.

Information

Some information about the two tiles:

first image (small, high resolution):

Driver: GTiff/GeoTIFF
Files: /home/Downloads/srtm_38_03/srtm_38_03.tif
Size is 6000, 6000

Coordinate System is:
GEOGCS["WGS 84",
DATUM["WGS_1984",
SPHEROID["WGS 84",6378137,298.257223563,
AUTHORITY["EPSG","7030"]],
AUTHORITY["EPSG","6326"]],
PRIMEM["Greenwich",0],
UNIT["degree",0.0174532925199433],
AUTHORITY["EPSG","4326"]]

Origin = (5.000000000000000,50.000000000000000)
Pixel Size = (0.000833333333333,-0.000833333333333)
Metadata:
AREA_OR_POINT=Area
Image Structure Metadata:
INTERLEAVE=BAND
Corner Coordinates:
Upper Left ( 5.0000000, 50.0000000) ( 5d 0' 0.00"E, 50d 0' 0.00"N)
Lower Left ( 5.0000000, 45.0000000) ( 5d 0' 0.00"E, 45d 0' 0.00"N)
Upper Right ( 10.0000000, 50.0000000) ( 10d 0' 0.00"E, 50d 0' 0.00"N)
Lower Right ( 10.0000000, 45.0000000) ( 10d 0' 0.00"E, 45d 0' 0.00"N)
Center ( 7.5000000, 47.5000000) ( 7d30' 0.00"E, 47d30' 0.00"N)
Band 1 Block=6000x1 Type=Int16, ColorInterp=Gray
NoData Value=-32768

second image (large, low resolution):

Driver: GTiff/GeoTIFF
Files: /home/Downloads/tiles/30N000E_20101117_gmted_mea075.tif
Size is 14400, 9600

Coordinate System is:
GEOGCS["WGS 84",
DATUM["WGS_1984",
SPHEROID["WGS 84",6378137,298.257223563,
AUTHORITY["EPSG","7030"]],
AUTHORITY["EPSG","6326"]],
PRIMEM["Greenwich",0],
UNIT["degree",0.0174532925199433],
AUTHORITY["EPSG","4326"]]

Origin = (-0.000138888888889,49.999861111111109)
Pixel Size = (0.002083333333333,-0.002083333333333)
Metadata:
AREA_OR_POINT=Area
Image Structure Metadata:
INTERLEAVE=BAND
Corner Coordinates:
Upper Left ( -0.0001389, 49.9998611) ( 0d 0' 0.50"W, 49d59'59.50"N)
Lower Left ( -0.0001389, 29.9998611) ( 0d 0' 0.50"W, 29d59'59.50"N)
Upper Right ( 29.9998611, 49.9998611) ( 29d59'59.50"E, 49d59'59.50"N)
Lower Right ( 29.9998611, 29.9998611) ( 29d59'59.50"E, 29d59'59.50"N)
Center ( 14.9998611, 39.9998611) ( 14d59'59.50"E, 39d59'59.50"N)
Band 1 Block=14400x1 Type=Int16, ColorInterp=Gray
Min=-81.000 Max=4712.000
Minimum=-81.000, Maximum=4712.000, Mean=284.827, StdDev=414.286
NoData Value=-32768
Metadata:
STATISTICS_EXCLUDEDVALUES=-32768
STATISTICS_MAXIMUM=4712
STATISTICS_MEAN=284.8270891059
STATISTICS_MINIMUM=-81
STATISTICS_STDDEV=414.28584583842

The image shows the two raster tiles. The orange box is just for orientation where the smaller tile is localized.

This image cutout shows the two different resolutions at the tiles border which I want to keep in the merged tile.

What did I try so far?

  gdal_merge: This utility will automatically mosaic a set of images. All the images must be in the same coordinate system and have a matching number of bands, but they may be overlapping, and at different resolutions. In areas of overlap, the last image will be copied over earlier ones.

I used gdal_merge.py in the following way:

/usr/bin/gdal_merge.py -o /home/Workspace/output.tif 30N000E_20101117_gmted_mea075.tif srtm_38_03.tif

which results in an image identical to 30N000E_20101117_gmted_mea075.tif with no change at all.

Using it with the images switched to receive an output with the high resolution like described in:

  -ps pixelsize_x pixelsize_y
  Pixel size to be used for the output file. If not specified the resolution of the first input file will be used.

using the command:

/usr/bin/gdal_merge.py -o /home/Workspace/output.tif srtm_38_03.tif 30N000E_20101117_gmted_mea075.tif

results in the following error:

0...10...20...30...40...50ERROR 2: gdalnodatamaskband.cpp, 244: cannot allocate 3456000000 bytes
[1]    19672 segmentation fault (core dumped)  /usr/bin/gdal_merge.py -o /home/Workspace/output.tif

Any ideas? Is it even possible to realize what I want to do?


texts:
I'm changing the stroke color of my features when they are selected like follows:

  var select = new Select({
    condition: click
  });
  select.on('select', function (e) {
    e.selected.forEach(function (feature) {
      feature.getStyle().getStroke().setColor([r, g, b, 1]);
    });
  });

It is working but when the map is translated, the deselected feature doesn't change its stroke.

I tried to use also the e.deselected property to set the default stroke color but it didn't work


texts: How to add Time info to GeoTiff served from Geoserver?

texts:
I have a geotiff file (somedata.tif) that I want to serve from geoserver wms associated with the time 10/22/2010 12:00:00.

I'm currently using a geoserver raster data source of type GeoTiff pointed to somedata.tif with a layer pointed at that data source.  Using this configuration, I can successfully serve the image, with no time information.

I'm looking for functionality that would be similar to mapserver where I set wms_timedefault=10/22/2010 12:00:00.

My desired end result is an http header like this:

Warning: 99 Nearest value used: TIME=2010-10-22T12:00Z

Can I add time information to a Geotiff data source in geoserver v2.0.2?  Does it belong in somedata.tif? Does it belong in my source or layer configurations?

My GIS vocabulary isn't 100% so if I need to clarify anything, please let me know.


texts: Can I use google app engine to work as a WMS or some sort of geoserver that I can use to connect openlayers?

texts: Has anybody ever heard or tried doing anything like this?


texts: Trying to delete identical values within an update cursor

texts:
I'm trying to use the DeleteIdentical_management tool within an update cursor. I'm pretty sure I'm going about this the wrong way, I'd like to see if someone could help.
I'm relatively new to python, so bear with me. Here is my bit of code:

#Delete Identical Features based on Field Name and Rig Number field, not deleting field names with null values.
FC = r'C:\Scripts\Drilling_Rig_Deduplication\Drilling_Rigs.gdb\MERGED_RIGS'
cursor = arcpy.UpdateCursor(FC)
for row in cursor:
   fieldValue = row.FIELD_NAME
   if row.FIELD_NAME is not None:
     arcpy.DeleteIdentical_management(FC, ["FIELD_NAME","RGNMBRCLC"])
   else:
    print "nothing left to delete"
del row, cursor

print "dedup process done"
print arcpy.GetMessages()

What I'd like to happen is for arcpy to ignore the rows where there is a NULL value in the FIELD_NAME field. Then have the delete_identical tool look through the rest of the rows and delete duplicate records. - Any ideas?


texts: How to get sum / total length of lines per polygon in QGIS?

texts:
I have a line vector layer which has multiple segments within individual polygons in a separate polygon layer. I have spatially joined the attributes so that each line segment now has the same ID as its corresponding polygon.

I need to calculate the lengths of the lines within each polygon according to the polygon ID. Is there an expression I can use to do this? I have already calculated the lengths of the individual segments but I need to get the sum of these according to each unique polygon ID.

I am a novice at complex SQL expressions within QGIS, ideally I would like a new column displaying the info.


texts: Join attributes by location very slow on QGIS 3.0

texts: "Join attributes by location" is very slow on QGIS 3.0. I have tried it with a CSV file and also with a shapefile (point layer) and a grid (made of polygons). On QGIS3.0 it takes about 3 hours to complete, but on QGIS2.18.16 only 7 seconds. The calculated results look exactly the same. Am I missing something? Here are the screenshots of the setup (QGIS3.0 in picture 1 and 2; QGIS2.18.16 in picture 3):


texts: Generating random locations nearby my location in PHP and MySQL

texts:
I'm using a method to generate a random locations around my location with given the initial position(lat,long) and radius this method is used  From the  reference post here,  but unfortunately the random locations not generated and display on Google Maps. In my code I had passed the value of lat and long of initial position and given radius around 100 meter around. The code as shown below. It seems I've missed some point in the code.

The code in PHP

            var lat=40.744656;
            var long=-74.005966;
            var southWest = new google.maps.LatLng(lat,long );
                    var r = 100/111300 // = 100 meters
                    , y0 = lat
                    , x0 = long
                    , u = Math.random()
                    , v = Math.random()
                    , w = r * Math.sqrt(u)
                    , t = 2 * Math.PI * v
                    , x = w * Math.cos(t)
                    , y1 = w * Math.sin(t)
                    , x1 = x / Math.cos(y0)
                    newY = y0 + y1
                    newX = x0 + x1
            // set multiple marker
            for (var i = 0; i < 20; i++) {
                // init markers
                var marker = new google.maps.Marker({
                    position: new google.maps.LatLng(newY, newX ),
                    map: map,
                    title: 'Click Me ' + i
                });


texts: create lines from points by position of the points and orientation

texts:
I want to create straight lines based on the points. I used function "Points to Line" which accomplish the task only when the points are sequentially ordered. If they are not ordered, then it creates a zig zag lines (I have included the figures for illustration - Line 2). I have hundreds of lines, and manually editing is not an ideal solution.

How can one create a straight line without manually editing the lines?

I believe one has to take into account of the orientation and the position of the points to address this problem.

[![enter image description here][3]][3]

I also tried running the EMST (Euclidean minimum spanning tree) function in the spanning tree tool. I could connect several lines in the same direction by using "Unsplit Line", but don't know what to do with the lines that connect parallel lines. These are unnecessary lines.