texts: Wrong coordinates after Usage of NumPyArrayToRaster(...)

texts:
My Project I am currently working at should enable cellvalue manipulation on raster files using ArcPy. My Script works fine except I am not getting the right Y-Values (northing values) for the modified raster image.
I will show you my code:

# Import arcpy module
import arcpy

# Check out any necessary licenses
arcpy.CheckOutExtension("3D")

# Local variables:
inputRaster = r'G:\Data\filegdb.gdb\ClipRaster2'
firstCoord = "33321792 6012729"
secondCoord = "33321809 6012763"
level = "5"
outputRaster = r'G:\Data\filegdb.gdb\RasterModTestData2'

# Get Raster Description:
desc = arcpy.Describe(inputRaster)
xmin = desc.extent.XMin
xmax = desc.extent.XMax
ymin = desc.extent.YMin
ymax = desc.extent.YMax
cellsize = desc.meanCellWidth
nodata = desc.noDataValue
sr = desc.SpatialReference
lowerleft = arcpy.Point(xmin,ymin)

xwert1 = firstCoord.split(" ")
UTMRechts1 = float(xwert1[0])
UTMHoch1 = float(xwert1[1])
xwert2 = secondCoord.split(" ")
UTMRechts2 = float(xwert2[0])
UTMHoch2 = float(xwert2[1])

# Converts the Raster to a NumPyArray
raster = arcpy.Raster(inputRaster)
array = arcpy.RasterToNumPyArray(raster,nodata_to_value=raster.minimum)

# Loop through the NumPy-Array and check if the point is inside the bounding rectangle
i = ymin
countX = 0
while i < ymax:
    j = xmin
    countY = 0
    while j < xmax:

        if (j >= UTMRechts1 and j <= UTMRechts2 and i >= UTMHoch1 and i <= UTMHoch2):
            array[countX,countY] += int(level)

        countY += 1
        j += cellsize
    countX += 1
    i += cellsize

newraster = arcpy.NumPyArrayToRaster(array, lowerleft, cellsize, cellsize, nodata)
arcpy.DefineProjection_management(newraster, sr)

if arcpy.Exists(outputRaster):
    arcpy.Delete_management(outputRaster)

newraster.save(outputRaster)

As You can see, the area bounding by a rectangle (lower left: 33321792 6012729, upper right: 33321809 6012763  //UTM 33N (ETRS89) )  should be modified with raise of its elevation by 5 metres. The following screenshot shows a extract from the result raster:

Unfortunately, the lower left corner of the modified area is at (33321792 6012708 metres UTM33N/ETRS89) and the upper right corner is at (33321809 6012743 metres UTM33N/ETRS89). So the easting of the coordinate fits perfectly to the requested extent, but the northing is displaced by about 21 metres to the south.
The input raster has one band, a cellsize of 0.5 m and 2254 columns and 2532 rows.

Additional Info: If I print the coordinates which are inside the bounding box, the values I want to have show up. So I guess, the mistake is in the retransformation from NumpyArray to the raster file. Has anyone experience with this problem?

Can anybody help me what to do to get the correct extent of the modified area or show me the mistake in my code?


texts: Field Calculator error: "calculated value is invalid"

texts:
I'm getting the following error when trying to calculate a field:

I'm just trying to calculate all the selected records to equal a value I specify, a simple two-letter zoning code.  So the field calculator expression is simply field = value (I'm not trying to make it equal another field or use a mathematical espression).  The field is string type, length 10, which should be long enough for two letters.  The field does not have a domain assigned to it.  I get the same error whether I'm in an editing session or not.  This should be easy - what's going on?


texts: QGIS Composer - Exported Map layer diagonal fill becomes solid (dense)

texts:
When I export my compiled map from QGIS to JPEG or PDF the fill pattern becomes solid rather than diagonal lines.  I am exporting at 600dpi when to JPEG.  But it is clear that the diagonal lines are still in the image but just highly condensed creating the appearance of a solid polygon feature. Which suggests to me that I am missing some export scale setting somewhere?

Also the point class features appears less clear in the exported image below.

Below is a screenshot of the map in composer:

When I Export the map to either JPEG or PDF the following is the result - note the change from a diagonal fill to a solid and reduced point feature clarity:

This is not a major problem but the idea is exporting the map that I spent a few hours compiling rather than a slightly different variant of the map.

I hope that it is a simple setting somewhere that I am unaware of that is condensing the fill pattern when exporting -- I am no GIS expert.


texts: Using for loops in R with Geospatial Modelling Environment

texts:
I am using the isectpolypoly command in GME and would like to iterate over several polygons ("proj_0_ect0207", "proj_1_ect0207"). My code is as follows:

for (i in 0:30)
{

isectpolypoly(in="F:\Users\p3\ResidentialShapefile\Britain_BNG.shp",poly=
("F:\Users\p3\new.gdb!"Pro_",[i],"_ject0207"), field="GRIDCODE", prefix=[i],
thematic=FALSE, proportion=FALSE, awm=TRUE, min=FALSE, max=FALSE, aws=FALSE);

};

I keep getting the error:

  Error: For loop commands must be contained within a pair of braces: { }.

I'm following guidance from http://www.spatialecology.com/gme/for.htm


texts:
I am trying to save edits to a raster file in QGIS, but;

When you edit a raster file and select "Save As"; no output file is generated.
My version is in Spanish; however "Save as raster layer" is displayed. Just as "Output"

How do I properly save/export a raster in QGIS?


texts: Trying to merge convert mult linestrings with breaks

texts:
I am trying to clean up a road network file. The roads are encoded as multlinestrings I want them as linestrings.

Normally I would use st_mergeline to convert from a multilinestring to linestring. But my issues is that some of the multilinestring include 1 or more holes/breaks in the linestring.

Where the two lines represent two  multilinestring I am just wondering if there is a sly way of inserting a bridge between the two strings.  This issue is come with  about 30% of my data.  I could resolve this issue by editing data with qgis.  But there must be a better way of sorting this problem out that hand editing all that data. Using functions like st_union result in a single geometry  object with it includes the break. I want to joint the two segments together.


texts: Extracting flat areas and saving them as lakes using QGIS?

texts:
I would like to extract lakes and sea from a DEM (SRTM) and save them as a separate layer.

Is this possible?

I have tried to extract flat areas with the Raster Calculator, but that did not work.


texts: Running MFmask in MATLAB: Error using imfill

texts:
I am trying to create a cloud and cloud shadow masking raster using MFmask (https://www.sciencedirect.com/science/article/pii/S0034425717303073) for a Landsat 8 scene. MFmask 1.1 is modified version of Fmask 3.3 for better detection of cloud shadow in mountainous region using DEM. Previously I run Fmask 3.3 standalone and that worked, but the cloud/cloud shadow detection wasn't good enough and I moved to MFmask. It has no standalone, thus it has been run through MATLAB, but I got following error during the running in MATLAB.

>> addpath('D:\INSTALL\Cloud detection\MFmask\MFmask-master\MFmask-master\MFmask_1_1_beta')
>> autoMFmask

MFmask 1.1 version start.
Cloud, cloud shadow, and snow detection for Landsat 8 imagery.
Cloud/cloud shadow/snow dilated by 3/3/0 pixels (default).
Cloud probability threshold of 22.50% (default).
Loading input data.
From DNs to TOA ref & BT
Detecting clouds. This might take some minutes so be patient.
Error using imfill
Expected input number 1, I1 or BW1, to be non-NaN.

Error in imfill>parse_inputs (line 165)
validateattributes(IM, {'numeric' 'logical'}, {'nonsparse' 'real','nonnan'}, ...

Error in imfill (line 123)
[I,locations,conn,do_fillholes] = parse_inputs(varargin{:});

Error in plshadow2 (line 27)
    data_nir_filled=imfill(data_nir);

Error in plcloud_ratelapse2 (line 228)
        Shadow=plshadow2( data_nir,data_swir,idlnd,mask,l_pt,...

Error in MFmask (line 5)
        =plcloud_ratelapse2(fullfile_path,im,demtype,cldprob,num_Lst);

Error in autoMFmask (line 108)
    clr_pct = MFmask(fullfile_path,norMTL.name,demtype,cldpix,sdpix,snpix,cldprob,num_Lst,num_near,gap_off); % newest version
    1.1

How can I handle these issues? Such errors could be found during running of Fmask as well.


texts: how to solve UnicodeDecodeError in arcpy?

texts:
I need to update datasource of Picture element of layout element of a mapdocument. but it is showing Unicode decode error.

    mxd= arcpy.mapping.MapDocument("CURRENT")
    for elm in arcpy.mapping.ListLayoutElements(mxd,"PICTURE_ELEMENT"):
         oldname = elm.name
    ...     if elm.name==oldname:
    ...         elm.name=oldname.replace(r"ä","ae")
    ...         print elm.name

after running the code I get following error message:

  Runtime error
  Traceback (most recent call last):
    File "", line 4, in
  UnicodeDecodeError: 'ascii' codec can't decode byte 0xc3 in position 0: ordinal not in range(128)


texts:
This seems like a simple question but I can't solve this issue.

I want to convert the following Polygons (see Image) into Raster using QGIS and used the Tools 'rasterize' from QGIS and SAGA as well as GRASS v.to.raster. They all give me the same result:

I want the raster to be only under the individual Polygons.
The problem is that I used the Polygon Layer as Output extent of the tool. But what settings can I use so that the raster is only under the polygon?

I also tried to clip this result with the tool "Clip raster by mask layer" but nothing happened. I got a new raster layer with the exact same extent.


texts: Is QGIS suitable to create maps with image popups?

texts:
I am a new user of QGIS and was wondering if I came across the right software to achieve my goal.

Long story short, I am trying to create a map with markups of various buildings. Once I click on a mark, I would like to see details of that property and possibly a photo.
Further more I would like to integrate this map on my website.

Can you please advise me if QGIS can help me achieve this? Or would you have a better recommendation? I can see that the software capabilities are far more advanced than the purpose of my goal, but thought of asking for your advice.

Many thanks.

George


texts: QGIS - creating cross section to meet contour line

texts:
How can I automate the process of creating a shapefile containing cross-sections of a river centreline that will meet to a set distance given by a contour shapefile?

Is there a plugin or a set of functions that will allow me to set the cross-sections to start perpendicularly from a centreline, continue perpendicularly for a set number of meters, then bend till they snap to another shapefile (contours). Something like the pic below:

I use QGIS 3.6 and the final purpose is to automate as much as possible the creation of cross sections, refine what needs to be refine and then export them to HEC-RAS.


texts: QField - create a form for multiple, 1:n photos

texts:
I have read the QField support on how to add a series of pictures to a feature. I've also read other threads and GIS Stack Exchange questions and answers on this: e.g., Setting QField for 1:n photos. However, I am still confused and unable to get this working properly and was hoping that someone could provide a detailed explanation that may help others in this process.
Step 1: I create the relation between the Referenced 'parent layer' (= Wetland Survey) and the referencing 'child layer' (= Wetland_pics):

I am not entirely clear on what 'uuid' actually is here. I know that it is a Universal Unique Identifier, but is this necessary? Why couldn't I use my site.code? I notice that using site.code()as a default value gives an error, where uuid() gives an identifier - so the system must recognize uuid for this function. This is not my primary question, but understanding this seems important in relation to the primary keys (not entirely sure of what they are).
Step 2. Building the form layout
a. Drag and Drop, set the uuid (parent relation):

b. Drag and Drop, set the Wetland Pictures (child relation):

Step 3. Setting drag and drop form in the child table ('Wetland_pics')  - set path to attachment widget type:

Step 4. Testing the form:

The form looks fine, but I am unable to select or type anything in the "Wetland Pictures" fields in QGIS. I'm trialing the form in QField beta windows version (second image) and you can see that I get the message "Cannot add child feature: parent primary keys are not available". This is what leads me to suspect that this has something to do with uuid.


texts:
I am trying to get the proper labels for pie charts on my map. I want to use the column names of the table (.csv file) that is used for the pie chart as labels for the pie chart on the map. The column names show up in my layers pane as categories, but I am having trouble making them the labels of the pie chart.

How can this be done?

I have added the .csv file as a layer through "add delimited text layer" and have not saved it as a shapefile.

Map with pie charts and the layer pane showing the labels I would like to use:

Format of my attribute table:


texts:
Following is my raster file and mask shapefile

files

I use QGIS 1.8.0 to clip the raster with shapefile by GDAL tools, but after I run warp in GDAL tools it's just give me no data in the output GeoTiff.


texts: Calculate field based on a set value if optional parameter is not specified

texts:
In my model, I have a Calculate Field which will populate the field with a certain integer (in this case, 0). I've added an optional integer model parameter Start ID, where the user can specify the integer if they want.

If Start ID is given a value by the user, this value must be used as the expression for Calculate Field, else Expression must be used. Since the tool can only take one expression as an input, would I have to create a copy of the Calculate Field with Start ID as the expression, and use Merge Branch with an if statement to test which branch must run i.e. the Start ID branch only runs if it has a value, else run the Expression branch?


texts:
I am using OS Translator to convert .gz OSMasterMap files to .SHP files for use in QGIS. I am now trying to add QGIS .QML files to add the styling within QGIS, I have followed the instructions in the Getting Started Guide (https://github.com/OrdnanceSurvey/OSMM-Topography-Layer-stylesheets/blob/master/Getting%20Started%20Guide%20-%20Styling%20OSMM%20Topography%20Layer.pdf) however when I add the layers using the .QML file the layer becomes greyed out and the features disappear (below).

I am using the Schema 9 stylesheets from https://github.com/OrdnanceSurvey/OSMM-Topography-Layer-stylesheets.

Does anyone know the reason for this issue or what I may have done wrong?

EDIT


texts: Anybody can help me to find the three color palettes of the following map?


texts: Relationship between attribute of vector data and raster variable in ArcGIS for Desktop?

texts:
I have a vector data of individual elements (buildings) with attributes (value, etc), and several raster grids with continuous variables. I would like to understand whether there is a relationship between the attributes of a building and the other variables - my understanding is that this would require  spatial multivariate regression.

This is a new topic for me.

Would anyone be able to advise where to begin such an analysis in ArcGIS for Desktop?


texts: Serve multiple map styles with mapnik and mod_tile from one server & database

texts:
I struggle to set up two different map styles for the mapnik-osm rendering tool chain. According to some help I already found on gis.stackexchange, I set up the renderd.conf for two map styles. But I can't get it to work. When I enter the URL for the second map style (http://myserverurl/print/0/0/0.png) the browser says "The requested URL was not found on this server." With the default map style it is no problem (http://myserverurl/osm/0/0/0.png).

I couldn't find any other setting in the mod_tile or renderd config files to enable the "print" map style to work. The renderd.conf looks like this:

[renderd]
stats_file=/var/run/renderd/renderd.stats
socketname=/var/run/renderd/renderd.sock
num_threads=10
tile_dir=/var/lib/mod_tile

[mapnik]
plugins_dir=/usr/lib/mapnik/3.0/input
font_dir=/usr/share/fonts/truetype/ttf-dejavu
font_dir_recurse=false

[default]
URI=/osm/
TILEDIR=/var/lib/mod_tile
TILESIZE=256
XML=/home/osm/openstreetmap-carto-4.24.1/style.xml
DESCRIPTION=Default map style
;ATTRIBUTION=&copy;OpenStreetMap and contributors, CC-BY-SA
;HOST=http://myserverurl
;SERVER_ALIAS=http://a.tile.openstreetmap.org
;SERVER_ALIAS=http://b.tile.openstreetmap.org
;HTCPHOST=proxy.openstreetmap.org

[print]
URI=/print/
TILEDIR=/var/lib/mod_tile
TILESIZE=512
XML=/home/osm/openstreetmap-carto-4.24.1/print/style.xml
DESCRIPTION=print map style
;ATTRIBUTION=&copy;OpenStreetMap and contributors, CC-BY-SA
;HOST=http://myserverurl
;SERVER_ALIAS=http://a.tile.openstreetmap.org
;SERVER_ALIAS=http://b.tile.openstreetmap.org
;HTCPHOST=proxy.openstreetmap.org


texts: NE Attribute editor not honouring custom selection

texts:
I have an application based on Network Engineer & ArcObjects.

A part of the functionality requires users to draw a rectangle an select features. It should "add" to the current selection.

I am using a custom selection method (incorporating QueryFilter & IFeatureSelection) and selecting features in the layers. Once I am done with the selection, I open the NE Attribute Editor to see the properties of selected records, but I don't see any of the newly selected features!

I tried calling the SelectionChanged() and MakePermanent() methods and even refreshed the selection but nothing helps...

Any ideas?

PS: I have yet not tried the built-in commands. Will need some time to check if it works (I am sure it will as it is "built-in", but I don't see a reason why the feature selection does not work).


texts:
Modis MYD03 and MOD03 files contain Land/SeaMask (uint8 with valid_range=0,7) and WaterPresent (uint8 with valid_range=0,8). Though the attribute fields define the valid range, they do not define the semantics of the actual values.

What do the values actually mean?
Where is a good place to search for these sort of information?


texts: Automate saving features by results of join attributes by location

texts:
In QGIS v3.8.3 I have one shapefile of all US Census Tracts and another of ~1,000 US Counties. There are multiple tracts in a county. So I did a many-to-one join in Join attributes by location to add a 'county' column to the Tracts file.

I'd like to export shapefiles for each county-based group of tracts. For example, one shapefile for the ten tracts in County 1, another shapefile for the twelve tracts in County 2...all the way to County 1,000. I'd prefer not to click Export > Save selected features as a thousand times.

Here's what I've tried and looked at so far, but I didn't think that any of them were specifically focused on batch exporting based on the value of a specific column (but I may be wrong, I'm not super Python-literate):
Batch save layers plugin
Exporting several files at same time in QGIS?

 QGIS batch conversion
 Automate SHP to JPG creation


texts: Save a copy and replace the last six characters (date format) from a list of user selected MXD documents

texts:
The following script has the following objectives:
User selected MXD files.  Here I use a multivalue parameter for FILE with a filter for MXD in the tool box settings.
Changing selected text based on Text Element Name (dynamic text).
Saving the mxd as copy : removes the last six characters (string reference to the date), and replaces these with the current date string (or a user specified input).

As it stands it saves a copy by simply adding a suffix of the date and does not remove the original ".mxd" so I get a file name that looks like "C:\Home\GIS_Projects\HabitatEnvironments_150122.mxd_150219.mxd" (note the 2 ".mxd"s)
If the original mxd was titled "C:\Home\GIS_Projects\HabitatEnvironments_140629.mxd" the new one (same folder-location) should be the same but with the current date "C:\Home\GIS_Projects\HabitatEnvironments_150219.mxd".

import arcpy
import os
import time
import datetime
from time import strftime

mxdString = arcpy.GetParameterAsText(0) #toolbox parameter workspace
#prepare for multivalue toolbox parameter
mxdList = mxdString.split(";")
arcpy.AddMessage(mxdList)
#set new text variables
#These are set parameters in text properties "Element name"
NewMapTitle = arcpy.GetParameterAsText(1)
NewProjectTitle = arcpy.GetParameterAsText(2)
CurrentDate = strftime("%y%m%d")
#THis is where it should loop through the ListInMXD and replace the text.
for mxdItem in mxdList:
    mxd = arcpy.mapping.MapDocument(mxdItem.strip("'"))
    for titleNew in arcpy.mapping.ListLayoutElements(mxd, "TEXT_ELEMENT", "MapTitle"):
        titleNew.text = NewMapTitle
    for projectNew in arcpy.mapping.ListLayoutElements(mxd, "TEXT_ELEMENT", "ProjectName"):
        projectNew.text = NewProjectTitle
    #Save a copy of each mxd with current date at the end
    #This is where it bugs.  I need to strip the last six characters and replace them with "CurrentDate".
    mxd.saveACopy(mxdItem + CurrentDate +".mxd")

del mxd


texts: I'm using QGIS, but not for a long time. I'm editing and saving shapes and geotiffs. And I'm happy how this works. I'm wondering if there's a way to import and export metadata ISO 19115, while saving and loading files. In a way that it is stored, together with the shapes or tiffs? I hate to loose the connection between data and metadata. Can this be done?


texts: How to get a set of all buildings of an OpenStreetMap section in form of single points / Multiple polygons as one point

texts:
The objective is to get a set of all buildings of an OpenStreetMap section in form of single points in a shapefile.

I downloaded the section from OSM.
Imported the section into QGIS (and later imported it into a PostGIS db using osm2pgsql)
To transform the polygons (representing buildings) into points I used the Polygon Centroids function

Unfortunately some of the buildings are represented by multiple polygons subsequently I did get multiple points for one building which I is not what I wanted.

Searching the web I found out that these polygons should be multipolygons but neither QGIS nor osm2pgsql does import these relations so that I could merge the polygons of a building

Any ideas how I could get single point list of the buildings?


texts: Understanding land cover RGB values in greyscale GeoTIFF

texts:
I apologise if I'm not using the right terms but I'm not a GIS expert.
This page has some GeoTIFFs (at least, I'm pretty sure they're GeoTIFFs) which I'd like to analyse and I have downloaded the first one:
http://data.ess.tsinghua.edu.cn/fromglc10_2017v01.html
When I open this file in an image viewer (e.g. PhotoShop) it opens as 8 bit greyscale. When I examine a specific pixel in, say, water (since that is visually obvious in the file) PhotoShop tells me the RGB values are (80,80,80). All other land cover values seem to follow the same pattern e.g. (29,29,29), (113,113,113) etc.
There is nothing on the previous page (http://data.ess.tsinghua.edu.cn/) that seems to give any indication of a "legend", or what type of land cover each RGB value in the GeoTIFF file maps to. I am assuming this is something commonly known or understood amongst GIS experts, but no amount of searching on this SE or google searching turns up any result. I even tried to look through the GeoTIFF standard but I couldn't see anything obvious there.
Is there somewhere I can get a breakdown of the RGB values in that file to common land cover types such as water, urban, forest etc.?


texts: Is it possible to have dynamic layer names in QGIS project?

texts:
I have a QGIS project with layers, based on PostGIS SQL queries that use current date. Since the data in the database changes, these queries return different data every new day.

Is it possible to dynamically change the name of the layer in layer tree panel, so it represents the current date? (i.e. layer name like Traffic on 24.01.2015, that changes every new day). The data for the caption should be taken from attribute of any feature in a layer - they all have the same value for field traffic_date.


texts: How to create a layer of directional/horizontal wells?

texts:
I have never needed to create a horizontal/directional shapefile (oil / gas well).  What do you need besides the lat/long, basic well info?  Direction/Length for the line?  My question is - What goes into creating this shapefile? This may seem like an ignorant question, but this is the first time I've had to deal with creating this shapefile.  I've always grabbed it off the State or County sites.

I have ArcGIS for Desktop Basic (formerly known as ArcView) 10.1, so any help / guidance is greatly appreciated.


texts: Cannot log on to ArcIMS Web Manager, first time logon

texts:
My current problem with arcims web manager dot net in arcims 9.3: after installing arcims and configuring it successfully and then installing web manager dot net, I cannot log on to web manager.  Does anyone have a solution to this problem?

Screenshot

More information:

OS: Windows xp sp3.
Arcims 9.3, with java jre5 u13, ServletExec5.0ISAPI, IIS6
arcims web adf dot net 9.3, integrated VB.Net 2008.

Thank you.


texts: arcpy: check for same record in cell above and below

texts:
After using the find duplicates in arcmap 10, I want to check if the feat sequence ("find duplicates" marks all records with the same attributes with the same number [in "feat_seq"] and creates a ROWID_ that has unique numbers for each record) is the same in the cell above and below and record and if it is the same then mark the duplicate/triplicate etc as "1", if not mark with 0.

  The logic we used in excel is
  =IF(OR(C1=C2,C3=C2),1,0)

  If c2 = c1 or c2 =c3, give the cell a value of 1, if not, give the
  cell a value of 0. Used to identify duplicate values that are in
  order.  If the cell above or below is the same as the cell in
  question, give the cell a value of 1.  If not, give the cell a value
  of 0.

The code I started to write was

import arcpy
rows = arcpy.UpdateCursor("TestDuplicateID")

for row in rows:
    rowCur=row
    rowPost= row.next()
    rowPre= # how can we get this?

    if rowCur.getValue("FEAT_SEQ")==row.getValue("FEAT_SEQ"):
        rowCur.setValue("Dups", "1")
    elif rowCur.getValue("FEAT_SEQ")==rowPre.getValue("FEAT_SEQ"):
        rowCur.setValue("Dups", "1")
    else:
        rowCur.setValue("Dups","0")

I can't seem to find a way to retain the PRE record and get the following error "Runtime error : unsupported operand type(s) for -: 'Row' and 'int'".

Thanks for your help


texts: Is Survey123 the only ArcGIS field app offering public forms?

texts:
I believe from my research that, of the various ArcGIS field apps, public offline forms are only available from the Survey123 app.
Is that correct?


texts:
At the moment the code I use creates a buffer extending 1000km from each geometry creating square buffer zones over polygons, lines and points. How can I amend the code so that it creates a 'sausage' buffer around the polyline instead of creating a square extent? Should I create another variable for biggerExtent?

This is the polyline current code which I need to amend from the others:

if (geometry.type === "polygon" || geometry.type === "multipoint" || geometry.type === "polyline") {
        var searchExtent = geometry.getExtent();
        var biggerExtent = new esri.geometry.Extent(searchExtent.xmin - 1000, searchExtent.ymin - 1000, searchExtent.xmax + 1000, searchExtent.ymax + 1000, new esri.SpatialReference({
                    wkid : 27700
                }));
        graphic1 = new esri.Graphic(biggerExtent, buffersymbol);
        buffergraphicsLayer.add(graphic1);
        bufferExtent = getBufferExtent();
        graphic2 = new esri.Graphic(bufferExtent, bufferextentsymbol);
        buffergraphicsLayer.add(graphic2);


texts: How to update raster metadata in postgis

texts:
How can I update the raster metadata in postgis?

 upperleftx | upperlefty | width | height | scalex | scaley | skewx | skewy | srid | numbands
------------+------------+-------+--------+--------+--------+-------+-------+------+----------
          0 |          0 |   600 |    900 |      1 |     -1 |     0 |     0 | 4236 |        1

to

 upperleftx | upperlefty | width | height | scalex | scaley | skewx | skewy | srid | numbands
------------+------------+-------+--------+--------+--------+-------+-------+------+----------
        137 |         56 |   600 |    900 |      1 |     -1 |     0 |     0 | 4236 |        1

I tried:

UPDATE vol SET (ST_Metadata(rast)).upperleftx=137
    FROM (SELECT (ST_Metadata(rast)).upperleftx FROM vol) as foo;

or

UPDATE vol SET (ST_Metadata(rast)).upperleftx=137

It is an idea, but this doesn't work.


texts: Retrieve geojson property and make layer visible but only based on some values of properties (OL3)

texts:
Please take a look at code: https://plnkr.co/edit/wF8VZnnLtRLHPkfrKGu5?p=preview
(OpenLayers v.3.8.2)
I would like to retrieve GeoJSON property "zoom".
I guess it would be something like this but dont know where to put it in a code:
var zoom = function(feature){
    var zoomGeojson = feature.get("zoom");
    console.log(zoomGeojson);
}

And then I would like to compare zoom level of the map to the GeoJSON property "zoom". If they would match, I would like to visualize only that geojson data. (for example: mapzoomlevel = 5, if mapzoomlevel = geojsonpropertyzoom visualize only data from geojson that contains property "zoom"=5)
Edit
Based on @Hicham Zouarhi answer.
I ran into two problems.
1.Code below throws me an error:
var tmpLayer= new ol.layer.Vector();
map.addLayer(tmpLayer);

source.forEachFeature(function(feature){
    var mapzoom = map.getView().getZoom();
    var geojsonzoom = feature.get("zoom");
    console.log('MapZoom ' + mapzoom);
    console.log('GeoJSONZoom ' + geojsonzoom);
    if(geojsonzoom == mapzoom){
        console.log('GeoJSON zoom property and map zoom are equal!');
        tmpLayer.getSource().addFeatures(feature); // error: Uncaught TypeError: Cannot read property 'addFeatures' of null(…)
    } else {
        console.log('Not equal!');
    }
});

Without that error it runs smoothly in a console, however it does not run in my code at all. I guess I am missing something or the order is not right.

I have reorganized code by breaking it into a smaller pieces: https://plnkr.co/edit/DeaVFGBK11jKCT3yHAVd


texts: How to (vertically) register a DEM to GPS readings?

texts:
I don't know if what I ask is possible or not but let's try...

I have two datasets:

a DEM as a point cloud in a relative coordinate system, constructed with 3D algorithm structure-from-motion.
20 ground control points with absolute elevations values (between 130 and 140 meters) taken with a GPS.
20 points of the point cloud representing the GCPs in the relative coordinate system (georeferencing targets visible on the point cloud).

I would like to constrain the elevation of my point cloud DEM with the GPS ground truth elevation, in order to have a DEM with absolute elevation calibrated to my study area.

At the moment, the study area is about 1 hectare and with slopes around max. 20% but the aim is to export the methodology in river catchment with a size of some kilometers-square and with max. slopes around 60%.

Any help is welcome, preferably in ArcGIS or R but I can deal with other softwares if necessary!

Thanks for the help.


texts: Label point features based on clock position?

texts:
I have recently been tasked with collecting data on the position in which sewer lines enter manholes for a local government Public Works department. To collect these points, field workers are using the collector app to drop a point on top of a manhole, and fill in a field that indicates the clock position that the lines are entering the manhole (Based on a North orientation).

So now that I have points that have a clock position field with values such as "3" if the line was coming in directly from the east, or "9" if the line entered from the west, I am wondering how I can label these points.

In essence, I am looking to assign label position based on the value of this clock position field.

Is this possible with ArcGIS Pro Maplex label engine?


texts: Adding Layers to a opened Qgis project via Python IDLE

texts:
I have a project opened in Qgis "qgis_training.qgs" and I am trying to add a layer to it via python2.7.6 IDLE using the following code:

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

# Get the project instance
project = QgsProject.instance()
# Print the current project file name (might be empty in case no projects have been loaded)
print project.fileName()

# Load another project
project.read(QFileInfo('/home/julierme/Desktop/qgis/projetos/qgis_training.qgs'))
print project.fileName()

bridge = QgsLayerTreeMapCanvasBridge( \
QgsProject.instance().layerTreeRoot(), canvas)

# Now you can safely load your project and see it in the canvas
project.read(QFileInfo('/home/julierme/Desktop/qgis/projetos/qgis_training.qgs'))

# Function addVectorLayer: Open and display a vector layer in QGIS
layer = iface.addVectorLayer("/home/julierme/Desktop/amazonia_sar/amapa/SHP/mask_covmingrad.shp", "mask_covmingrad", "ogr")
if not layer.IsValid():
    print "Layer failed to load!"

# Save the project to the same
project.write()

I am getting the error message below while  running it:

Traceback (most recent call last):
  File "/home/julierme/pyqgis/scripts/script2.py", line 17, in <module>
    bridge = QgsLayerTreeMapCanvasBridge( \
NameError: name 'QgsLayerTreeMapCanvasBridge' is not defined

How can I solve this issue and achieve my goal?


texts: Embedding Images in map layouts and documents (*.mxd)?

texts:
I wanted to place some images in my map layout but I found that images are not embedded in the MXD files so now, after I deleted them, I've lost the graphics and have to do them again.

Is here any way to embed the image in the MXD files itself? If not, where would you suggest I place the images; the geodatabase that comes with the map(s)?


texts: Loading Raster and vector data In MapGuide using Maestro

texts:
I am new to Mapguide , working on it for my summer Internship.I have Installed it perfectly.
Now the task is to load vector data and Raster data in map Guide (mainly have to load *.shp files and *.tif files). I have find not much resource that provide exact steps. Also I don't know where to put data files ? In Repository or somewhere else ? How to link these file ?
Can any one help me in this?


texts: I can't figure out how to batch deselect layers in ArcGIS Pro. I have added ~200 images to a project map, and now want to deselect without deleting all of them because I applied a common symbology. This was super easy in ArcMap, could select and right click, but that doesn't seem like an option in ArcGIS Pro.


texts: Preset min/max values in symbology when importing a raster into QGIS

texts: I have 200+ .tif files I need to import into QGIS and visually inspect. Due to the nature of the data, I need to manually set the min/max values to 0 and 0.8 , respectively. No matter what, the min/max values always need to be the same in order to view the data. Is there a way to automatically have QGIS import the rasters with those min/max values already set so I don't need to manually set each one in the symbology tab when added as a layer?


texts: Seeking database that identifies adjacent postal addresses and/or zoning of postal address?

texts:
Can someone point to:

a database that identifies adjacent postal addresses (for example, if I provide an address, the database provides address of houses on either side of house address)?
a database that lists zoning of addresses, so as to distinguish commercial, house, and apartment buildings?

I realize it's unlikely I'll get either of this data apart form large urban centers.

This is for a web app, and my target would be major US cities. I usually use a php-javascript stack, but I can also work with python. I would prefer a free database with api access over other options, but I'm open. The idea of the web app would be to collect user addresses and then compare those addresses (say to identify a neighbor).


texts: What's the command to open a project file and navigate to certain coordinates?

texts:
Is it possible to run a command line which opens QGIS, loads a particular project file and navigates to a set of coordinates?

If so then how would I do this?


texts:
I'm trying to change the default symbol (font, halo etc) for labels of a layer that I am adding to an mxd. Below is what I have so far, but the template layer font style didn't show up in the map. arcpy.mapping.updateLayer didn't work either.

citySym = r"C:\CitySymboly.lyr"
cityFc  = r"C:\\City"
cityMap = arcpy.mapping.Layer(cityFc)
#
#......other code not relevant
#
cityname = cityFc + subsetname
addCity = arcpy.mapping.Layer(cityname)
cityLayer=arcpy.mapping.AddLayer(df, addCity, "TOP")
inputCity = arcpy.mapping.ListLayers(mapdoc,cityLayer,df)[0]
#arcpy.mapping.UpdateLayer(df, inputCity, citySym, True)
arcpy.ApplySymbologyFromLayer_management (inputCity, citySym)
inputCity.showLabels = True


texts: Geoprocessing in CRS WGS84 and UTM Zone in QGIS

texts: Why do I can't do union or do another geoprocessing two shapefiles with different CRS in QGIS? For example, I want to do union two shapefiles, one shp with CRS WGS84 and another shp with CRS WGS84/UTM 47N. It can't be happened and a message box will appears. Even though this process can be done in ArcGIS. I'm using QGIS 2.14 Essen. Any suggest?


texts: Improve render speed for QGIS canvas (hardware)

texts:
I'm using QGIS to build complex maps from OSM-Data (GPKG or PostGIS). I'm wondering, how I should upgrade my system, to speed up rendering in canvas:

Intel I5 4x3,2GHz, 16GB DDR3, SSD

Or what about using RAM-Disk or RAID?

In the GIF-Animation you can see real time rendering at changing scale:
https://media.giphy.com/media/MUlt8lb1iuzm80alyE/giphy.gif
2 or more seconds for every canvas move or scale are just too much for a smooth work.

I already use parallel rendering for all CPU and simplification:


texts: Create conditional to check if raster has value

texts:
I have a problem with this code:

ac = ExtractByMask(raster, "area.shp")
arcpy.RasterToPoint_conversion(ac, "point.shp", "VALUE")

The problem is that sometimes the raster has no value.
(It's a flow accumulation and the precision of my raster is not enough for what I am asking but that's not the point).
So the RasterToPoint can't work.

So my question is:
How can I make a condition like: "if my raster has no value, don't do the RasterToPoint"


texts: How to merge connected lines with same direction (PostGis)

texts:
I am looking for a way to join all lines that are connected,  have the same direction, i.e. start point of one is end point of the other, and share some common attributes. In fact, I want to connect oneway streets having an common start/end point and the same oneway direction. I found a way to identify pairs of lines that fulfill these condition:

select a.osm_id, a.name,a.oneway,
       b.osm_id, b.name,b.oneway, from planet_osm_line a, planet_osm_line b
   WHERE (st_startpoint(a.way)=st_endpoint(b.way) or st_endpoint(a.way)=st_startpoint(b.way))
         and a.oneway=b.oneway ;

I am now looking for a way to merge the geometries of these lines.


texts: gdal.Warp -cutline option (geometry object instead of vector file)

texts:
I am new to ogr and vector data manipulation.

I am interested to cut an image (in this case a DEM file) using gdal.Warp() with the '-cutline' option to used the borders of a rectangular polygon to cut the underlying image. I saw good examples here

I wonder if it is possible to provide a "polygon object" instead of a "vector file"(.kml, .shp etc.) as the input to "cutlineDSName"?


texts:
Hi i've been setting up a script to assign values depending on one type of polygon intersecting with another polygon. I have two feature tables - one with building footprints and one with polygons that each have a specific value that needs to be tied to the building these polygons intersect with. So far it is set up by first creating a new table (named temp) with only the buildings that intersects with the polygons:

SELECT buildings.geom as geom_build, polygons.geom as geom_poly, polygons.value, buildings.id
into temp
From buildings, polygons
Where st_intersects(buildings.geom, polygons.geom);

and then the values are divided into sections using a case, that will later be tied to the buildings:

ALTER TABLE temp ADD column IF NOT EXISTS "interest" integer

and

UPDATE temp set interest =
(Case
    When temp.value <= 5 THEN 1
    When temp.value BETWEEN 6 AND 10 THEN 2
    When temp.value BETWEEN 11 AND 15 THEN 3
    When temp.value BETWEEN 16 AND 20 THEN 4
    When temp.value >= 21 then 5
End);

These 1, 2, 3, 4 and 5 are the values of interest and they are paired back to the buildings table by matching buildings id. There can be multiple polygons under one building, and it seems it is the value of the last (in the iteration) polygon that intersects with the building that is assigned.

My wish is to always assign the lowest value from however many polygons that may intersect to the building instead of getting the value that came from the last polygon in the iteration. Or in other words, if more than one condition in the case were to be true isolated from eachother, I want some condition returns to have higher priority than others.

For example if a building intersects with three polygons, where the first polygon has the value 2, the second has the value 7 and the third one has the value 12. If the third polygon was the last to go through the case, then the building would end up with the case return value '3', as the first two returns would be overwritten. Is there a way to make sure a return of '1' will never be overwritten by a higher number?

update:

I solved it by using the following sql lines:

Select id, MIN(interestvalue)
into tempTwo
from temp
GROUP BY id;

And then matching this min interest value back to the buildings using the ID.


texts: Searching for c# code to convert from UTM to WGS1984 (and back)

texts:
I'm looking for code in c# that converts coordinates from UTM to WGS1984 (and preferably backwards as well). I need to convert coordinates from only one zone (36N).

I've found this and this, but they throw me 150-200 meters from the real point (The allowed diversion I'm aiming  for is 5 meters tops).

EDIT: Found this as well, but did not have an opportunity to check this for accuracy.

thanks in advance.

EDIT: Found this as well, but did not have an opportunity to check this for accuracy.


texts: How to add a layer for the ocean from imported 'bil' files

texts: I am very new to QGIS. I downloaded data from earth explorer for east of florida. I am trying to create a decent landscape for my flight simulator "Condor" as I am a glider pilot. So, I merged the bil file, then warped it and then clipped it. I added various vector layers from QuickOSM. I downloaded the water from the Earth Explorer site and added the layer. But it doesn't show up when I deselect the clipped layer. I attached 2 pics, one of the clipped and the other without it. The one with the clipped layer has a black line that's the coast. And the other is the intercoastal water way. How do I add the small portion in the clipped for the ocean? Clipped view of bil file


texts: Counting overlapping polygons from one Shapefile using ArcGIS Desktop?

texts:
I have a shapefile with 16,400 polygons. Each polygon shows the expansion of a bird species for the whole world.

Now I have to count the overlapping polygons. I tried it with union and dissolve (count the union), but the union is not working for so many polygons.

Then I tried to clip continents, but this is also not working because of the huge number of polygons. Moreover I tried this method , also without success.

Therefore I'm asking you I there is a way to count overlapping polygons if 16400 polygons are in one shapefile?

I'm working with 10.0 and can work with 10.2. An ArcPy solution is also wonderful.

At the moment I am thinking about creating a fishnet and iterate over the rows of the shp with the 16400 polygons and write 1 to a value field of a fishnet cell if the polygon is in this cell and than take the next row (polygon) and if this is also in the fishnet cell count +1.

But I don't know if this is a good solution and how to realize it. Or I have to learn R to use this approach.

The result: It should be a shape where you have new polygons out of the overlapping ones and a field where the overlaps are counted.

So in the end there should be a shapefile where you can see how many bird species are found at the same place.


texts:
I am trying to request an WFS url within the circle. With the following url:

https://localhost:81?request=GetFeature&version=1.0.0&service=WFS&typename=ferrets:ab&outputformat=json&cql_filter=DWITHIN(geom,Point(-232855.14975877, 7787813.215242146),100,kilometers)

But I am getting an error:

Could not parse CQL filter list. Encountered "," at line 1, column 36. Was expecting one of: "-" ... <INTEGER_LITERAL> ... <FLOATING_LITERAL> ... Parsing : DWITHIN(geom,Point(-232855.14975877,7787813.215242146),100,kilometers).

The same I am trying to implement using the OpenLayers I am facing there error too.
function featureRequest(featureType, postalCode) {
        return new ol.format.WFS().writeGetFeature({
          srsName: 'EPSG:3857',
          featureNS: 'http://www.opengeospatial.net/ferrets',
          featurePrefix: 'ferrets',
          featureTypes: [featureType],
          outputFormat: 'application/json',
          filter: ol.format.filter.dwithin('geom', 'Pont(-232855.14975877, 7787813.215242146)', 1000, 'kilometers')
        });
    }

Using OpenLayers also I am not able to get jsondata. I am not sure what I am doing wrong here.


texts: Displaying symbology outside map canvas window extent?

texts:
Say we have equally distributed x,y mesh vector file.  Each coordinate has a map-unit based diameter circle symbology.  If we zoom into a region of the map canvas far enough, some centroids exist outside of the map canvas view, thus they will not be displayed.

Is there any means to set a buffer surrounding the map canvas extents equal to the map-circle symbol diameters to prevent border coordinates from being omitted from render?

(yellow highlighted region in image attached image below)

Furthermore, implementing such into a pyqgis plugin?

(Hopefully) creating such a buffer would ensure coordinates just outside the map-view would still maintain symbol display.

SOLUTION TEST.1:

The following snip was an (unsuccessful) adhoc test.  Likely as I am improperly using the QgsRenderContext() Class?

    (init method)
    ....
    self.cellDim=None
    self.iface.mapCanvas().extentsChanged.connect(self.updateMapBuffer)
    ....

def updateMapBuffer(self):
    try:
        if self.cellDem==None:
            self.aquireGridDim()
        xMinBuffer = self.iface.mapCanvas().extent().xMinimum()-(self.cellDim*10)
        xMaxBuffer = self.iface.mapCanvas().extent().xMaximum()+(self.cellDim*10)
        yMinBuffer = self.iface.mapCanvas().extent().yMinimum()-(self.cellDim*10)
        yMaxBuffer = self.iface.mapCanvas().extent().yMaximum()+(self.cellDim*10)
        newExtent = QgsRectangle(xMinBuffer,yMinBuffer,xMaxBuffer,yMaxBuffer)
        QgsRenderContext().setExtent(newExtent)
    except AttributeError: pass

Probably also important to note - this test gives miniDumps and subsequent QGIS crash 9/10 times when opening.  Obviously not the proper means.


texts: How to reproject spatial data using free libraries?

texts:
How can I use free libraries to transform spatial data?

For example,
I want to change the projection of a Shapefile within the code of my C# web application. How do I do that?


texts: extract elevation values along a line using Google Elevation API

texts:
I note that there is the Google Elevation API.

I see they have an example at the bottom showing an elevation profile for a line. Is there a premade service someone has put up on a website for me to submit a line and it will return elevation values at preset distances along the line?

I have downloaded the SRTM data for the area I'm interested in and can do this using ArcGIS, but it would be good to compare it to what Google Elevation API has as well.

I see someone has done this with Mathematica, but I have never used this program and don't understand how easy it would be for a beginner to produce what they have done on the blog. General comments about this are welcome.


texts: Im using arcmaps 10.1 and i have to imput pen points off a buch of old maps into arcmaps i have everything in order in arc maps but need help displaying the utm zones, basemaps, blocks, possibly sub-blocks as well. Having a shapefile with that all in it would be nice or if theres a simple way to get this to display in data view that would be appreciated.


texts: Creating a tile layer for a polygon for use in QGIS Atlas?

texts:
I want to create a series of plans at a set scale that covers a large area defined by a polygon.  The polygons are a buffer for a water network, so is very irregular with sections that are long and skinny.

Is there a way to create a set of rectangular tiles that covers the polygon, that I can then use the Atlas feature in QGIS 3,8. Ideally the rectangular tiles could be rotated to better accommodate the long skinny sections, but this is not critical.

A copy of the polygon I want to tile over is below.  There would be a large number of plans, as I want to generate the plans at say 1:20,000 scale to enable viewing of detail.  I would use the plans to print a number of plans with different themes, which I would periodically regenerate.


texts: Can you animate zooming without triggering navigation

texts:
In my current project I'm updating my app's GPS system from ArcGIS runtime 10.2 to 100.7. From testing, everything seems to work but the zooming is not smooth whatsoever. In fact currently when you set the AutoPanMode to Navigation it is horrendously choppy and borderline painful to watch it zoom to your position. I've created a test application with the 10.2 references and setting AutoPanMode to Navigation smoothly pans over to your location but never zooms unless you call some Zoom method afterwards. I'm guessing that's why it used to look so smooth?

Currently this is what I've found creates the same functionality

await MapView.SetViewpointAsync(new Viewpoint(e.Position.Extent), TimeSpan.FromSeconds(1));
await MapView.SetViewpointScaleAsync(300);

The problem is when I do this, the AutoPanMode gets set to Off and the only workaround is to set it back to Navigation right after

await CustomMapView.SetViewpointAsync(new Viewpoint(e.Position.Extent), TimeSpan.FromSeconds(1));
await CustomMapView.SetViewpointScaleAsync(300);
MapView.LocationDisplay.AutoPanMode = LocationDisplayAutoPanMode.Navigation;

This is not really ideal since it feels a little hacky. Additionally in UWP it causes a bunch of threading issues that crash the app (I can fix that but it looks really gross). Is there a way I can animate my zooming without having to do this?

What had to be done to fix this:

MapView.LocationDisplay.InitialZoomScale = 0;
MapView.LocationDisplay.IsEnabled = true;
MapView.LocationDisplay.AutoPanMode = LocationDisplayAutoPanMode.Navigation
await Task.Delay(2500);
await MapView.SetViewpointScaleAsync(500);


texts: How to produce spatial grid from raster?

texts: I need a Spatial Grid as a master grid for diverse thematic maps. How do I produce a Spatial Grid from a raster discarding all NA pixels?


texts: SLD filter questions using the other two conditions

texts:
I would like to specify the color and size in accordance with the speed and scale.
This is a step of respectively 10.

    <Rule>
      <Name>rotate_ps1-1</Name>
      <Title>arrow1</Title>
      <Abstract>arrow_rotate_ps</Abstract>
      <MinScaleDenominator>0</MinScaleDenominator>
      <MaxScaleDenominator>5000</MaxScaleDenominator>
      <ogc:Filter>
        <ogc:And>
          <ogc:PropertyIsGreaterThanOrEqualTo>
            <ogc:PropertyName>speed</ogc:PropertyName>
            <ogc:Literal>0.0</ogc:Literal>
          </ogc:PropertyIsGreaterThanOrEqualTo>
          <ogc:PropertyIsLessThan>
            <ogc:PropertyName>speed</ogc:PropertyName>
            <ogc:Literal>0.05</ogc:Literal>
          </ogc:PropertyIsLessThan>
        </ogc:And>
      </ogc:Filter>
      <PointSymbolizer>
        <Geometry>
          <ogc:Function name="startPoint">
            <ogc:PropertyName>geom</ogc:PropertyName>
          </ogc:Function>
        </Geometry>
        <Graphic>
          <Mark>
            <WellKnownName>ttf://Wingdings 3#139</WellKnownName>
            <Fill>
              <CssParameter name="fill">#5F00FF</CssParameter>
            </Fill>
          </Mark>
          <Size>100</Size>
          <Rotation>
            <ogc:PropertyName>direct</ogc:PropertyName>
          </Rotation>
        </Graphic>
      </PointSymbolizer>
    </Rule>

The code is a part of my sld code.
When in that way, my sld code is too long.
Is there any other way?


texts: Why is the dynamic text on data driven pages not displaying a value from the grid attribute table?

texts:
Relating to my question here, I've been trying to return a value from the attribute table using the following code:

  <dyn type="page" property="lu1p"/>

However, even if I try to use a different attribute heading (i.e. lu2p), it only displays "[empty]". Is there something which I'm missing?

The grid index page is actually a set of irregular study areas (i.e. not perfectly square), but it hasn't effected the use of Data Driven Pages. I have looked at the help files for this but have not found any solution.

Edit: As requested, here is a screenshot:

(http://i.imgur.com/8FXwS.jpg) available here also

Edit: I'm still having this issue so if anyone could provide an alternative answer, I'd really appreciate it. I tried again in a new mxd and it worked once but then it stopped worked when I re-opened it.


texts:
I am reclassifying 43 rasters in Grass r.reclass function and am getting the following error:

Processing algorithm 1/43…
Algorithm r.reclass starting…
Input parameters:
{'GRASS_RASTER_FORMAT_META': '',
'GRASS_RASTER_FORMAT_OPT': '',
'GRASS_REGION_CELLSIZE_PARAMETER': 0.0,
'GRASS_REGION_PARAMETER': None,
'input': 'NDWI_20160619',
'output': <QgsProcessingOutputLayerDefinition {'sink':C:/QGIS_Projects/Reclass_Results/reclass_NDWI_20160619.tif, 'createOptions': {}}>,
'rules': 'C:/Users/Steffan/Desktop/reclassrules.txt',
'txtrules': '-1 thru 0 = 0 LAND \n0 thru 1 = 1 WATER'}

g.region n=0.0 s=0.0 e=0.0 w=0.0 res=100.0
r.reclass input=rast_5c34701952bf13 rules="C:/Users/Steffan/AppData/Local/Temp/processing_5b7498a321e24a2e95e9076e4c4b05b3\5c3470487e14f5" output=outputac479b03fcc243a28f02487b7db35111 --overwrite
g.region raster=outputac479b03fcc243a28f02487b7db35111
r.out.gdal -t -m input="outputac479b03fcc243a28f02487b7db35111" output="C:\QGIS_Projects\Reclass_Results\reclass_NDWI_20160619.tif" format="GTiff" createopt="TFW=YES,COMPRESS=LZW" --overwrite
Starting GRASS GIS...
Traceback (most recent call last):
File "C:/PROGRA~1/QGIS3~1.4/apps/qgis/./python/plugins\processing\algs\grass7\Grass7Algorithm.py", line 413, in processAlgorithm
Grass7Utils.executeGrass(self.commands, feedback, self.outputCommands)
File "C:/PROGRA~1/QGIS3~1.4/apps/qgis/./python/plugins\processing\algs\grass7\Grass7Utils.py", line 386, in executeGrass
feedback.reportError(line.strip())
TypeError: reportError(self, error: str, fatalError: bool): not enough arguments

Batch execution completed in 0.23 seconds


texts:
Folium is the python wrapper of leaflet.

It is possible in folium as well in leaflet to use markercluster when the number of gps grows.
But markercluster have a limit of 50k gps , recently i found that the owner of leaflet have created another clustring algorithme for huge spatial data called supercluster that can handle 6 millions of gps.

my question is: Is it possible to use the leaflet supercluster on folium maps ?  if yes how to use it on folium  ?


texts: Setting No Data value in a TIF using GDAL

texts:
I am working with various TIF-files using a Python Script and GDAL Library.
Calculations within the array work fine and TIF-output worked too. I just want to Set the NoData-value by getting it from another TIF with correct NoData-Value but NoData-Value doesen't change.
May the input-files already be corrupt or is there an error in my code that I don't see?

from osgeo import gdal

intwifile = 'TWI_normed.tif' # File where NoData-Value needs to be Set
inndvfile = 'TWI.tif' # File, where NoData-Value is to be picked from
corrndTWI = 'NDVCorr_TWI.tif' # Output-file

driver = gdal.GetDriverByName("GTiff")

ndraster = gdal.Open(inndvfile)  #read raster to get ND-value
ndband = ndraster.GetRasterBand(1) # read band to pick ND-value
ndval = ndband.GetNoDataValue() # pick ND-Value

ras = gdal.Open(intwifile) # read raster where ND-value needs to be changed
#print(type(ras))
xlen, ylen = ras.RasterXSize, ras.RasterYSize # get cols and rows to write output file
band = ras.GetRasterBand(1) # get band information to write in output file

dtype = band.DataType # get data type for output-file

ndcorrected = driver.Create(corrndTWI, xlen, ylen, 1, dtype) # create output raster
data = band.ReadAsArray() # read input band as Array
transform = ras.GetGeoTransform() # get transformation parameters from input layer
crs = ras.GetProjection() # get coordinate reference from input layer
ndcorrected.SetGeoTransform(transform) # set transform for output
ndcorrected.SetProjection(crs) # set crs for output
ndcorrected.GetRasterBand(1).WriteArray(data) #write data from input to output
ndcorrected.GetRasterBand(1).SetNoDataValue(ndval) # set no data value
ndcorrected.FlushCache()


texts: points2one plugin subsitute for QGIS 3.4.4 version

texts: Points2one plugin is not compatible with QGIS 3.4.4. It's quite useful for our works related to forestry. Is there any other substitute plugin or new upgrade version of points2one plugin?


texts: ArcGIS - Is it possible to bypass an inaccurate WMS bounding box provided by server metadata?

texts:
In ArcGIS Desktop, is there a way to bypass / ignore an inaccurate WMS bounding box provided by server metadata?

Indeed, for a specific WMS layer I am using, the bounding box provided by server metadata does not match the underlying data's extent (some data exist out of the "official" bounding box delivered by metadata request). Thus, I would like to prevent ArcGIS from using the inaccurate metadata so that I can see all of the data.


texts: Why can't I create a feature class in ArcCatalog when using a .MDB file?

texts:
I am trying to use an access database with geographic information in ArcCatalog. I created the database as a .accdb and backsaved to a .mdb format. I can create a 'shapefile', but not a 'File and Personal Geodatabase feature classes'. The error I get is:

  There was an error trying to process this table.

  The table was not found. [GDB_ObjectClasses]

  The Microsoft Jet database engine cannot find the input table or
  query 'GDB_ObjectClasses'.  Make sure it exists and that its name is
  spelled correctly.

I can bring the table into arc, create a shapefile, but when I try to use any tools (specifically kriging) I get a similar error.

I don't see any connections, relationships, anything in the access table itself. I've tried bringing the data out into Excel, copying and pasting in a new excel with no formatting and then bringing it back into access and I get the same errors. I don't think it is the dataset because of this. I am using 10.5.1. Any suggestions?


texts:
I load a map with data and collect information in a JSON file.  I want a tooltip appear when mousing over a town on the map.

   var vecLayer = new OpenLayers.Layer.Vector("vector", {
        styleMap: styleMap,
        protocol: new OpenLayers.Protocol.HTTP({
            // url: "app/data/summits.json",
            //url: "app/data/chargeIndicateur.php",
            url: "../web/app_dev.php/json/indicateur/"+gIndicateur+"/"+gEchelle+"/"+gPeriode,
            format: new OpenLayers.Format.GeoJSON(),
            idProperty: 'Issam'

        }),

, {
// dock: 'right',
xtype: 'toolbar',
id    : 'teste2',

items: [{
text: 'Nom',
labelStyle: 'font-weight:bold;',
}]
} ,


texts:
I want to show a number of GPX files in a TileMill map. Here's the process that I thought would work:

gpsbabel -i gpx -f "$myfile" -o kml,lines=1 -F "$myfile.kml"
Add as layer in Tilemill, style with line-width: 2; and voila.

But no. Only points are being rendered, not the lines between them.

The generated .kml files contain <linestring> elements, which I would have thought is the key:

    <Placemark>
      <name>Path</name>
      <styleUrl>#lineStyle</styleUrl>
      <LineString>
        <tessellate>1</tessellate>
        <coordinates>
          146.153490,-38.211270

Is the problem a shortcoming at the TileMill end, or something in how I'm converting the files? Unfortunately none of the other formats supported by TileMill (GeoJSON, Shapefiles...) are supported by GPSBabel.

Other options? I'm using PostGIS, so a way to load them straight into the database would also be fine.


texts: How can you map counties to congressional districts?

texts: Given county A, I need to know what congressional districts are within it.


texts: Ideal file format for QGIS that allows mixed geometry types?

texts:
I have a spreadsheet for an archaeological project that is hosted on Google docs. This allows multiple contributors to edit and update the table as revisions and field data comes in. The table is not very complex, just site attribute information and the spatial information in WKT. Some of the sites are points and some are polylines.

Currently, when I want to map the data I export the table from google docs as a csv and then import into QGIS. In the csv import dialog one must specify the geometry type, so I must run the import process twice for each geometry type and then end up with two layers in QGIS, one for point features and the other for poylines. Although this process works well enough, it gets confusing for spatial operations and analysis because the project is now split across two separate layers and tables.

Is there some way to import data in QGIS from a csv with WKT and transform it into a file type that will allow the data to remain organized by attribute (ie project), and not by geometry/structure?

I've tried this in PostGIS/postgresql, but it appears that the project data would be split into two tables based on geometry, if I use QGIS to import the csv and then export to the PostGIS database. I've tried to import a csv without spatial geometry set into QGIS and then export to PostGIS but can't seem to tell the database to then recognize the WKT field. At any event, I have no idea if PostGIS is even a solution to this problem.


texts:
I created a layer to visualize cable-laying and want to export the sum of the cable-lengths to an excel file when I'm done so i created an attribute for it.

How do I link an attribute to the auto-calculated $length?

I am using QGIS 3.0


texts: Is it recommended to change Geoserver admin username?

texts:
The docs and a bunch of other questions on here are circulated around changing the default password for admin on GeoServer. I know it's highly recommended to change the default password for the admin account for security reasons, but what about the admin username?

Is it standard practice to change the admin username for security purposes?


texts: QgsVectorLayerCache causes QGIS to crash when called within def run(self)

texts:
This is in reference to the similar question regarding QgsVectorLayerCache .
The original code is taken from here

I'm trying to implement this method within the def run(self) method.

 def run(self):

    new_dialog = QDialog()

    new_dialog.resize(270, 210)

    layer = iface.activeLayer()
    canvas = iface.mapCanvas()
    feature_list_view = QgsFeatureListView(new_dialog)

    vector_layer_cache = QgsVectorLayerCache(layer, 1000)

    # for layer in QgsMapLayerRegistry.instance().mapLayers().values():
    #     if layer.name() == 'Pins':
    #         # pins_layer = layer
    #         self.pinLayer = layer
    #         self.provider = self.pinLayer.dataProvider()

    attribute_table_model = QgsAttributeTableModel(vector_layer_cache)

    attribute_table_model.loadLayer()
    attribute_table_filter_model = QgsAttributeTableFilterModel(
        canvas,
        attribute_table_model
    )
    feature_list_model = QgsFeatureListModel(attribute_table_filter_model)
    feature_list_model.setDisplayExpression('"Pins"')
    feature_list_view_delegate = QgsFeatureListViewDelegate(feature_list_model)
    feature_list_view.setModel(feature_list_model)
    # feature_list_view.setItemDelegate(feature_list_view_delegate) # Crash if uncommented. TODO: fix issue

    new_dialog.show()

Debugging it with PyCharm has revealed that this line here causes QGIS to crash:

vector_layer_cache = QgsVectorLayerCache(layer, 1000)

How can I put this in a proper function without causing a crash?
QgsVectorLayerCache should still be availabe in QGIS 2.18 according to the  API docu


texts: Raster Calculator Workaround in Arcpy for Suitabiliy Analysis

texts:
This is my first attempt at a complex stand alone python code. I am self teaching the arcpy module, well at least attempting to.

I am having a great deal of difficulty getting my site suitability analysis code to run. It consistently crashes at the raster calculator step. I understand that the raster calculator tool is difficult to manage in arcpy, however I am equally unfamiliar with the method using map algebra.

The outputs/inputs are rasters. Everything is converted to rasters in prior steps. I get an error code: 000989 Python Syntax error: parsing error invalid syntax, referencing the line of code for raster calculator.

My code is included below. Can someone help me with a workaround for this issue?

# -*- coding: utf-8 -*-
# ---------------------------------------------------------------------------
# Megan E Brown
# ModelCode.py
# Created on: 2014-04-21 23:51:03.00000
#   (generated by ArcGIS/ModelBuilder)
# Usage: ModelCode <Scratch_Workspace> <Trails> <Trail_Proximity> <DEM> <Min_Elevation__m_> <Max_Slope____> <Final_raster>
# ---------------------------------------------------------------------------

# Import arcpy module
import arcpy, arcview
#from arcpy.sa import *
#env.workspace = "C:\\temp"

# Check out any necessary licenses
arcpy.CheckOutExtension("3D")
arcpy.CheckOutExtension("spatial")

# Script arguments
Scratch_Workspace = arcpy.GetParameterAsText(0)
if Scratch_Workspace == '#' or not Scratch_Workspace:
    Scratch_Workspace = "C:\\temp" # provide a default value if unspecified

Trails = arcpy.GetParameterAsText(1)
if Trails == '#' or not Trails:
    Trails = "Z:\\GIS500\\FinalProject\\FPData\\WELP_trails.shp" # provide a default value if unspecified

Trail_Proximity = arcpy.GetParameterAsText(2)
if Trail_Proximity == '#' or not Trail_Proximity:
    Trail_Proximity = "200 Feet" # provide a default value if unspecified

DEM = arcpy.GetParameterAsText(3)
if DEM == '#' or not DEM:
    DEM = "Z:\\GIS500\\ned\\ned" # provide a default value if unspecified

Min_Elevation__m_ = arcpy.GetParameterAsText(4)
if Min_Elevation__m_ == '#' or not Min_Elevation__m_:
    Min_Elevation__m_ = "170" # provide a default value if unspecified

Max_Slope____ = arcpy.GetParameterAsText(5)
if Max_Slope____ == '#' or not Max_Slope____:
    Max_Slope____ = "10" # provide a default value if unspecified

Final_raster = arcpy.GetParameterAsText(6)
if Final_raster == '#' or not Final_raster:
    Final_raster = "Z:\\GIS500\\FinalProject\\FPData\\tempdata\\yay" # provide a default value if unspecified

# Local variables:
Percent_Slope = "C:\\temp\\Percent_Slope"            #slope output
Slope_LE = "C:\\temp\\Slope_LE"                      #"output of less than equal"
Calc = "C:\\temp\\Calc"                              #output from raster calculator, input to final step Equal to
Elevation_GE = "C:\\temp\\Elevation_GE"              #"greater than equal output"
Buff = "C:\\temp\\Buff"                              #"output of buffer"
buffer_raster = "C:\\temp\\buffer_raster"            #"output of feature to raster making buffer into polys"
Near_trails = "C:\\temp\\Near_trails"                #"output from not equal reclassing the buffer raster for raster calculation"
Cell_size_30m = 30
Value_0 = "0"
Value3 = "3"

# Process: Slope
arcpy.Slope_3d(DEM, Percent_Slope, "PERCENT_RISE", "9.0021749575298E-06")

# Process: Less Than Equal
arcpy.gp.LessThanEqual_sa(Percent_Slope, Max_Slope____, Slope_LE)

# Process: Greater Than Equal
arcpy.gp.GreaterThanEqual_sa(DEM, Min_Elevation__m_, Elevation_GE)

# Process: Buffer
arcpy.Buffer_analysis(Trails, Buff, Trail_Proximity, "FULL", "ROUND", "NONE", "")

# Process: Feature to Raster
arcpy.FeatureToRaster_conversion(Buff + ".shp", "FID", buffer_raster, Cell_size_30m)

# Process: Not Equal
arcpy.gp.NotEqual_sa(buffer_raster, Value_0, Near_trails)

# Process: Raster Calculator
arcpy.gp.RasterCalculator_sa(((Slope_LE) + (Elevation_GE) + (Near_trails)), Calc)

# Process: Equal To
arcpy.gp.EqualTo_sa(Calc, Value3, Final_raster)


texts: How to Make my Legend Service for ArGIS Server 9.3 Map Service SSL?

texts:
I am developing an application using map services served by ArcGIS Server 9.3. The application must be completely https since it is shared through a secure Intranet. The application itself, all of the javascript libraries, and all of the map services I call are https. The browser I have to use is IE 11 and I am not able to upgrade to Server 10.

The map services do not load because I receive the network error that:
http://utility.arcgis.com/sharing/tools/legend?soapUrl=..... is not SSL enabled

I've previously read the following,

  "If the layers are version 10 or lower the legend is created using the ArcGIS.com legend service. In order to use the ArcGIS.com legend service your map service needs to be publicly accessible and your application must be able to access ArcGIS.com."

Is there any way around this? Also note - I do not actually need a legend - so even if there is a way to disable the legend service to be called, that would work!!


texts:
I successfully used the COGO traverse tool to create a polygon for a piece of property described in a legal description. My issue arose because it did not line up with data I already have of the area. I did not expect it to line up perfectly, but I did expect it to be closer than it was. The projections for all of the data in the map are the same so I'm not sure if it is because there is a specific projection used for legal descriptions or what.

I am using ArcMap 10.2.

The data I already have for the area is vast, but I am using a buildings feature class (polygon features), and a quarter section feature class (polygon features) from the datasets in our database. I also have a georeferenced image in my map that shows the property boundary surrounding the particular apartments I am looking at. This image is based on the legal description and comes from the engineers who worked on the apartment complex. Again, all data is set to the same projection.

My steps include:

Creating a polyline shapefile and a polygon shapefile.
In my editor toolbar I went to editing options and changed the "Direction Type" to  Quadrant Bearing, and the "Direction Units" to Degrees Minutes Seconds. I checked the "Ground to Grid Connection" box, then I made the "Distance Factor" 1 since my projection is in State Plane US Feet.
I started an editing session and began using the COGO traverse tool to create a polyline feature that according to the legal description commences at the corner of section 28 (remember I have quarter section map, which has the referenced section) and ends at the "True Point of Beginning."
From the True Point of Beginning I use the COGO traverse tool again to create a polygon for the property boundary based on the legal description.

If more details about my steps are still needed please let me know, but for now I am trying to keep it simple to avoid boring or dissuading the reader from reading.

Also, for information purposes, this project was given to me to see if it was possible to use GIS to map Legal Descriptions. So far I can't report yay or nay because yes I successfully mapped using a legal description, but no it did not end up where it should. It is very very close, but the boundary I COGO'ed does not cover all the apartments in my building feature class that it needs to. It just needs to move to the south about 100 ft. I am having a hard time believing the legal description is right and our building feature class is wrong because that would mean all the data we have in our database is wrong. I am also not convinced that the legal description is wrong. So that leaves me with thinking I made a technical user error. I can't make any recommendations until I figure out what the issue is and whether it can be fixed.

This is a screenshot of what I have been working on. The polygon is supposed to surround the smaller buildings in the middle.

This is a screenshot of the legal description I have been using.


texts: How to return feature.shape from server by using WCF and ArcGIS Silverlight API?

texts:
I am using the ArcGIS Silverlight API. I all I need to do is to get the geometry from a feature by objectid to zoom to it:

getFeature(objectid);

So I used WCF service to get the geometry but I received the following error:

Type 'System.__ComObject' cannot be serialized

Can anyone help?


texts: ArcGIS for Javascript Template available?

texts:
ESRI has a stack of templates in the ArcGIS Online section. But when you download them it turns out that in order to make them run, you have to first create a webmap via their WebMap Viewer application. And even then the map that you then host only references the webmap application you create on the ESRI Portal.

Coming from a FlexViewer background where I could use the Viewer app to create an application and then all files and folders were generated locally and I could edit or tweak as needed I was hoping for something similar with the WebMap Viewer.

I am looking for a blank webmap using the ArcGIS API for Javascript complete with widgets, place to insert basemaps, operational maps, etc. Does such a map template exist?


texts: Create GeoServer SQL View with variable number of arguments ( double precison [ ] )

texts:
I was trying to create a PostGIS procedural function which will perform a particular operation depending upon the user's provided coordinates. But since the user can provide any number of coordinates, I had to make the procedural function for all the possibilities (a painful approach). Example SQL query is below

select geom from test_4points('ways', 28.94603, 41.00764, 28.95402, 41.01789, 28.96145, 41.00522, 28.96650, 41.01191) -- Four points SQL
select geom from test_5points('ways', 28.94603, 41.00764, 28.95402, 41.01789, 28.96145, 41.00522, 28.96650, 41.01191, 28.96650, 41.01191) -- Five points SQL
select geom from test_6points('ways', 28.94603, 41.00764, 28.95402, 41.01789, 28.96145, 41.00522, 28.96650, 41.01191, 28.96650, 41.01191, 28.96650, 41.01191) -- Six points SQL

And in order to deal with each of these functions I had to create SQL View in my Geoserver for each one of them, like the following:

SELECT geom FROM test_4points('ways', %x1%, %y1%,%x2%,%y2%,%x3%,%y3%,%x4%,%y4%) ORDER BY seq -- for Four points
SELECT geom FROM test_5points('ways', %x1%, %y1%,%x2%,%y2%,%x3%,%y3%,%x4%,%y4%,%x5%,%y5%) ORDER BY seq -- for Five points
SELECT geom FROM test_6points('ways', %x1%, %y1%,%x2%,%y2%,%x3%,%y3%,%x4%,%y4%,%x5%,%y5%,%x6%,%y6%) ORDER BY seq -- for Six points

But now, I have create a modified SQL procedural function which can accept any number of coordinates, like the following:

select geom from test_allpoints('ways', 28.94603,41.00764,28.95402,41.01789,28.96145,41.00522,28.96650,41.01191) -- Four Points SQL, Notice that same function "test_allpoints" is taking any number of pts
select geom from test_allpoints('ways', 28.94603,41.00764,28.95402,41.01789,28.96145,41.00522,28.96155,41.01715,28.96650,41.01191) -- Five Points SQL
select geom from test_allpoints('ways', 28.94603,41.00764,28.95402,41.01789,28.96145,41.00522,28.96155,41.01715,28.94309,41.02114,28.96650,41.01191) -- Six Points SQL

How to create an SQL view in Geoserver in order to interact with this SQL 'test_allpoints()' functions query? Since there is no limit of the users provided coordinates, how to declare the arguments in the SQL View.

The "test_allpoints()" pl/pgsql function is of the following type:

test_allpoints(character varying, double precision[])

UPDATE1:
An elaborate version of my test_allpoints() plpgsql function is here, but my problem is not related with the pl/pgsql fun but the geoserver.

Regards!


texts: Converting each cell value of raster into floating-point representation using QGIS?

texts:
I have several bands of raster from the Landsat8 satellite, and I need to convert them into a floating-point representation. In ArcMap this is done through:

Spatial Analyst Toolbox -> Math -> Trigonometric -> Float.

How do I accomplish this in QGIS?


texts: Command line tool for reading .shp, .prj and converting to Google Maps lat lngs

texts:
I have a data file exported from an ARCGIS system (parcel data from local government).  Its a zip file containing 7 files:  f.dbf, f.prj, f.sbn, f.sbx, f.shp, f.shp.xml, f.shx.

I have a simple ruby script that uses the rgeo and rgeo/shapefile gem to read the .shp file and extract points.  The points don't look like lat, long eg.

936934.017177403 695547.1331697404

I would like to convert the points to lat longs that I can upload to google maps or similar.

I think the solution involves reading the .prj and .shp file, then telling the tool to project into EPSG:3867.

Is there a simple script or commandline tool that handles conversions like this.  eg:

convert f.shp f.prj -EPSG:3857 > new.shp

I would prefer to use a ruby / rgeo solution.  But anything unixy that can be scripted would work.


texts: Print Composer QGIS 2.6, how to export image?

texts:
When I finish map in print composer and click on Export to an image format and after I go to folder where I saved it there is no image.
If do same with Export as PDF ,I get PDF file but can't open it and it says that file is 0 Mb in size and every time I do export I lose all stuff in composer so I have to do all over again. Is this error in QGIS or something else?

How do I get JPEG image file from my map out of the composer?


texts:
I need to get MODIS NDVI (MOD13A3), but output from EarthExplorer or Daac2Disk HDF format is a numerical range is -2000 to 9996.

But I need for futher work basic range -1 to 1.

Is there any conversion method?


texts: Printing how many features satisfy a QgsFeatureRequest (QGIS 3)

texts:
I have a QgsFeatureRequest() that filter all the features by an attribute.
Is there any python function that could tell me how many features were filtered?
I know I could create a list with this features and then count them. But I want a better way.


texts: ArcGIS Network Analyst tool for Android?

texts:
I am trying to include the Network Analyst tool into an Android application.
I have published the tool on ArcGIS Server and I integrate its URL address to use it.
On the same time, I import a MapServer which display my network with all the layers needed to make NA work.

I tried to create a GeoProcessor and to input the tool in it and then execute it, but it doesn't work and I haven't found any clues on the web...

Does anyone has an idea about how to deal with this ?


texts: Performing Geoprocessing SQL Queries in QGIS from a PostGIS Database

texts:
This may be a dumb question, but is there way to run a PostGIS geoprocessing SQL query in QGIS? Specifically, I am thinking of running SQLs such as ST_Union, ST_Split etc etc.

I have tried the following SQL in the QGIS PostGIS Query Editor without any luck, while the query worked in PostGIS.

SELECT covtype, ST_Multi(ST_Union(f.the_geom)) as sGeom
INTO testLCC
FROM lcc2000_001k AS f
GROUP BY covtype;

SELECT populate_geometry_columns
('testLCC'::regclass);

ALTER TABLE testLCC ADD CONSTRAINT pkid PRIMARY KEY(covtype);


texts: Combining non-overlapping binary rasters with ArcGIS Raster Calculator?

texts:
The issues I keep seeing depend on unique combinations of values.

I have three rocktype raster layers for area and they are edge-to-edge (they do not overlap). They each have "VALUE" fields and "COUNT" fields. There's one row of data for each attribute table where "VALUE"=1 where the rock type exists. I guess this is like mask.

What I need to do is combine these three rasters into one with raster calculator using this logic.

rocktype1==1, output value=1

rocktype2==1, output value=2

rocktype3==1, output value=3

Here's what I've tried:

Con("rocktype2"==1, 2) | Con("rocktype1"==1, 1) | Con("rocktype3"==1, 3)
Con("paleo"==1, 2, Con("precamb"==1, 1, Con("tert_quat"==1, 3)))
Combine(["rocktype2","rocktype1","rocktype3"])
rocktype1 + rocktype2 + rocktype3

I continuously get a blank output raster. What am I doing wrong here?

I'm using ArcMap 10.3.


texts: How to "compile" ways and relations from OSM Planet

texts:
I need to get political boundaries, house polygons and some other data from OSM. To do that I downloaded OSM Planet dump, converted it into .pbf format (cause its fast and doesn't take a lot of space). Then, created a streamer that imported nodes, ways and relations into MongoDB.

Now I need to compile ways from nodes (in each way instead of reference to nodes, put their coordinates) and then compile relations from them both.

I tried to do that using Redis (key-value in-memory database) as compact as I could (id as key + lat,lon as value). But couldn't even import nodes - it stopped on 14% (I have 128 GB RAM)

My next try was to do the same thing with MongoDB but it's too slow.

How can I achieve this?
Are there any tools like osmconvert or osmfilter?

I'm sure people do that somehow cause there are a lot of projects based on OSM, but I have no idea how.


texts: Why are font and line sizes different in QGIS's Print Composer and how do I normalize them?

texts:
Why does my map look so different when I look at it in the QGIS print composer? The widths of lines change, the size of markers change, the size of fonts change.

Is there any way to have the two always in sync? My goal is generated a printed map so it's frustrating to have to load the Print Composer and then go back and tweak everything in the regular view. Everything is always smaller in the editor view and bigger in the print composer view.

I feel like there's something obvious I'm missing or not getting.


texts: Is there a application, or QGIS plugin to generate database diagrams (specifically for PostGIS) similar to the ESRI Diagrammer?


texts: GUI for managing GeoPackage (and other geodata forms)?

texts: I know that the browser in QGIS has quite a bit of capability for managing GeoPackages and other geodata, but I was wondering if there was a lightweight, standalone tool similar to ArcCatalog?


texts: This concerns post-processing via base station correction in Pathfinder Office. I understand the base station is a known location and the GPS receiver's location is corrected relative to the base. I'm trying to better understand how the post-processing process works. Does the correction application use the base station's location data to do a bias correction of the uncorrected data or does it redo the trilateration somehow?


texts: How to get elevation from Google Earth in Python

texts:
I'm new with the Google Earth Engine. I'm trying to get the elevation of 3 nodes in a list [[lan, lat],[lan, lat],[lan, lat]] and save it in a variable elevation. Can anybody show me how or suggest a good answer?
dataset = ee.Image('CGIAR/SRTM90_V4');
elevation = dataset.select('elevation');
slope = ee.Terrain.slope(elevation);
elevation_list=[....some long tat values]
elevatioval=Map.setCenter(elevation_list);

the output should be like this
[evation1,elevation2,elevation3...]


texts: How to burn a different value for each polygon in a json file using gdal.RasterizeLayer?

texts:
Is there a way to rasterize each polygon instance with its own burn value? So what I am trying to do is create a mask of the polygons in mask.shp using the projection and geotransform information from raster.tif. Each polygon mask needs to have a unique value. My current approach burns all of them as 255. Here is my code:

vector_layer = "mask.shp"
raster_layer = "raster.tif"
target_layer = "mask.tif"

# open the raster layer and get its relevant properties
raster_ds = gdal.Open(raster_layer, gdal.GA_ReadOnly)
xSize = raster_ds.RasterXSize
ySize = raster_ds.RasterYSize
geotransform = raster_ds.GetGeoTransform()
projection = raster_ds.GetProjection()

# create the target layer (1 band)
driver = gdal.GetDriverByName('GTiff')
target_ds = driver.Create(target_layer, xSize, ySize, 1, gdal.GDT_Byte)
target_ds.SetGeoTransform(geotransform)
target_ds.SetProjection(projection)
source_ds = ogr.Open(vector_layer)
source_layer = source_ds.GetLayer()

ds = gdal.RasterizeLayer(target_ds, [1], source_layer, burn_values=[255])

target_ds = 0


texts: Finding points where another point layer does NOT intersect a polygon layer

texts: I am using PostGreSQL to deal with 3 layers. I have two point layers(points A and points B) and a polygon layer. Some of the points in both A & B fall within the same polygon. I want to find all of the points in points A where points B do not intersect with the polygon layer. How can I intersect B and the polygons while only needing the points from point layer A?


texts: Where to put the save command in a Python with loop and if?

texts:
I have code to change the path of a shapefile with a specific name.
When it finds the same name, the path of the shapefile is changed.

This is the important part of the code:

for mxdFile in arcpy.ListFiles("*.mxd"):
    mxdPath = env.workspace + "\\" + mxdFile
    mxd = arcpy.mapping.MapDocument(mxdPath)
    dfList = arcpy.mapping.ListDataFrames(mxd, "*")
    print mxdPath
    for df in dfList:
        for lyr in arcpy.mapping.ListLayers(mxd, "", df):
            if lyr.name == layerfile:
                lyr.replaceDataSource(Folder, "SHAPEFILE_WORKSPACE", NewShapefile)
                print mxdPath," - changed"
    mxd.save()
del mxd

My problem is that the script saves all the MXD in the same folder of the .py file. Not just the ones it changed.

If I put the mxd.save() under the same align of the if lyr.name == layerfile: it saves just the ones it changes. But, there is a but....

But, the script stops saving through the looping for the dataframes. It changes the first one, then saves. And then, it seems to change in all of DF, because at the it really prints mxdPath," - changed" for each df. But the save just happens at the first DF changing...

How can I fix this?


texts: QGIS join by location not returning for features

texts:
I have eight contours/features in one shapefile (contour.shp) that i want to find the number of population within each contour, this is in another shp.

I am using join attributes by location to do the population counts, both shp's have the same CRS.

The results are fine for all but one of the contours/features. For that that contour a feature is not returned at all in the exported shapefile (not even a null), leaving only seven contours/features with counts.

I really need the counts in that 8th contours/feature, is this a problem with input shapefile (contour.shp)? if so what can solve the problem?

Keep all records was selected. The shapefile is polygon (I think). What are the implications either way.


texts: Trying to Install on MAC OS Lion and getting Python error message

texts:
Brand new to this and struggling to install QGIS on MAC OS X with Lion, I am not a techy - but noticed in the install read me that /usr/bin needs to be added to the GDAL.frameworks path.
Failing to see how and where to add this?
Full error message attached - any advice greatly appreciated and thanks for your time in advance
Couldn't load plugin GdalTools due an error when calling its classFactory() method

Traceback (most recent call last):
  File "/Applications/QGIS.app/Contents/MacOS/../Resources/python/qgis/utils.py", line 138, in startPlugin
    plugins[packageName] = package.classFactory(iface)
  File "/Applications/QGIS.app/Contents/MacOS/../Resources/python/plugins/GdalTools/init.py", line 32, in classFactory
    from GdalTools import GdalTools
  File "/Applications/QGIS.app/Contents/MacOS/../Resources/python/qgis/utils.py", line 283, in _import
    mod = _builtin_import(name, globals, locals, fromlist, level)
  File "/Applications/QGIS.app/Contents/MacOS/../Resources/python/plugins/GdalTools/GdalTools.py", line 39, in
    raise ImportError( error_str )
ImportError: No module named osgeo [python-gdal]

Python version:
2.7.1 (r271:86832, Jul 31 2011, 19:30:53)
[GCC 4.2.1 (Based on Apple Inc. build 5658) (LLVM build 2335.15.00)]

QGIS version:
1.7.4-Wroclaw Wroclaw,

Python path: ['/Applications/QGIS.app/Contents/MacOS/../Resources/python', '/Users/Rupert/.qgis/python', '/Users/Rupert/.qgis/python/plugins', '/Applications/QGIS.app/Contents/MacOS/../Resources/python/plugins', '/System/Library/Frameworks/Python.framework/Versions/2.7/lib/python27.zip', '/System/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7', '/System/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/plat-darwin', '/System/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/plat-mac', '/System/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/plat-mac/lib-scriptpackages', '/System/Library/Frameworks/Python.framework/Versions/2.7/Extras/lib/python', '/System/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/lib-tk', '/System/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/lib-old', '/System/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/lib-dynload', '/System/Library/Frameworks/Python.framework/Versions/2.7/Extras/lib/python/PyObjC', '/Applications/QGIS.app/Contents/Resources/python/plugins/fTools/tools']


texts: I'm fairly new to Mapserver and I was able to get US County boundaries to plot. Is there a way to get the County Name to plot in the center of the county? I've searched around and couldn't find any examples of doing this. Perhaps theres a better alternative that I haven't thought of. Thanks


texts: Database to get New York City frontages?

texts:
I'm looking for a database that gives the length of frontage of an address or gives the information necessary to calculate this for New York City (NYC) or for another large US city.

So far I haven't been able to find anything like this.

Does anyone know of a data set that would include this information?


texts: Editingt pixels of raster (geotiff) file using QGIS

texts:
I have a raster file which I want to edit. My raster is the output of a classification (on Google Earth Engine), but some pixels in the edges of some features are wrongly classified (inside the red polygon) and I want to erase them. See example:

I know I can create a mask drawing polygons to cover the undesired pixels, but this is a huge image with a lot of imperfections, and it would take forever to fix it.
Is there some technique with QGIS that allows to 'erase' these pixels as if it was photoshop using some brush?


texts: Copying non-spatial data from Oracle database to file geodatabase?

texts:
I have a model (shown below) that runs fine minus one small portion (outlined in red). I have tried various ways to remedy this but to no avail. The issue I'm having is this: The top portion of the model runs smoothly, as it is pulling data from an Oracle dB that has geometry (via the 'Make Query Layer' tool), and therefore I'm able to use the 'Copy Features' tool in order to store it in a File Geodatabase. This, as you may know, is because the 'Copy Features' tool only accepts a Feature Class or a Layer that has spatial attributes as input. Ergo, when it comes time to run the portion below it, I cannot use 'Copy Features', as it will not accept the 'Make Query Layer' of something from oracle that is not spatial, due to it simply being a table. So in its place, I am using 'Make Table View', which successfully accepts my 'Make Query Layer' as input. However, the error message seen below appears when it comes time to run itself. The rest of the model then takes the result of the 'Make Table View', and runs 'Make XY Event Layer' (as the table being pulled has an X/Y column for coordinates). From there, I then run the 'Copy Features' tool to import it into a File Geodatabase. Oddly enough, I also receive an error (listed below) if I try to simply go straight from the 'Make Query Layer' output and use it as the input for the 'Make XY Event Layer'. Trying the tool 'Copy Rows', also receives an error.

Is anyone aware of any other way in which you can pull non-spatial data (e.g. a table) and then somehow transfer it to a File Geodatabase Table?

I've searched through the other tools that seemed intuitive enough to look for, but had no such luck.


texts: Projection and transform lat,long coordinates in geo in OpenLayers

texts:
I am working on a project that from the ads-b data puts on the OpenLayers map the position of the intercepted planes ....
I encounter 2 problems:

The coordinates are in latitude and longitude and receive the following error:

Collection.js: 292 Uncaught e {code: 58, name: "AssertionError",
message: "Assertion failed. See
https://openlayers.org/en/v5.3.0/doc/errors/#58 for details."}

Two plane icons are displayed but there are many more ...

I think it is a problem with the projection of coordinates because the coordinates should somehow fall in Romania or at least near Romania and they fall in southern Africa ....

The code that takes json from ads-b:
$.post( "data.php",
                 {
                      get_data: 1
                 },
                 function( data ) {
                    var obj = JSON.parse(data);
                    var aircrafts = obj.acList;
                    //console.log(aircrafts);
                    for(var i = 0; i<aircrafts.length; i++ ){
                        if(aircrafts[i].Long || aircrafts[i].Lat){
                            var coordinates = [aircrafts[i].Long, aircrafts[i].Lat];
                            console.log(coordinates);
                            addMarker(coordinates);
                        }
                    }
                    })

The code with which I try to display on the map the icons at the position in latitude and longitude:
function addMarker(coordinates) {

            proj4.defs('EPSG:31700','+proj=sterea +lat_0=46 +lon_0=25 +k=0.99975 +x_0=500000 +y_0=500000 +ellps=krass +towgs84=28,-121,-77,0,0,0,0 +units=m +no_defs');
            ol.proj.proj4.register(proj4);

            console.log(coordinates);
            var marker = new ol.Feature(new ol.geom.Point(coordinates).transform('EPSG:4326', 'EPSG:31700'));
            console.log(marker);
            var zIndex = 1;

            marker.setStyle(new ol.style.Style({
              image: new ol.style.Icon(({
                anchor: [0.5, 36],
                anchorXUnits: "fraction",
                anchorYUnits: "pixels",
                opacity: 1,
                src: "img/airport.png",
                zIndex: zIndex
              })),
              zIndex: zIndex
            }));

            vectorSource.aircraft.addFeature(marker);
            map.addLayer(aircraft);
          }

Exemple coordinate:
[29.746857, 42.511047]

!!!!!UPDATE!!!!
I was able to display all the intercepted planes in ads-b:
The line had to be moved: map.addLayer(aircraft); when taking over json
$.post( "data.php",
                 {
                      get_data: 1
                 },
                 function( data ) {
                    var obj = JSON.parse(data);
                    var aircrafts = obj.acList;
                    //console.log(aircrafts);
                    for(var i = 0; i<aircrafts.length; i++ ){
                        if(aircrafts[i].Long || aircrafts[i].Lat){
                            var coordinates = [aircrafts[i].Long, aircrafts[i].Lat];
                            console.log(coordinates);
                            addMarker(coordinates);
                        }
                    }
                    map.addLayer(aircraft);
                    })


texts: How do I topologically clean "pointy bridges" using GRASS v.clean?

texts:
I have a shapefile with hundreds of pointy areas connected at the tip, forming a "pointy bridge". This creates a self-intersection and duplicate node. I'm not exactly sure what this type of error is called, so I do not know which v.clean parameters to use. How do I fix this using v.clean? I do not care about splitting the bridge vs. merging the bridge - just want to get rid of the error.

The last thing I tried was running v.clean with these parameters, and it did nothing:

break:0
 bpol:0
 rmdupl:0
 rmbridge:0
 rmdac:0


texts: Change layer datasource in qgis project with python

texts:
I have different qgis project with some layer and their datasource like ../my_path/my_file.shp.
I need to insert this shp file to a remote db. For do it i use shapefile from PGeo module and it's work fine.
But i need change layer datasource in qgs project with dbname= , host etc  where schema_name is a name of qgis project, table_name is a layer id, and i need to maintain all layer's property. I tried to use:

 uri = QgsDataSourceURI()
 uri.setConnection("localhost", "5432", "mydb", "postgres", "postgres")
 uri.setDataSource(str(schema_name),str(layer.id()),'the_geom')
 layer = QgsVectorLayer(uri.uri(), str(layer.id()), "postgres")
 QgsMapLayerRegistry.instance().addMapLayer(layer)
 layer.reload()
 print layer.source()
 if not layer.isValid():
    print "Layer is not valid"
 project.write(QFileInfo(str(dir)+'/'+project_name))

But i always got "Layer is not valid" and nothing is changed on my qgis project.
Is there another way to change layer datasource and save changes to the project? I see this solution but it doesn't work for me