texts: DSAS tool compatibility with ArcMap 10.4?

texts:
Using ArcGIS Desktop 10.4 I am attempting to quantify rate of shoreline change using the USGS's DSAS tool for arcmap. I am able to complete every step up until I run the statistics in which I get the error message:

C:\Users\mfazio\AppData\Roaming\USGS\DSAS\calc\DSASCore.calc.exe  -e -d "C:\Users\mfazio\AppData\Roaming\USGS\DSAS\data\20170912_101754.Tansect_S100.in.xml" -o "C:\Users\mfazio\AppData\Roaming\USGS\DSAS\data\20170912_101754.Tansect_S100.DSASCore.out.xml" -c "SCE,NSM,EPR,LRR,WLR,LMS"
[9/12/2017 10:18:40 AM] : ********EXCEPTION INFO BEGIN*********
[9/12/2017 10:18:40 AM] : System Exception:
[9/12/2017 10:18:40 AM] : Exception Source:
[9/12/2017 10:18:40 AM] : ESRI.ArcGIS.Geodatabase
[9/12/2017 10:18:40 AM] : Exception Target Site:
[9/12/2017 10:18:40 AM] : ESRI.ArcGIS.Geodatabase.ICursor Search(ESRI.ArcGIS.Geodatabase.IQueryFilter, Boolean)
[9/12/2017 10:18:40 AM] : Exception Message:
[9/12/2017 10:18:40 AM] : The server threw an exception. (Exception from HRESULT: 0x80010105 (RPC_E_SERVERFAULT))
[9/12/2017 10:18:40 AM] : System Exception:
System.Runtime.InteropServices.COMException (0x80010105): The server threw an exception. (Exception from HRESULT: 0x80010105 (RPC_E_SERVERFAULT))
   at ESRI.ArcGIS.Geodatabase.ITable.Search(IQueryFilter QueryFilter, Boolean Recycling)
   at DSAS.RateCalculation.readStatsFromFile(String outputTableName, String filename)

I'm wondering if this is an issue multiple people are having or just me.


texts: co-registering a vector file to a georeferenced raster image

texts:
i would like to co-register a vector file to a raster image by a simple translation in envi.

Is it possible to do this in envi?

I have tried using the georeferencing tool in arcgis but the vector layer cannot be selected for georeferencing!


texts: Creating PDF mapbook document that will contain map for each county in counties feature class using ArcPy?

texts:
Through a plethora of trials and research, I finally was able to get a script to work.

My original question encompassed me needing to take an .mxd and create a PDF mapbook document that will contain a map for each county in the counties FC. Then, I needed to change the title of each map to the name of the county being shown. Lastly, I need to scale the map a bit so it shows an area just a bit larger than the actual extent.

The following script works great. I hope this will help others in the future;

import arcpy, os
from arcpy.mapping import *

countiesPath = "C:\\Data\\COUNTIES.shp"
pdfPath = "C:\\Data\\COCounties.pdf"
mxdPath = "C:\\Data\\Counties.mxd"
mxd = MapDocument(mxdPath)

dataFrame = (mxdPath, "COUNTIES")[0]
listLayers = ListLayers(mxd)

finalPDF = PDFDocumentCreate(pdfPath)

df = ListDataFrames(mxd, "Layers")[0]

for layer in listLayers:
    if layer.name == "COUNTIES" :
        with arcpy.da.SearchCursor(layer, "COUNTY") as cursor:
            for row in cursor:
                countyName = row[0]
                whereClause = "COUNTY='"+ countyName + "'"
                arcpy.SelectLayerByAttribute_management(layer, "NEW_SELECTION", whereClause)
                df.zoomToSelectedFeatures()
                arcpy.SelectLayerByAttribute_management(layer,"CLEAR_SELECTION")
                df.scale *= 1.5
                titleItem = ListLayoutElements(mxd, "TEXT_ELEMENT", "title")[0]
                titleItem.text = "Map of: " + countyName + " COUNTY"
                countyPDF=r"C:\\Data"+countyName+".pdf"
                ExportToPDF(mxd,countyPDF,"PAGE_LAYOUT",resolution=600,df_export_width=1100,df_export_height =800)
                finalPDF.appendPages(countyPDF)
                os.remove(countyPDF)

finalPDF.saveAndClose()
del finalPDF
del mxd


texts:
I'm trying to use GDAL calc to convert values in a .tif file, I've done this before with GDAL calc but for converting to a single value. Now I'm trying to convert multiple values in this fashion:

0.0 - 0.2 -> 1

0.2 - 0.4 -> 2

0.4 - 0.6 -> 3

etc

I've tried this:

    calculation='((0.0<A<0.2)*1)+((0.2<A<0.4)*2)+((0.4<A<0.6)*3)+
    ((0.6<A<0.8)*4)+((0.8<A<1.0)*5)+((1.0<A<1.5)*6)+((1.5<A<2.0)*7)+
    ((2.0<A<2.5)*8)+((2.5<A<3.0)*9)+((3.0<A<9999)*10)'

    call([sys.executable, 'C:\Program Files (x86)\GDAL\gdal_calc.py', '-A',
    inraster, "--outfile=" + outraster, '--calc='+calculation,'--
    NoDataValue=-9999','--overwrite'])

This results in the following error:

    0 .. evaluation of calculation [thecalculation] failed

Any ideas? Maybe I need to split up the steps of converting but this could add significantly to the run time.


texts: How to handle/ avoid intersect errors when clipping the difference between two layers?

texts:
Using QGIS + Grass, I've added a shapefile of coastline and islands to my map.

I then buffer this shapefile to create a coastal zone of a known width. At this stage the buffer consists of the offshore (buffered) area + the area occupied by the land/ islands. As I am only interested in the offshore part, I then "clip" the land shapefile from the buffer.

I have tried this using both Vector>Geoprocessing Tools>Difference and GRASS v.overlay.xor. Either approach results in geometry errors (lots of them).

e.g.,

segments 0 and 24 of line 43 intersect at xx.xxx, yy.yyyy

This appears to be due to narrow connections between land areas (yellow in the picture below) resulting in intersecting nodes after the ocean (blue) is clipped.

My land data comes from an external source, and I don't have a lot of scope to simplify or adjust this.

What is the best way to handle these errors, or to avoid them in the first place?

I'm not a GIS expert, so there's a good chance that I'm approaching this problem from the wrong end. Would appreciate any tips or advice to avoid these errors.


texts: Performing math in Python parser of ArcGIS field calculator?

texts:
I never used python but now I need a code to do this calculation:

I'm doing this:

math.atan(math.fabs( !Y2! - !Y1!)/math.fabs( !X2! - !X1!) )*(180/math.pi)


texts: Equivalent of excel Alt+Enter in attribute table?

texts:
On Excel, you can do "carriage return", that is, jump a line, using Alt+Enter within a cell.

When I copy a cell from excel in my attribute table, this format is respected. I can find no way to write it directly within QGIS table. I tried Alt+Enter, Control+Enter, Shift+Enter... nothing works.

Has anybody got a clue ?


texts: Using Google Maps image as layer with Matplotlib Basemap?

texts:
I have produced a heat map of assaults in Sydney by suburb. However, my map looks a bit plain and I'd like to do something fancier - I'd like to use a Google Maps image as a background layer - I haven't figured out how to use a Google Maps image as a layer just yet and I'd like your help. I would like something that looks like this: http://www.ncomva.se/guide/Images/Layers.PNG As you can see, a Google Maps image is used a base layer. Then on top of the Google Maps image they put a see through heat map. Basically, I want to know is if this is possible in Python - I've tried to figure out how to do this but whatever I've tried just isn't working. I'm new to Python so this is difficult! If this is possible, can you please direct me to code that would be helpful?

I've also included some code below (but not the full code as it's too long) that gives you an idea how to use basemap.

import pandas as pd
import matplotlib
import matplotlib.pyplot as plt
import matplotlib.cm as cm
from matplotlib.colors import Normalize
from matplotlib.collections import PatchCollection
import matplotlib.artist
import matplotlib.colorbar
from mpl_toolkits.basemap import Basemap
from shapely.geometry import Point, Polygon, MultiPoint, MultiPolygon, GeometryCollection
from pysal.esda.mapclassify import Natural_Breaks as nb
import pysal.esda.mapclassify as mapclassify
from descartes import PolygonPatch
import fiona
from itertools import chain
import numpy as np
from shapely.geometry import asShape, mapping

#### http://pandas.pydata.org/pandas-docs/stable/generated/pandas.read_excel.html
data = pd.ExcelFile('data/sr16-14330 Jane Baseby Monthly counts of non-DV assault.xls')
data = data.parse('Suburb', skiprows=4, skip_footer=5, index_col=None, na_values=['NA'])

# load in the shapefile - use the file SSC_2016_AUST or "suburb_subset" for only a subset of suburbs
shp = fiona.open('data/ssc_2016_aust_shape/suburb_subset.shp')
bds = shp.bounds  # extract shapefile boundary
shp.close()

extra = 0.1  # define a variable called extra which we will use for padding the map when we display it (in this case I've selected a 10% pad)
#ll = (bds[0], bds[1])  # define the lower left hand boundary (longitude, latitude)
#ur = (bds[2], bds[3])  # define the upper right hand boundary (longitude, latitude)

NSEW = {'S': -33.892078720829694, 'E': 151.24185030131247, 'W': 151.18691866068747, 'N': -33.846468362714106}

ll = (NSEW["W"],NSEW["S"])
ur = (NSEW["E"],NSEW["N"])

coords = list(chain(ll, ur))  # concatenate the lower left and upper right into a variable called coordinates
w, h = coords[2] - coords[0], coords[3] - coords[1]  # define variables for the width and the height of the map

m = Basemap(
    # set projection to 'tmerc' which is apparently less distorting when close-in
    projection='tmerc',

    # set longitude as average of lower, upper longitude bounds
    lon_0=np.average([bds[0], bds[2]]),

    # set latitude as average of lower,upper latitude bounds
    lat_0=np.average([bds[1], bds[3]]),

    ellps='WGS84',

    # set the map boundaries. Note that we use the extra variable to provide a 10% buffer around the map
    llcrnrlon=coords[0],
    llcrnrlat=coords[1],
    urcrnrlon=coords[2],
    urcrnrlat=coords[3],

    # resolution of boundary database to use. Can be c (crude), l (low), i (intermediate), h (high), f (full) or None.
    resolution='f',

    # alpha = 0.5,

    # don't show the axis ticks automatically
    suppress_ticks=True
    #,epsg=4269
)

m.readshapefile('data/ssc_2016_aust_shape/test', 'assaults', color='none', zorder=2)


texts: ArcGIS Python Toolbox Parameter Filter.List ValueList Delete/Remove for Reconciling Versions

texts:
Does anyone know if and how you can remove a value from a parameter filter list?

In my case, I am using a python toolbox in ArcGIS 10.2.2 and want a drop down list of all the versions for the person reconciling to choose from. But I do not want to include dbo.DEFAULT and DS.Draft (our QAQC version, child of Draft).

param1.filter.type = "ValueList"

param1.filter.list = [v.name for v in arcpy.da.ListVersions(r"\\gis01\DatabaseConnectionFiles\DataStewardsOnly\GIS_Vector_DS.sde")]

I tried using -=("dbo.DEFAULT", "DS.Draft") and treating it like a list, but am just stuck on how to do this, since I read elsewhere that these act like tuples, and not like lists. And not sure if this must be done in one line of code or can be altered or if I need to use UpdateParams, but seems like that is for updating based on user input.

I am very new to python and scripting.


texts:
In my database of points with lat and long. After clipping the data to get only the ones in the specific state, my ID column that I had before disappeared. How can I clip my data and still keep the GUID that I had before?

I am clipping directly from Microsoft SQL server. On the server the GUID column is type uniqueIdentifier.


texts:
Goal: display map with few routes, which have colored segments. Also displays few waypoints with text + gives option to download few gpx routes and additional text.

I see that gpx format has possibility to define traces, routes, waypoints where you can define color for segments.

Question 1: what is the easiest solution to display such map? Is there some existing service (maybe something like http://umap.openstreetmap.fr) ?

Question 2: if it isn't and I have to do it with Leaflet, is there some plugin that helps me to display colored traces? I was searching through them but I did not see it directly mentioned in the documentation. I saw plugins displaying colors based on other things.


texts: Calculating LST from Landsat 8 in Google Earth Engine?

texts:
I'm trying to write a code for Land Surface Temperature (LST) from Landsat 8 images in google earth engine. I'm using the code of this question as a guide and I succeeded at intermediate steps but I couldn't solve the previous issues of the script nor make it work.

Here's the link to the code:
https://code.earthengine.google.com/2b2c4e1ece6ac2a201f8e854c34cd577

Any suggestions/ideas for getting the final LST product?
My final goal is to obtain land surface temperature (from landsat 8), so if you know another way to get there that would be helpful too.

Here's my code:

//cloud mask
function maskL8sr(image) {
  // Bits 3 and 5 are cloud shadow and cloud, respectively.
  var cloudShadowBitMask = (1 << 3);
  var cloudsBitMask = (1 << 5);
  // Get the pixel QA band.
  var qa = image.select('pixel_qa');
  // Both flags should be set to zero, indicating clear conditions.
  var mask = qa.bitwiseAnd(cloudShadowBitMask).eq(0)
                 .and(qa.bitwiseAnd(cloudsBitMask).eq(0));
  return image.updateMask(mask);
}

//vis params
var vizParams = {
  bands: ['B5', 'B6', 'B4'],
  min: 0,
  max: 4000,
  gamma: [1, 0.9, 1.1]
};
var vizParams2 = {
  bands: ['B4', 'B3', 'B2'],
  min: 0,
  max: 3000,
  gamma: 1.4,
};

//load the collection:
{
var col = ee.ImageCollection('LANDSAT/LC08/C01/T1_SR')
.map(maskL8sr)
.filterDate('2018-01-01','2018-12-31')
.filterBounds(geometry);
}
print(col, 'coleccion');

//median
{
var image = col.median();
print(image, 'image');
Map.addLayer(image, vizParams2);
}

// NDVI:
{
var ndvi = image.normalizedDifference(['B5', 'B4']).rename('NDVI');
var ndviParams = {min: -1, max: 1, palette: ['blue', 'white', 'green']};
print(ndvi,'ndvi');
Map.addLayer(ndvi, ndviParams, 'ndvi');
}

//

//select thermal band 10(with brightness tempereature), no BT calculation
 var thermal= image.select('B10').multiply(1000);
 Map.addLayer(thermal, imageVisParam, 'thermal');

// find the min and max of NDVI
{
var min = ee.Number(ndvi.reduceRegion({
   reducer: ee.Reducer.min(),
   geometry: geometry,
   scale: 30,
   maxPixels: 1e9
   }).values().get(0));
print(min, 'min');
var max = ee.Number(ndvi.reduceRegion({
    reducer: ee.Reducer.max(),
   geometry: geometry,
   scale: 30,
   maxPixels: 1e9
   }).values().get(0));
print(max, 'max')
}

//fractional vegetation
{
var fv = ndvi.subtract(min).divide(max.subtract(min)).rename('FV');
print(fv, 'fv');
Map.addLayer(fv);
}

/////////////

  //Emissivity

  var a= ee.Number(0.004);
  var b= ee.Number(0.986);
  var EM=fv.multiply(a).add(b).rename('EMM');
  Map.addLayer(EM, imageVisParam2,'EMM');

  //LST c,d,f, p1, p2, p3 are assigned variables to write equaton easily
  var c= ee.Number(1);
  var d= ee.Number(0.00115);
  var f= ee.Number(1.4388);

var p1= ee.Number(thermal.multiply(d).divide(f));
var p2= ee.Number(Math.log(EM));
var p3= ee.Number((p1.multiply(p2)).add(c));

  var LST= (thermal.divide(p3)).rename('LST');

  var LSTimage = ee.Image(LST)

  Map.addLayer(LSTimage, {min: 0, max: 350, palette: ['FF0000',
  '00FF00']},'LST');


texts:
I have this file http://www.meteofrance.com/mf3-rpc-portlet/js/datas/zones_AVDEPT73.json where the "zone" field is put in an area HTML element (so does "shape").

I have a little script that takes the zone and transform it into GeoJSON. Even if the GeoJSON is not at the right place here is the script:

from typing import Dict

def parse_geom(zone: str) -> Dict:
    geom = {'type': 'Polygon', 'coordinates': None}
    coords = zone.split(',')
    # we want the area to be projected in EPSG:4326. But lat need to be bound
    # in -180;+180 and long -90;+90.
    coords_int = [float(x)/100 for x in coords]
    geom['coordinates'] = [[coords_int[x+1], coords_int[x]] for x in range(0, len(coords_int)-1, 2)]
    geom['coordinates'] = [geom['coordinates']]
    return geom

Here is the result (in geojson.io)

Here is the entire script to downloads the file:

https://gist.github.com/RemiDesgrange/fd29aadb9493b1d7a88b7eb246a93425

I tried to rotate it via QGIS. But I don't understand what I see. It seems to me that this is a fairly simple math problem...


texts: OpenLayers-v5: How to get WFS features before vector layer is visible

texts:
I expected a WFS call to be made when a Vector Layer was added to the map, but it doesn't happen until the layer is marked as visible.

Here is how I'm currently creating my vector layers. The visible param is initially marked as false via layer.checked

function createVectorLayer (layer, zIndex) {
  const source = new VectorSource({
    format: new GeoJSON(),
    url: function (extent) {
      return `${process.env.GEOSERVER_URL}/wfs?service=WFS&version=${process.env.WFS_VERSION}&request=GetFeature&typename=${layer.name}&outputFormat=application/json`
    }
  })
  return new VectorLayer({
    source,
    style: VectorStyle,
    opacity: layer.opacity,
    visible: layer.checked,
    zIndex
  })
}

I then add them to the map via the standard method:
map.addLayer(vectorLayer)

The problem is that until I update the layer's visibility via
vectorLayer.setVisible(true) the wfs call made to the url defined in VectorSource won't be made and therefore the layer has no features which is causing problems.

Again, I assumed would be that the wfs call would be made when the Layer was added to the map and not on the condition of its visibility.

I haven't been able to find anything in the docs that explains this, so I could use some help/information.

How can I get a vector features (via wfs call) before the layer is marked visible?

  app-specific information

I'm creating a feature search that sets the map's center to that feature's coordinates
but, if the layer is not 'visible' yet (the user just opened the map and everything is at default settings), the layer has no features to search from

I can then set the layer as visible on search, but then the wfs call is made in the background and I have no way to listen to it's return. I could listen to a layer's change event but that's irregardless of the main question of:
how (if possible) can I get a layer's features via wfs before the layer is visible


texts: Binarization of water body extraction in Google Earth Engine?

texts:
I'm a new learner of Google Earth Engine.

What I am doing right now is to extract water bodies from Landsat 8 using such water index as NDWI(normalized difference water index) in GEE.

Is there any available function in the binarization of the water bodies to obtain the optimal threshold (such as the OTSU)?

I found it difficult to perform the binarization after calculating the NDWI.


texts:
I have a roads table and I want to break the lines in the intersections like the v.clean tool from GRASS (I'm talking about how the break tool in the parameters works).

How can I do the same process in PostGIS?


texts: Batch project a folder of shapefiles to one projection with string appended to each output

texts:
I'm trying to create a relatively basic batch reprojection tool with Python by projecting all shapefiles in a folder to NAD83 UTM Zone 11N.

I want them to go to a specified folder that I have defined as a variable. I want the outputs to maintain the same name but have "_projected" appended to the final outputs. I also have defined a variable that is the file path to a vector dataset whose spatial reference will be used in the reprojection. I want to do all this with a For loop.

At this point, the outputs go to this file path below and the wrong name, which is not what I want:

"C:\\WCGIS"

How to make it go this file path below and to have the string, "_projected" appended to each newly projected shapefile?:

"C:\\WCGIS\\Project2Output"

Right now I'm trying this part of the code I'm posting just below to have the outputs go to the folder called Project2Output and it keeps on giving me names like Project2OutputCityBoundaries, Project2OutputFerries, Project2OutputStateRoads, etc, and to "C:\WCGIS".

        outPath = outLocation + sFile + "_projected"

        arcpy.Project_management(sFile, outPath, outProjection)

I will eventually be turning the variable paths targetData and out projection as arcpy.GetParameterAsText variables for my Python tool but just want to see this work in my IDE first. Here is all my code below:

import arcpy
import os

arcpy.env.overwriteOutput = True
targetData = "C:\\WCGIS\\Project2Data" # Shapefiles that need to be projected
outLocation = "C:\\WCGIS\\Project2Output" # Where the projected shapefiles will be stored
outProjection = "C:\\WCGIS\\Project2Data\\CountyLines.shp" # Shapefiles with correct target projection: NAD 1983 UTM Zone 10

# Get a list of all feature classes in the Project2Data folder
arcpy.env.workspace = targetData

try:
    for sFile in arcpy.ListFeatureClasses():
        dsc = arcpy.Describe(sFile)
        if dsc.spatialReference.Name == "NAD_1983_UTM_Zone_10N":
            print "Skip this shapefile because it already has the correct projection:" + sFile
        else:
            outPath = outLocation + sFile + "_projected"
            arcpy.Project_management(sFile, outPath, outProjection)
            print arcpy.GetMessages()

except arcpy.ExecuteError:
    print(arcpy.GetMessages(2))

except Exception as ex:
    print (ex.args[0])


texts: My fusion table map with search not working

texts:
My maps are not working. It shows the search box but not the map:

// ]]></script>

&nbsp;

 Postcode
  (with space e.g. HP11 2QW) <input id="search-string_0" type="text"> <input
   onclick="changeMap_0()" type="button" value="Search">


texts: Combining popup info from overlapping polygons in Mapbox GL JS

texts:
I've tried several solutions, but I'm having trouble populating a popup with information from overlapping polygons. When there are polygons that overlap, I want to drill down and pull the names from all polygons at the click point. Any help on this would be awesome.

 // When a click event occurs on a feature in the states layer, open a popup at the
location of the click, with description HTML from its properties.

map.on('click', 'testlayer', function (e) {
    new mapboxgl.Popup()
        .setLngLat(e.lngLat)
        .setHTML(e.features[0].properties.name)
        .addTo(map);
});

The jsfiddle here


texts: Selecting more than one feature from raster attribute table using ArcGIS Desktop?

texts:
How do I select more than one features from raster attribute table?

For example:

"lc_2001" = 21, 22, 23, 24 AND "lc_2006" = 11, 31, 41, 51


texts: How to display dynamic text from a joined table attribute

texts:
I'm using ArcMap 10.4 with dynamic pages setup. Everything I've found online says that if you join a table to the index layer you should be able to add the data from one of the table's fields via dynamic text (example: How to change the dynamic information in data driven pages). Nothing is working.

Does anyone have any suggestions?

I can't just add a field to the index since it's used in a bunch of other files and I can't get access to unlock it. The join works fine, I can see the attributes in the attribute table of the index layer, but no joy on getting them into dynamic text in anyway.

currently using: <dyn type="page" property="attribute" field="CONCEPT_NAME" domainlookup="true" />

This just yields a text box saying [EMPTY]


texts: Stretchingt geometry (shapes) up to desired point in any direction using FME?

texts: How to stretch geometry -line, up to desired point in any direction in FME?


texts: Is there a QGIS 3 setting to completely turn off plugins?

texts:
I would like to turn off all plugins without uninstalling them.  I usually do this while testing to rule out an issue being caused by a plugin.  I am looking for the QGIS equivalent to the ArcMap Administrator setting BlockAddIns being set to 4:

0—Load all add-ins regardless of whether or not they have digital signatures.
1—Only load and use add-ins that are digitally signed by a trusted certificate authority.
2—Only load and use add-ins that have been published by ESRI.
3—Only load add-ins from the administrator folders (described below) and all add-ins published by ESRI.
4—Do not load or execute any add-ins in this application.


texts:
I'm trying to install landsat-util in an Ubuntu AWS instance. First, I created a Python 2.7 environment called landsat_py2 to avoid conflicts and the installation process seems to be fine. Nonetheless, when I run the help: landsat -h, the following appear:

(landsat_py2) imendieta@wef-oceans:~$ landsat -h
Traceback (most recent call last):
  File "/home/imendieta/.conda/envs/landsat_py2/bin/landsat", line 7, in <module>
    from landsat.landsat import __main__
  File "/home/imendieta/.conda/envs/landsat_py2/lib/python2.7/site-packages/landsat/landsat.py", line 29, in <module>
    from .image import Simple, PanSharpen, FileDoesNotExist
  File "/home/imendieta/.conda/envs/landsat_py2/lib/python2.7/site-packages/landsat/image.py", line 16, in <module>
    import rasterio
  File "/home/imendieta/.conda/envs/landsat_py2/lib/python2.7/site-packages/rasterio-0.32.0-py2.7-linux-x86_64.egg/rasterio/__init__.py", line 8, in <module>
    from rasterio._base import eval_window, window_shape, window_index
ImportError: libgdal.so.1: cannot open shared object file: No such file or directory

[UPDATE: After running $ conda install -c conda-forge gdal rasterio jpeg=8d matplotlib=1.5.1 the problem resolves and Landsat-util runs flawlessly.]


texts: Will NMEA support the GNxx string in future?

texts:
At the moment there is no problem when using GPS only devices, NMEA string with information about speed, DOP value etc starts with GP.. prefix. But with Glonass included NMEA string has different prefix GN.. and information about DOP and speed are not recognized in Qgis. Is it planned to update this also?
Position is ok, in both cases.
regards,
radek


texts: Since I have little coding skills, I'm looking for an open source bus route planner that consumes GTFS data (That I have to create myself), All I could find was OpenRouteService, and although it has some pretty nice features, the bus routing informations are very lacking (no bus numbers, stations name, transit...), Looking for something similar to Google Transit.


texts:
ArcMap 10.2.

I want to export map in draft mode. Without layers. Only few rectangles displaying main part of map, legend, scale bars, north arrows, margins and guide lines.

E.g. something like this:

Is it possible? I've tried to export by Export Map but I got normal map.


texts: How to do kernel density estimation with a physical barrier/ boundary in QGIS?

texts:
Does anyone know of a programme that will allow kernelling that takes a boundary into account in its calculation rather than simply masking areas that are impossible?

So far I have found:
GME add on for Arc10 (Hawthorne Beyer) - I get an error everytime I specify a boundary shapefile.  I've tried lots of different shapefile types and boundary complexities etc.  It works fine when I don't specify a boundary.

AdehabitatHR package in R (Calange 2011) - this works well but the boundary you specify has to be very simple - line segments 3x kernel bandwidth in length, and not too tortuous.  For my data this is a big oversimplification.

So I'm wondering if any other software can do it, GRASS or QGIS for example.

Thanks


texts:
My goal is to display rotated labels for roads on my openlayers based webmap.
The problem is actually solved (with a custom SVG solution), but i have some problems with the styling of the labels.
As I using the custom SVG js, 'labelOutlineColor' and 'labelOutlineWidth' (which I use as styleMap parameters) stop doing anything whith the labels. Do I have to implement those parameters into my custom SGV js? And if is that so, how can I do it?

Here is the CustomSGV.js I use:

OpenLayers.Renderer.CustomSVG = OpenLayers.Class(OpenLayers.Renderer.SVG, {

drawText: function(featureId, style, location) {

    var resolution = this.getResolution();

    var layer = this.map.getLayer(this.container.id);

    var feature = layer.getFeatureById(featureId);

    location = (feature.attributes.centroid ? feature.attributes.centroid : location);

    var x = (location.x / resolution + this.left);

    var y = (location.y / resolution - this.top);

    var label = this.nodeFactory(featureId + this.LABEL_ID_SUFFIX, "text");

    var tspan = this.nodeFactory(featureId + this.LABEL_ID_SUFFIX + "_tspan", "tspan");

    label.setAttributeNS(null, "x", x);

    label.setAttributeNS(null, "y", -y);

    label.setAttributeNS(null, "pointer-events", "none");

    if (style.angle || style.angle == 0) {

        var rotate = 'rotate(' + style.angle + ',' + x + "," + -y + ')';

        label.setAttributeNS(null, "transform", rotate);

    }

    if (style.fontColor) {

        label.setAttributeNS(null, "fill", style.fontColor);

    }

    if (style.fontFamily) {

        label.setAttributeNS(null, "font-family", style.fontFamily);

    }

    if (style.fontSize) {

        label.setAttributeNS(null, "font-size", style.fontSize);

    }

    if (style.fontWeight) {

        label.setAttributeNS(null, "font-weight", style.fontWeight);

    }

    var align = style.labelAlign || "cm";

    label.setAttributeNS(null, "text-anchor",

        OpenLayers.Renderer.SVG.LABEL_ALIGN[align[0]] || "middle");

    if (this.isGecko) {

        label.setAttributeNS(null, "dominant-baseline",

            OpenLayers.Renderer.SVG.LABEL_ALIGN[align[1]] || "central");

    } else {

        tspan.setAttributeNS(null, "baseline-shift",

            OpenLayers.Renderer.SVG.LABEL_VSHIFT[align[1]] || "-35%");

    }

    tspan.textContent = style.label;

    if(!label.parentNode) {

        label.appendChild(tspan);

        this.textRoot.appendChild(label);

    }

},

CLASS_NAME: "OpenLayers.Control.CustomSVG"
});


texts:
I need to hide a layer with a function, use this script:

function hideLayer(layerName) {
        layerName.setVisibility(false);
}

I have this error: 'map' is null or undefined

my code page:

var map;
var lon =9.10;
var lat =45.46;
var zoom =10;
function init(){
    epsg4326 = new OpenLayers.Projection("EPSG:4326")
    map = new OpenLayers.Map("map", {displayProjection:epsg4326,
    controls: [
        new OpenLayers.Control.Navigation(),
        new OpenLayers.Control.MousePosition(),
        new OpenLayers.Control.LayerSwitcher()
        ]
    });
var newLayer = new OpenLayers.Layer.OSM("Mappa locale", "Tiles/${z}/${x}/${y}.png", {numZoomLevels: 19, alpha: false, isBaseLayer: true});
 map.addLayer(newLayer);
 var layer = new OpenLayers.Layer.Vector("Mezzi disponibili", {
   strategies: [new OpenLayers.Strategy.BBOX({resFactor: 1.1})],
   protocol: new OpenLayers.Protocol.HTTP({
       url: "./sinottico_mezziliberi.txt",
       format: new OpenLayers.Format.Text()
        })
    });
map.addLayer(layer);


texts: On which TextBox event should I call ICompletionNotify.SetComplete?

texts:
I have a IToolControl derived control which contains a c# TextBox.
In the IToolControl::OnSetFocus function documentation, it states that you need to call ICompletionNotify::SetComplete() when it's finished to let the application know that the control should lose focus.

On another IToolControl ComboBox control, I called the SetComplete() function on the DropDownClosed event.

On which TextBox event should I call ICompletionNotify.SetComplete()?


texts: Getting unique values from multiple fields as matched using PyQGIS

texts:
There are Code and Color fields in my layer like this: (Sample 6 rows)
Code  Color
----  -----
A     Red
B     Blue
A     Red
A     Red
B     Blue
C     Green
...   ...

All As are Red, all Bs are Blue and so on.
I can get the unique values for each field:
layer = iface.activeLayer()

code_index = layer.fields().indexFromName("Code")
code_unique = layer.uniqueValues(code_index)

color_index = layer.fields().indexFromName("Color")
color_unique = layer.uniqueValues(color_index)

print(code_unique)
print(color_unique)

# OUTPUT
# {'A', 'C', 'B'}
# {'Blue', 'Red', 'Green'}

As you can see, the order of two outputs doesn't match. (When you try, you might get an ordered result, but it depends and is not reliable)
How can I get unique value pairs like [('A', 'Red'), ('B', 'Blue'), ('C', 'Green')] using PyQGIS?
I couldn't find groupby-like method. It can be done by means of groupby(["Code", "Color"]).first() in GeoPandas.


texts: How to cut lines in a shapefile following a GPX-source?

texts:
I have a OSM-Shapefile containing all roads of a given area. I need an overlay cut out of this source-Shapefile that contains only some routes that I want to highlight on my map. I have some GPX-Tracks that show those routes, but they aren't precise enough to use them directly as an overlay.
How can I use those GPX-Tracks to select the corresponding roads of the original shapefile, cut them out and copy them in a new shapefile? Is it possible to set a certain "buffer area" for selecting lines also if they don't exactly overlay with the GPX-file?
I'm using QGIS normally, and possibly I'm searching for a opensource solution.
Thank you!


texts:
I am facing this problem 3th day and I really don't know where is the problem, so please help me.

Here is my REQUEST: http://geoportal.cuzk.cz/WMS_SM5V_PUB/service.svc/get?SERVICE=WMS&VERSION=1.3.0&REQUEST=GetMap&LAYERS=SM_STP_HRBITOV&STYLES=&CRS=EPSG%3A102067&BBOX=-895000,-865000,-623000,-620000&WIDTH=1600&HEIGHT=610&FORMAT=image/png&

Here is Capabilities file: http://geoportal.cuzk.cz/WMS_SM5V_PUB/WMService.aspx?service=WMS&request=GetCapabilities&

Can somebody please explain me why blank map is returned?


texts: Removing triangles that fall into the concave but exterior part of polygons in ArcGIS Desktop?

texts: I have created a TIN for all the polygons in a block and then converted it to a feature layer in order to be able to delete the triangles that fall in the polygons. Now I want to remove triangles that fall into the concave but exterior part of polygons. But I don't know how to do it in ArcGIS. This condition makes triangles with all three vertices located on the same object to be removed from consideration.  The screenshot of part of my data is as the first picture. And I want to remove exterior white triangle as shown in the second picture. And My final purpose is to form a skeleton graph for each street block, like the third picture. In other words, I want to implement Thiessen polygon for polygons instead of points but I don't know how to do it. If I would be able to implement Thiessen polygon for polygon objects, it would be used instead.


texts: What does it mean to QI and how is it done?

texts:
I'm new to Esri's Arc* interfaces, and trying to figure out how to register an existing table using the ArcObjects API. The following question answers this to some extent:

Is there any way to write a snipet of code with function of "SDELayer -o register" command

  Maybe the IClassSchemaEdit.RegisterAsObjectClass method will help (see the GeoDatabase OMD). You should be able to QI it from an ITable, opened in the usual fashion (e.g. IFeatureWorkspace.OpenTable(in Name: String) since you've already built your GDB).

But there is some terminology used that I don't understand. Specifically, what does it mean to QI something, and how is it done?

This other question's answer defines QI, somewhat vaguely, as Query Interface (hopping) but goes no further.

Navigating ArcObjects Object Model Diagrams (OMDs)?

  It really helps if you know UML. In the class that I had in college the professor called it QI (Query Interface) hopping. The Getting To Know ArcObjects book has a nice break down of how to work with the diagrams.

I'm using Java - is this a .net thing?


texts:
Is it possible to convert a list of pixel locations (txt file) into lat long?

What I have is:

(left_x:  137   top_y:  146)
(left_x:  175   top_y:  236)
(left_x:  182   top_y:  153)
(left_x:  185   top_y:  166)

What I would like is to convert them to lat long so I can generate a points layer.

Currently using QGIS, and am new to Python.


texts: Does arcpy have any tools for transforming field values in a fieldmapping?

texts:
I built my fieldmappings and, for the most part, everything is working. I recall reading somewhere that I could also use fieldmappings to transform field values during the import process. Unfortunately, I can't find the link where I read this. For example, in field ABC, replace value "1" with "TRUE" or if field ABC is null, replace it with DOG.

Can I do achieve this using tools within arcpy? Or, do I need to build this from scratch?


texts: Exporting Geometry Collection as Shapefile on Earth Engine

texts:
I have a FeatureCollection that I cannot export as a shapefile on Earth Engine —  I get the following error: Error: Shapefiles cannot contain multiple geometry types; found 'LineString', 'Polygon'.
This is the export code:
Export.table.toDrive({
  collection: table,
  description:'vectorsToDriveExample',
  fileFormat: 'SHP'
});

And upon using table.geometry(), I get the following output:
GeometryCollection
type: GeometryCollection
geometries: List (10 elements)
0: LineString, 3 vertices
1: LineString, 3 vertices
2: LineString, 3 vertices
3: LineString, 2 vertices
4: Polygon, 11 vertices
5: Polygon, 1484 vertices
6: Polygon, 7113 vertices
7: Polygon, 4 vertices
8: Polygon, 4 vertices
9: Polygon, 10205 vertices

I suppose it would be necessary to convert the LineString to Polygon, but I haven't found a way to do that.


texts: Assigning vector data attributes from other adjacent data

texts:
I am supposed to do a presentation on how to do this, and I'm not exactly sure what it is asking.

The question is: "how can you assign vector data attributes from other adjacent data?"

I am in a geology field class and am learning GIS and mapping my field sites as well. This is supposed to be something I am able to learn through the help menus and/or forums, but the question itself is not clear to me and they won't answer me. Anyone have an idea what this would most likely be asking?


texts: Changing DEM with 10 meter cell size into 30 meter cell size in ArcMap?

texts:
I have a 10 meter cell size Oregon DEM that I would like to match my Northern California DEM thats at a 30 meter cell size.

How would I go about changing the Oregon DEM cell size from a 10 meters to 30 meters using ArcGIS Desktop 10.5?


texts: How can I make the topology checker plugin ignore very tiny errors?

texts:
I'm a QGIS2.0 user.
Topology checkers is a very valuable for me as a data creator.
I use this plugin for checking data quality.

But, this plugin finds too much errors include very, very tiny errors.
I have tried to set a "tolerance" for the rule but it does not work well.

Please tell me how to use this function.


texts: Geoserver - raster layer title in getfeatureinfo template

texts:
Relatively new to GeoServer, moving from MapServer and I'm trying to generate a GetFeatureInfo template for a raster layer.

I'm trying to pull in the layer title via "${feature.type.title}", which seems to work fine for vector layers, however, when attempting this with the (multiband) raster layer, it seems to return "Layer: " followed by the layer name.

Is there any way to access the title for this raster layer?


texts: How to filter the graph on which I want to find the shortest path

texts:
I have a large graph that contains edges of Europe. I want to calculate the shortest path between two points on this graph with pgRouting that have a function pgr_trsp for shortest path calculations.

This pgr_trsp needs a sql query as a text that filters the edges of the graph for a specified shortest path calculation. Since the calculation is very slow if I give the whole Europe graph to the pgr_trsp, I want to give only the specified territory that is used by the shortest path.

But I have to be careful, because pgr_trsp will calculate the shortest path only on the filtered graph. If this filtered graph does not contain edges that are necessary for finding the real shortest path, then I will get a false shortest path.

So the trade-off is to filter the graph to contain as few edges as it is possible, but do not miss any edge that is necessary for calculating the real shortest path.

My edges are classified into 6 classes regarding to its importance:

0 : roads of international importance (most important highways)
1 : roads of national importance (highways and most important roads)
2 : roads of regional importance (making larger settlements and parts of settlements accessible)
3 : local roads of high importance (making minor settlements accessible)
4 : local roads
5 : local roads of minor importance (unpaved roads, roads in 30km/h zones, back roads etc.)
6 : restricted roads (roads that are not accessible for normal traffic)

I made a quick demonstration about the roads class(0,1,2,3). On the picture below you can see that how the road class(0,1,2,3) graph looks:

It is a continuous graph of course, but here I show only smaller circles due to the performance of the visualization tool.

All of these colorful edges belong to class(0,1,2,3). You can see on the very left part of the picture a green circle with several blue markers. This is a big city where the edges (class(0,1,2,3)) are very dense: average distance between edges about 1.5km. They are much rarer in the country where the average distance between two edges about 20km.

It means that if I put a point anywhere on the map, I can find a more important road within 1.5km in the city, and 20km in the country. So if I find the closest edge for the point, I need to involve every roads in 1.5km diameter circle around the closest edge in the city and 20km diameter circle in the country into shortest path calculation. I make this step for the second point as well, and after that I involve the class(0,1,2,3) roads between the two points around the line that connects them.

Now it is clear that if I use a 1.5 km diameter in the country it is very possible that I will not find the link to the roads class(0,1,2,3) that are between the two points around the connecting line. If I use 20km diameter for finding the link between the surroundings of the points and the higher classified roads, then I will involve too many edges in the city and the calculation will be slow.

I can complicate the query if I try to find whether a point is in a city or not, and belongs to this use a longer or shorter diameter. But now it is difficult to decide whether a point is in the city or out of it. And it seems to me that this solution is not the best way.

Does anybody know how the Titans (big companies) do this? What is the best practice? You can suggest me a book, or article as well. But if you know the best solution, you are welcome of course! :-)


texts: How can I divide a single postgis linestring with SRID 4326 into segments (preferably an array of linestrings) with no segment longer than a given distance? I've tried st_segmentize(), but it spits out another single linstring.


texts:
A client has asked me to provide a dataset in seamless table format. Could anybody describe to me what a seamless table is?
I understand that it is a format used by mapinfo, but can I create a seamless table in some other way, by FME for example?
Thanks.


texts: QGIS print composer to image export don't match

texts:
In QGIS I have set up a map in print composer but when I export to image or pdf the saved image does not match what is in the print composer.

Details - I am using an openstreetmap as my background map and I have a group of around 360 points witch I have uploaded from a cvs file. On the print composer the points are on one side of a river and after export they are on the other (6-7 km difference)

Any suggestions would be helpful. It seem to be working fine earlier this week.

Map as seen in print composer:

Map when exported:


texts:
I am attempting to create Supervised Land Cover classification maps using the classes (Mangrove, Non-Mangrove, Barren, Urban, and Water) for four period of interests (1988, 2000, 2017 and 2020) using ArcMap 10.8 Maximum Likelihood Classifier (MLC). I am using Landsat images acquired from USGS Earth Explorer. The objective is to evaluate the extent of changes in each of the classes for each period.
I have been successful at classifying the images for 1988 and 2000 using the 1988 signature file. However, when attempting to run the MLC via a Supervised Classifier on the remaining images I get the following error messages pictured in the screenshot below.


texts: Understanding nested list of asPolygon in PyQGIS

texts:
In Getting list of QgsPoint from polygon Layer in PyQGIS, a user explains how to get the coordinates of a polygon and gives the following code:

layer = iface.activeLayer()
feature = layer.getFeatures().next()
polygon = feature.geometry().asPolygon()
n = len(polygon[0])
for i in range(n):
     print polygon[0][i]

What is polygon[0] in this case and what would be polygon [5] for example?

In the for i in range(n) part it seems that i represents the vertex number of the polygon and its coordinates. But I don't get the first polygon[0].


texts:
I am looking for a way to calculate sequential numbers in points feature along line.

I've got polygon 1 feature intersected with other layer of polygons. Then I extract points from the intersection, and get points1 layer. After that I extract the points from the vertices of polygon (1) (points2 layer). Also I create points along line by distance (points3). Then I merege all layers together (points1 + points2 + points3). Then I want to give them a sequential numbers going up along the border of polygon 1 by order.

My goal is to add points along the line that will be in ascending order, in favor of an orientation map for contractors in the field.

Is there a way to calculate it by field calculator?

I am using Arcmap 10.6.


texts: How to add a halo to a dynamic label of an attribute in Atlas?

texts:
How can I apply this for a dynamic label of an attribute in atlas? I'm using QGIS 2.18.10

How to add a halo to a text in the composer?


texts: Convert Java object to JSON with jackson

texts:
I am trying to read the data from a shapefile and once read, I want to convert it to a JSON file.

To read the data I'm using the geoTools library, and to write to JSON, I use the geojson-jackson library (https://github.com/opendatalab-de/geojson-jackson).

I have managed to correctly read the information data extracted from the shapefile, but I cannot get the JSON file to be well generated.

This is part of my code:

Read Data from .shp

public void esShapefileParcel(){

// ShapeFile Directory
                String archiveshpParcel = directory+"\\"+constant.getTmpParcel()+"\\"+constant.getTmpParcel()+".shp";

                Parcel parcel = new Parcel();

            try {

                File shapeFile = new File (archiveshpParcel);
                FileDataStore store = FileDataStoreFinder.getDataStore(shapeFile);
                SimpleFeatureSource featureSource = store.getFeatureSource();
                SimpleFeatureCollection collection = featureSource.getFeatures();
                SimpleFeature fParcel;

                SimpleFeatureIterator itrParcel = collection.features();

                while(itrParcel.hasNext()){

                   parcels.clear();

                   fParcel = itrParcel.next();
                   parcel.setExternalReference((String)fParcel.getAttribute("REFCAT"));

                   parcels.add(parcel);

                   ObjectMapper Obj = new ObjectMapper();
                   String jsonStr = Obj.writeValueAsString(parcels.toString());
                   System.out.println(jsonStr);

                }

        } catch (IOException e) {
            e.printTrace();
        }

            }
        }

And, my Parcel object is:

public class Parcel{

    protected String type;
    private String externalReference;

    Parcel(){};

    Parcel(String type, String externalReference) {
        this.type = "parcel";
        this.externalReference = externalReference;
        }

    public String getType() {
        return type;
    }

    public String getExternalReference() {
        return externalReference;
    }

    public void setExternalReference(String externalReference) {
        this.externalReference = externalReference;
    }

    @Override
    public String toString() {

             return "  \"type\" : \"" + type + "\""+
               ", \"externalReference\" : '" + externalReference + '\'' ;

    }

But this, only returns a invalid Json:

  "[  \"type\" : \"null\", \"externalReference\" : '3545013XG2734N']" "[
  \"type\" : \"null\", \"externalReference\" : '8401550XG1670S']" "[
  \"type\" : \"null\", \"externalReference\" : '8401522XG1680A']" "[
  \"type\" : \"null\", \"externalReference\" : '9293903XG1699C']" "[
  \"type\" : \"null\", \"externalReference\" : '9293907XG1699C']"


texts: Distinguishing lines that intersect from lines that touch?

texts:
How do I distinguish between these cases in ArcGIS 10?

Case 1: Both endpoints of a line touch another line
Case 2: Both endpoints dangle over the lines it intersects

I am looking at the Trim Line function but that's not what I want (destructive).

The real world use case is to distinguish between service roads connecting both roadways of a freeway, and other cases of roads intersecting with freeways.


texts: Google Maps Ground Overlay: How to calculate bounds?

texts:
I'm working with the Google Maps API for the first time and I'm trying to append a simple image of a Building on top of my current map view.

This seems to be very straight forward when I look at the GMapsAPI here …
https://developers.google.com/maps/documentation/javascript/overlays#GroundOverlays

However, I have no idea how to calculate the bounds and the exact position of my overlay.

As you can see in the documentation the position of the ground overlay is done via LatLng

var imageBounds = new google.maps.LatLngBounds(
    new google.maps.LatLng(40.716216,-74.213393),
    new google.maps.LatLng(40.765641,-74.139235));

In my case I have no idea how to calculate those bounds. What are those values exactly? Why are there two seperate LatLng positions needed?

I have a simple .png of a building that I want to layer on top of an exact position. I get this position by entering the destination in Google Maps and reading the Lat and Lng values from the embed-code.

Is there a better way of doing this? Any ideas, tipps or tricks on how I could do this?

Or is it possible to calculate those bounds when I know the center of the image-position and the dimesions of my image file? (My image is 50x50px and a .png)


texts: Extraction of line segments in front of the street

texts:
I have a layer of lines that represents a layer of plots. I want to eliminate all the inner segments to have only the segments that go out into the street. Thus, I want to determine the sum of the segments that lead to the road.
I started like this:

We removed the double segments with the function "Delete duplicate geometries"
I used dissolve on the group field to form a single line for a single plot.
I turned the lines into polygons.
I dissolved again on the resulting polygons to remain only with the contour of the plots, without the inner lines.

I used "Extract by location", using the touch parameter to extract the face to the street (only the lines facing the street)
<img src="https://i.stack.imgur.com/qHUNy.png" alt="
" />

This is the end result. Do you have any idea why these lines are broken and do not return the outline of my plot? Do you have any idea if there is a difference between layers? Or if not, you know some other way to extract these lines, but extract it from the original file, to keep the tabular data and segment characteristics.

After buffering the layer, it looks like this:
I set the tolerance of 0.00000001. However, I think it works even smaller.


texts: Basic script for Raster Files - PyQGIS 3

texts:
I am trying to make a script, to run at the toolbox,  that receives a raster file and has an identical output to it, following the idea from here.

I couldn't find a complete script that works unchanged so I can understand how to manipulate rasters with PyQGIS 3.
The most complete I found was the answer from root676 here,
but I couldn't execute the script.

So I would like guidance on how to set up a simple script for input and output of raster files, without using native QGIS or GDAL tools, to not pre-define the output file name.

EDIT:

When I run the script mentioned, I get:

  The following layers were not correctly generated.C:/Users/cs279069/AppData/Local/Temp/processing_26353d6a350741f99b38886692ba0306/35605f7cacd647be810cbfc73ee98186/OUTPUT_RASTER_A.tifC:/Users/cs279069/AppData/Local/Temp/processing_26353d6a350741f99b38886692ba0306/d89846f8fbb8400da1e7891ff94270f0/OUTPUT_RASTER_B.tifYou can check the 'Log Messages Panel' in QGIS main window to find more information about the execution of the algorithm.

The Log Messages Panel is empty.
Probably because there is no calculation on it.

Clarifying the question:

The input raster can be any raster of any size, because I'll use this in several different scripts, just to rename the output file.
When I create some script, I use the native tools of QGIS, and that generaty an output with the tool predefined name. Like "Reprojected", "Clipped (mask)".

So basically I need something like this: OUTPUT_RASTER = INPUT_RASTER * 1

So the output will be identical from the input and will have the display name that was defined in the def initAlgorithm section.

More info:

I´m trying to run this code.

def initAlgorithm(self, config=None):
    self.addParameter(QgsProcessingParameterRasterLayer(
        self.INPUT_RASTER,
        self.tr("Input Raster"), None, False))

    self.addParameter(QgsProcessingParameterRasterDestination(
        self.OUTPUT_RASTER,
        self.tr("Output Raster"),
        None, False))

def processAlgorithm(self, parameters, context, feedback):
    input_raster = self.parameterAsRasterLayer(parameters, self.INPUT_RASTER, context)
    output_raster = self.parameterAsOutputLayer(parameters, self.OUTPUT_RASTER, context)

    output_raster = input_raster

    results = {}
    results[self.OUTPUT_RASTER] = output_raster
    return results

That´s the error given:

  The following layers were not correctly generated.C:/Users/cs279069/AppData/Local/Temp/processing_779dd8aa0c284112b878829e6b56c348/4bacded9efb04c6c8c9502e6572f5c1d/OUTPUT_RASTER.tifYou can check the 'Log Messages Panel' in QGIS main window to find more information about the execution of the algorithm.


texts: How can I make use of the $scale function in QGIS data defined properties?

texts:
In QGIS 2.0 I can use with success the $scale function in the labeling expressions.

When I tried to use it for symbolizing in the data defined properties it has a value of 0.

Is there a way to make it perform like it does in labeling (getting the scale)?

As I saw in the expression string builder help that

$scale function
Returns the current scale of the map canvas.
Note: This function is only available in some contexts and will be 0 otherwise.

is there a source where I can read more about this e.g. in what contexts it is available?
Thanks, Tudor


texts: How can I rotate a map layer in QGIS so that center of map always have north upwards?

texts:
I'm working in the polar regions and make maps in different scales. The convention is to show e.g. the Antarctic continent with the prime meridian upwards on the map and polar stereographic projection. In this way, south is toward the center of the map and east is clockwise (and opposite for the Arctic region).

However, when working on higher resolution, I want to have north facing upwards on the map, both on the canvas and in the composer, even if the layers and my project are polar stereographic projection.

The rotate function for the canvas is really handy, and I know that I can also rotate a map in the composer and use an expression. This is how I do it today, but it's complicated as I have to change the rotation every time I adjust the extent of the map, due to the converging meridians near the poles.

I could imagine a couple of ways to simplify my work flow:

A dream scenario would be if the center of the view was always be rotated so that true-north faces upwards. Just as working in a Mercator projection.
Or, how should I write an expression in the print composer to rotate the map to 360 - lon_c where lon_c is the longitude of the center of the map extent. My data is in UTM, but can I access the coordinates without writing a complex formula to convert meters to degrees?
A last option would be to pass the canvas rotation to map rotation in the composer. Can it be scripted?

How can I achieve any of those options? I don't want to change the projection, just rotate the view so that north at center of the map is facing upwards in print composer.


texts: r.slope.aspect grass gis datum format on dx and dy analysis

texts:
I'm performing terrain analysis due to photovoltaic installation on QGIS with r.slope.aspect plugin. I'm interested in evaluate differently N-S slope and E-W slope based on the PV system I would install.
Processing dx and dy tool I don't understand in which format the output raster is formatted.
Choosing degrees or percent the results are the same. For example the output can be -1.17 NS (percent) what does it mean -1.17? 117%.


texts: gdalwarp cuts off top and right edge of globe from gfs data

texts:
I am using data from this GFS forecast: https://nomads.ncdc.noaa.gov/data/gfs4/201905/20190502/gfs_4_20190502_1800_000.grb2

Specifically the bands 102 and 103 containing wind data.

I first combine the u and v components (band 102 and 103) using:

gdal_calc.py --overwrite -A gfs_4_20190502_1800_000.grb2 --A_band 102 -B gfs_4_20190502_1800_000.grb2 --B_band 103 --outfile=temp.tiff --calc="sqrt(A*A + B*B)"

Then I'm reprojecting the data to an orthographic projection using:

gdalwarp -overwrite -t_srs "+proj=ortho +lat_0=14 +lon_0=80 +x_0=0 +y_0=0 +R=6371000 +units=m +no_defs" -of GTiff -ts 2000 0 "temp.tiff" "output.tiff"

But the end result for that is a globe that's chopped off at the top and right by a few degrees. like here:

Can't figure out if this is something in the data or something in my commands/projection that's wrong?


texts: Selecting intersecting elements that have same value in another attribute and merging them into one using ArcPy

texts:
I have a polyline shapefile "rivers" containing three columns:

[Start,  End,  Meas]
[1,     2,     3]
[2,     3,     3]
[2,     4,     2]
[4,     5,     3]
[5,     6,     3]

What I want to do is inside a for loop test all elements in the shapefile (large number)
If Start of one element = End of another element & Meas is the same the join the 2 elements.

The result would be:

[Start,  End,  Meas]

[1,     3,     3]
[2,     4,     2]
[4,     6,     3]

Any ideas how to do this (merge 5 elements into 3 elements) in arcpy?

I mostly want these 2 elements to be merged spatially. Don't really care if the Start and end columns are updated by the script correctly..

A small sample of my attribute table:

What I am trying to do is a for loop.. I think that I have to use cursors, but I am new at this.. My for loop is like this but cant assign to literal (Of course maybe I am completely wrong about this): `>>>

... for i in "Channels2020":
...     temp3=["NODE_A"]
...     temp3a=["ORDER"]
...     temp3b=["SEGMENT_ID"]
...     for k in "Channels2020":
...         temp4a["NODE_B"]
...         temp4a=["ORDER"]
...         temp4b=["SEGMENT_ID"]
...         if temp3b!=temp4b & temp3==temp4 & temp3a==temp4a:
...             ["SEGMENT_ID"]=temp3b


texts: Change True/False boolean value to "Occupied" label

texts:
Change the label for a boolean: Instead of "True" and "False" I want the label to display "Occupied."

I have tried many iterations of this expression.  Nothing has worked.

Case when "Field Name" = "TRUE" THEN "Occupied"
Else Null
End


texts: How can I write the output of a ModelBuilder process to my attribute table?

texts:
Using ArcGIS 10.2, I've put together a Model Builder model that counts the number of wind turbines in a buffer zone:

When I'm done, I can see how many turbines have been selected at the bottom of the attribute table panel like this:

How can I use Model Builder to capture this output (e.g., the number of turbines "1040") and save it to an existing (empty) field in my attribute table?


texts: Importing 14 bit raw TIFF Thermal data from FLIR Zenmuse XT DJI drone camera

texts:
I am trying to import Thermal data from a FLIR Zenmuse xt radiometric lens thermal camera.

Flir does not include a conversion with the camera and does not respond to my emails regarding the conversion.

On the Pix4D website they include the conversion for the camera which is

"0.04 * [Thermal Image] - 273.15"

This equation works for the most part but some images can be as much as 30 degrees F off.

Currently working in QGIS to apply the Pix4D conversion using raster calculator.

I have attached the TIFF as well as the .RJPEG for viewing in Flir tools.

Please let me know if you have any questions or if there is anything else I need to add.

The general temperature range I am looking for is 32.8-106.6 F

Thermal images can be found here
https://drive.google.com/drive/folders/0BzQ1355jgmF9Y0kycmdlcWZiLVE?usp=sharing


texts: Filter values to include features within a specific max/min range

texts:
I have a table of values that have a MIN and MAX integer value, i.e.

I would like to select all values that range from Age_MIN >= 105 to Age_MAX <= 110.  There are features in the table that would fall into that range, but I'm not sure how to query them.

When I do a definition query with Age_MIN >= 105 to Age_MAX <= 110 I get this result:

As a visual representation of what I would like to do, please see below image:

Here, the features (horizontal bars) in blue would be selected, and the yellow bars would be unselected as they don't fall within the range of 100 to 110.

Any ideas?


texts: Improving this script - GPX reading and processing with Python OGR

texts:
I am trying to plot every GPX trace I have in a folder, by reading it with Python-OGR, clipping it by a bounding box, and plotting

Curently I have a working script, but it seems very uggly, and I suspect there should be a more straightforward way to achieve the same results.

Specific questions in comments:

#coding:utf-8

import glob
import ogr
import json
import numpy
import matplotlib.pyplot as plt

# 1) Is there a more direct or idiomatic way to create this rectangle?
ring = ogr.Geometry(ogr.wkbLinearRing)
ring.AddPoint(-51.22, -30.027)
ring.AddPoint(-51.22, -30.047)
ring.AddPoint(-51.18, -30.047)
ring.AddPoint(-51.18, -30.027)
ring.AddPoint(-51.22, -30.027)

poly = ogr.Geometry(ogr.wkbPolygon)
poly.AddGeometry(ring)

for filename in sorted(glob.glob('../../../04 GPS/REGISTROS/*/*.gpx'), reverse=True)[:1]:

    # 2) Are all those "if not X: continue" statements the right way to check?
    datasource = ogr.Open(filename)
    if not datasource:
        continue

    layer = datasource.GetLayerByName('tracks')
    if not layer:
        continue

    feature = layer.GetNextFeature()  # Shouldn't I get *every* feature?
    if not feature:
        continue

    geometry = feature.GetGeometryRef().Intersection(poly)
    if not geometry:
        continue

    # 3) Exporting was the only way I found
    # to get numeric values. Is there any other?
    j = json.loads(geometry.ExportToJson())

    # 4) one of returned objects have only "geometry" key...
    if not 'coordinates' in j:
        continue

    for segment in j['coordinates']:
        seg = numpy.array(segment, ndmin=2)
        plt.plot(seg[:,0], seg[:,1], lw=0.5, alpha=0.1)

plt.axis('equal')
plt.show()


texts: Memory limit eror counting every pixel in an image with Google Earth Engine

texts:
I am trying to count the number of pixels in my large AOI (roughly the size of the state of Iowa) for a Landsat derived image.  To get a rough pixel count, I could divide the area of my ROI by 900 (30x30m scale) but I prefer GEE's ee.Reducer.count(), which include partial pixels in the ROI if the pixel's centroid is there.
Printing the pixel count to the console returns Dictionary (Error) Computation timed out.
var pixel_count = ndvi_mosaic.reduceRegion({
  reducer: 'count',
  geometry: aoi,
  scale: 30,
  maxPixels: 1e20,
});
print(pixel_count)

Exporting pixel_count to a CSV also gives an error Error: User memory limit exceeded.
Export.table.toDrive({
  collection: ee.FeatureCollection([
    ee.Feature(null, pixel_count)
  ]),
  description: 'pixel_count',
  fileFormat: 'CSV'
});

I've also tried setting tileScale up to 16 and bestEffort: true, with no luck.  Any suggestions that don't involve having to reduce the size of my image?


texts: Generating a raster layer with intersecting line lengths as cell values

texts:
I have a polyline vector layer and a template raster dataset. For this polyline vector layer I would like to generate a new raster dataset based on the template raster dataset, but in which each cell contains a value corresponding to the length of the line segment that crosses that cell. All other cells should contain a value Nodata. Heres a stupid example I made in Paint:

I am using python based tools including rasterio, geopandas, shapely and rasterstats. I would like to avoid using lower level libraries such as gdal and ogr if I can help it, but I will resort to them if it makes the difference between life and death.

What I am currently doing is the following:

Extract template raster dataset by mask using polyline dataset
reassign all values in extracted dataset a unique value
polygonize raster dataset based on following function

def polygonize_raster(dataset):

    # Read the dataset's valid data mask as a ndarray. Dataset is a rasterio read object open for reading
    mask = dataset.dataset_mask()

    array = dataset.read(1)
    generator = rasterio.features.shapes(source=array, mask=mask, transform=dataset.transform)
    # Extract feature shapes and values from the array
    geom_list = []
    for geom, value in generator:
        # Print GeoJSON shapes to stdout
        geom = shapely.geometry.shape(geom)
        geom_list.append(geom)
    return geom_list

generate a new list with line lengths intersecting each polygon based on the following functions:

def calc_length_in_cell(polygon_geom, linestring_geom):

    length = linestring_geom.intersection(polygon_geom).length
    return length

def calc_length_in_cell_list(polygon_geom_list, multilinestring):

    #loops through all polygon geometries, calculates cell length inside and returns list with all lengths
    lengths = []
    for i,polygon_geom in enumerate(polygon_geom_list):
        lengths.append(calc_length_in_cell(polygon_geom, multilinestring))
    return lengths

add the list of lengths to a geodataframe containing the polygonized raster (little boxes corresponding to each template raster cell).
burn the above geodataframe into a new raster with the line lengths as burn value

I have 3 issues with this workflow:

Smaller tests with square rasters containing about 9 cells with unique values seem to polygonize just fine, and each cell gets its own polygon. However for large datasets this stops being the case, as evinced by the image below
 despite the fact that adjacent cells have different values
the calc_length_in_cell_list function above is a killer. It really blows up for large datasets since it has to perform an intersection for every polygon. Heavy stuff.
the polygonized raster is technically invalid geometry as can be determined by this question.

I would like to either 1) come up with a different workflow or 2) tweak the current one such that this can be accomplished for large datasets with ease, and that the resulting raster with intersecting polyline lengths is correct.


texts:
I would like to enable a Login form for Geoserver. This would means that, when a client makes a request using Leaflet or OL3, with an index file etc. The client would be prompted to log in before the layer would load. Any Advice, I have consulted the security of Geoserver but its like i don't understand a thing.
Geoserver 2.10


texts: Convert dijit.TitlePane to dijit.form.ToggleButton

texts:
I have an app that I'm building using the javascript API and have created 3 titlepanes in the header that open up to reveal information (basemap gallery, legend, layers). It functions exactly how I want it to, now I would like to make it "look pretty" by giving each pane its own image, where when you click on the image, it opens the pane in the same way.

I have tried overriding the CSS however, that gives all the titlepanes the same image. It looks as though I'll need to convert the titlepanes to a togglebutton. Reading through some of the samples it looks like I'll need to add additional code beyond the html change. Is there a best way to go about this edit?

EDIT:
Here is an example of the html that creates the titlepane the "basemapGallery" div is the information inside the pane that I would like to make a drop down or any other type of button where I can control the image.


texts: Get specific parameter in CRS information with Python gdal

texts:
I have a tif image containing its CRS. I would like to extract the parameter "central_meridian". I have the following piece of code:

from osgeo import gdal,osr
ds=gdal.Open(r'hrs0000bf5c_07_if177l_trr3_CAT_scale_trial_p.img.tif')
prj=ds.GetProjection()
print prj

srs=osr.SpatialReference(wkt=prj)
print srs.GetAttrValue('PARAMETER',0)
print srs.GetAttrValue('PARAMETER',1)

The output I get is the following:

    PROJCS["Mars_Npolar_Stereo",GEOGCS["GCS_Unknown",DATUM["Unknown",SPHEROID["S_Unknown",3396190,169.8944472236118]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433]],PROJECTION["Polar_Stereographic"],PARAMETER["latitude_of_origin",90],PARAMETER["central_meridian",160.0712280273438],PARAMETER["scale_factor",1],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]]]

latitude_of_origin
90

I would like to fetch the second parameter "central_meridian" but i don't know how to do it? Any ideas?

##EDIT

If trying to open a list of files given by a list like:

from osgeo import gdal,osr
from gdalconst import GA_ReadOnly

f = open('list.txt')
for line in iter(f):
    print line
    ds=gdal.Open(line,GA_ReadOnly)
    prj=ds.GetProjection()
    srs=osr.SpatialReference(wkt=prj)
    print line + " " + srs.GetProjParm( osr.SRS_PP_CENTRAL_MERIDIAN)
f.close()

I get the following error:
ERROR 4: `frt0000703c_07_if187l_trr3_CAT_scale_trial_p.img.tif
' does not exist in the file system,
and is not recognised as a supported dataset name.

I tried different ways to do it and went through documentation but no luck. Any ideas why it does not recognizes it as a filename?


texts: Where can I find the shapefiles for Canadian postal codes?


texts:
I ran this query successfully, but the column highwayid isn't updated as expected:

update unique_latlongs set highwayid =
  (SELECT h.osmid FROM rn.highways h WHERE h.geom &&
  ST_Expand(ST_Transform(st_setsrid(st_makepoint(lon, lat),4326), 3414), 50)
  ORDER BY ST_Distance(h.geom,
  ST_Transform(st_setsrid(st_makepoint(lon, lat),4326), 3414)) limit 1);

What do you think could be wrong in this query? In unique_latlongs, lat and lon are columns.

EDIT: Perhaps I should explain what I was trying to do: I have lat and long values in the unique_latlongs table. I also have the column highwayid, which I want to update to the id the the closest highway to each respective lat long in unique_latlongs. The highways are store in the rn.highways table, with the h.geom being the geometry value of the highway.


texts:
I want to run Jupyter Lab (Notebook) with ArcPy. I have ArcMap 10.7 with Python Version 2.7.15, Numpy=1.9.3, matplotlib=1.5.2.

I followed the instructions from the following link:

https://gisday.wordpress.com/2016/07/18/setting-up-anaconda-pysal-with-arcgis-python-environment/
and I have downloaded Anaconda 2 32-bit.

I can only run the Anaconda Powershell Prompt as an Administrator, but this option does not exist for the Anaconda prompt (I do not know whether this makes a difference...).

The problem I have occurs when running the command:

  conda create -n arc107 python=2.7.15 numpy=1.9.3 matplotlib=1.5.2
  pyparsing xlrd xlwt pandas scipy ipython ipython-notebook
  ipython-qtconsole

I get the following error:
PackagesNotFoundError: The following packages are not available from current channels:
then it lists: matplotlib=1.5.2, ipython-notebook, ipython-qtconsole

I deinstalled Anaconda and reinstalled it again, but it does not work. I also did not get the thing running with Anaconda 3, 32-bit.


texts: How to get in point geometry to pinpoint accurately where geometries have mistakes?

texts:
Basically I have some line and some polygon tables which have invalid some geometries and apparently PostgreSQL returns errors when trying to return reports which rely on intersecting these features
To help people fix these mistakes I'm trying to come up with a view which points out the mistakes.
I could easily just do:
    SELECT id,name,geometry
    FROM   polygontable
    WHERE  ST_IsValid(geometry) is not true

The problem is that people will have to check the entire polygon when it's probably a single node that is wrong.
How to actually return the point location of the mistake?


texts: Check/Uncheck Trim Polygon option in Network Analyst

texts: Can anyone tell me when to check/uncheck Trim Polygon option in Service Area Layer Properties dialog box while doing Network Analyst problems? By default Trim Polygon values comes as 100 m. Shall we need to change it according to our problem?


texts: Copying field names, aliases and datatype to Excel through clipboard using ArcGIS Pro?

texts:
I am making an inventory of the field names, aliases and their data types to get an overview of my entire geodatabase. I want to therefore copy them to excel.

It works on ArcMap.

However, I cannot copy aliases easily.

The format I want is in ArcGIS Pro's design table but I cannot cop it.

Any suggestions?


texts: Automating citation of data source in QGIS print composer?

texts:
I have 20 thematic maps with different sources like CITY PLANNING OFFICE, ENVIRONMENT OFFICE, AGRICULTURE OFFICE etc.

Since I have many maps, I want to automate the data source citation of my thematic maps by simply having a code or process. I want that each map layout, the data source will change based on the map that is being viewed. This will lessen the burden because I am not having to manually type the data source for each map layout.

I saw someone discussed about variables and it got my hopes high but it didn't work for me because the layer variables cannot be used in text label in map layout. Only the Global and project variables were working.


texts:
I created a model in QGIS 3.6 with three inputs parameters (two vector layer and one vector field).
When I run the model, the inputs are not sorted in a correct way: the vector field is requested before the vector layer.

For this simple example, it doesn't matter, however for more complex models with severals inputs, I fear the user could feel not confortable.

Is it possible to choose the order of the inputs parameters ?

Here are two screenshot.
"Conduite" = pipe in french
"regard" = manhole in french

I want the user to select firstly "regard" layer, then "ID_regard" field within "regard" layer.


texts: Convert Polygon to Multipolygon with GeoPandas

texts:
I'm trying to create a vectors grid from a single geometry, from PostGIS's table to another PostGIS's table; both tables are in the same DB. The input geometry is MultiPolygon and I've previously created with SQL the output table:
CREATE TABLE IF NOT EXISTS comuni_campania_copy (
    id SERIAL PRIMARY KEY,
    comune VARCHAR(255) NOT NULL,
    comune_id INTEGER NOT NULL
);
SELECT AddGeometryColumn('comuni_campania_copy', 'geom', 32633, 'MULTIPOLYGON', 2);

I'm be able to create with GeoPandas a vectors grid from a shapefile and save it as another shapefile. I've modified my script to use this functionality with PostGIS.
from shapely.geometry import MultiPolygon
from sqlalchemy import *
import geopandas as gpd
import numpy as np
import psycopg2

# Connection data
con_params = create_engine('postgresql://postgres:1234@localhost:5432/postgres', echo=False)

# SQL parameters
sql_params_in = "SELECT id as comune_id, geom, comune FROM comuni_campania WHERE id = 3"

def vector_grid_from_db(sql_params: str, con_params: psycopg2.extensions.connection, dx: int, dy: int) -> None:
    # Read vector
    vector_in = gpd.GeoDataFrame.from_postgis(
        sql=sql_params,
        con=con_params
    )

    # Read bounding box
    bbox = vector_in.bounds

    # Generate columns and rows
    cols = list(np.arange(bbox.minx[0], bbox.maxx[0], dx))
    rows = list(np.arange(bbox.miny[0] + dy, bbox.maxy[0] + dy, dy))

    # Create vector grid
    polygons = []
    for x in cols:
        for y in rows:
            polygons.append(MultiPolygon([(x, y), (x + dx, y), (x + dx, y - dy), (x, y - dy)]))

        grid = gpd.GeoDataFrame(geometry=polygons, crs=vector_in.crs)

    epsq_code = int(str(vector_in.crs)[5:])

    grid.to_postgis(
        'comuni_campania_copy',
        con_params,
        if_exists="append",
        index=False,
        dtype={
            "geometry": Geometry(
                "MULTIPOLYGON",
                srid=epsq_code
            )
        }
    )

vector_grid_from_db(sql_params_in, con_params, 300, 300)

Using my script I see this error:

Traceback (most recent call last):   File "create_grid_to_db.py", line
71, in
vector_grid_from_db(sql_params_in, con_params, 300, 300)   File "create_grid_to_db.py", line 45, in vector_grid_from_db
polygons.append(MultiPolygon([(x, y), (x + dx, y), (x + dx, y - dy), (x, y - dy)]))   File
"/home/maxdragonheart/DEV_FOLDER/MIO/GISPython/devenv/lib/python3.8/site-packages/shapely/geometry/multipolygon.py",
line 62, in init
self._geom, self._ndim = geos_multipolygon_from_polygons(polygons)   File
"/home/maxdragonheart/DEV_FOLDER/MIO/GISPython/devenv/lib/python3.8/site-packages/shapely/geometry/multipolygon.py",
line 181, in geos_multipolygon_from_polygons
N = len(exemplar[0][0]) IndexError: invalid index to scalar variable.


texts: Classifying LiDAR ground points using laspy?

texts:
I have a classified LiDAR point cloud (.las) with misclassifications. The walls of steep surfaces that should be ground (class 2) have been misclassified as vegetation (classes 3, 4, 5). I am attempting to reclassify these erroneous vegetation points as ground.

My approach is to:

Find all single-return vegetation points within my data set.
Iterate through each of the above points and get a list of all its neighbors within an 80 meter radius.
If one of those neighbors is a ground point with a higher elevation, then I reclassify the original point from vegetation (3, 4, or 5) to ground (2).

I have tried to implement this workflow using a combination of laspy, numpy, and scipy.spatial:

import laspy
import numpy as np
from scipy.spatial import cKDTree
import matplotlib as plt

#read in my original, misclassified LAS file. Create a point_records array
in_file = laspy.file.File(r"misclassified.las", mode = "r")

point_records = in_file.points

#find all points that and single-return vegetation points (class 3, 4, or 5)

#note: num_returns is not included in the point_records above. So I am
#creating new arrays here.
class_arr = in_file.raw_classification
return_arr = in_file.num_returns

single_veg_pts = np.where(np.logical_and(np.logical_or(class_arr == 3,
                                                       class_arr == 4,
                                                       class_arr == 5), return_arr == 1))

#get the filtered points' records
single_veg_pt_records = in_file.points[single_veg_pts]

#create an x, y, z array of all points to build a 3D KDTree
point_records_xyz = np.array((point_records['point']['X'],
point_records['point']['Y'], point_records['point']['Z'])).transpose()

ctree = cKDTree(point_records_xyz)

#loop through my single vegetation point records to compare their X, Y, Z
#against my 3-D KDTree matrix. When a nearest neighbor is found...
for idx in single_veg_pt_records:
    neighbors = ctree.query_ball_point(np.array((idx['point']['X'],
    idx['point']['Y'],idx['point']['Z'])).transpose(), 8200)

for neighbor in neighbors:
   #if the neighbor point is class 2 and above idx, change the class of idx
   #in point_records.

The final for loop is where I am hitting a wall.

How do I compare idx to all of its neighbors. How can I  modify the original set of points (point_records) using the index value of idx?


texts:
The size of mdf file is much bigger than its gdb file.

As all the content of mdf file is exported (or imported) to a gdb, the size drops dramatically. The size of mdf is 8GB while its corresponding gdb is 200MB!!!

How can I “clean” the mdf file so that it gets as small as the gdb?


texts:
I have a Uniform Grid of 1KMx1KM squares as a shapefile with population data in each grid in a specific column as an integer.  I would like to extrude each polygon based on the population and export the whole dataset using a python library.  ArcScene can do this, but ArcPy can't hook into that program.

Can someone recommend a python package to do this?

It doesn't have to be dependent on ArcGIS Desktop.

I was looking at matplotlib, but I can't figure out how to create 3D bar graphs.


texts: Adding Map Key to map or console in Google Earth Engine?

texts:
Is there a way to print (or plot) a map key in the console (or on the map) in google earth engine ? What would be an approach?

var img= ee.Image('COPERNICUS/S2....................');
var ndvi = im.normalizedDifference(['B8', 'B4'])
var ndvi_brg = {min:-1, max: 1, palette: ['blue','#8B0000','green']};
Map.addLayer(ndvi.clip(rect), ndvi_brg, 'ndvi',true);

I'm ploting an NDVI in a specific area and would like to print the map key for the color palette 'ndvi_brg' as a visual clue
(Like something in the following picture:)


texts: Layer properties dialog in my QGIS plugin

texts:
I am developing a plugin which should display some information like geometry of the layer, feature count etc. Now those are currently displayed in QMessageBox.information. I want it to be displayed similar to the QGIS Layer properties Tab.

What changes should be done in the code?


texts: Iterate through TOC elements in ArcGis Server 9.3

texts: I'm developing website with a custom tool in ArcGis Server 9.3 using Visual Studio 2008 and the ADF components, the thing is that when I click on my custom tool it calls a class, this class needs to check each layer available to see if it is enabled and visible (this properties changes when you zoom out or zoom in), does anyone knows how to accomplish this? Thanks.


texts: Converting .vrt (with many tiles) to a single .tif using gdal_translate creates an empty raster

texts:
I can't get gdal_translate to convert a .vrt file with ~1400 100x100 pixel tiles (5 bands) to a single .tif with 5 bands. The results are always an empty raster brick – each band is empty...
Here is the command:
gdal_translate -of GTiff biomass.vrt biomass.tif

Here is the input (.vrt):

And here is the result (.tif):

UPDATE
This may be a versioning issue. The command works on my local machine with GDAL 3.1.4, but fails on my ubuntu aws box with GDAL 2.2.3.


texts: return the text entered in a QTableWidget

texts:
I am creating a plugin in Qgis and I need to enter values to a QTableWidget through the interface, but I cannot find a function that returns them to me.
I do not do it with a LineEdit because I need the number of spaces to fill out equal to those of a listWidget.
try to do it with QTableWidget itemAt but the return is a QTableWidgetItem.
the code is the following:

def Input_shape (self):
    filefolder= QFileDialog.getExistingDirectory(self.dlg, "", "/desktop")
    filename = glob.glob(filefolder+"/"+"\\*.shp")
    self.dlg.comboBox_2.clear()
    self.dlg.listWidget.clear()
    self.dlg.comboBox_2.addItems([filefolder])
    self.dlg.listWidget.addItems(filename)
    self.dlg.tableWidget.clear()
    self.dlg.tableWidget.setRowCount(len(filename))
    self.dlg.tableWidget.setColumnCount(2)

def run(self):
    if self.first_start == True:
        self.first_start = False
    self.dlg = SubsueloDialog()
    self.dlg.boton_entrada.clicked.connect(self.Input_shape)
    filefolder= self.dlg.comboBox_2.currentText()
    filenames = glob.glob(filefolder+"/"+"\\*.shp")
    limites= self.dlg.tableWidget.item(1,1)

The interface looks like this:


texts: Has anybody else experienced having issues opening QGIS or .gqz project files on macs?

texts:
I have an old Macbook pro with QGIS 3.10 installed. I have not been able to open any QGIS projects I have saved because of this error -> There is no application set to open the document “test project save.qgz”.
I was thinking this was most likely a problem of where I had QGIS installed on my computer. I checked and it is under Macintosh HD -> Applications. I have another Applications folder if I go to my user profile from Mac HD which does not seem to have QGIS or R or several other of my apps I regularly use. Instead that folder has many games I play. I am the only user of my computer and so there is only my profile and the default guest one which is empty.
I dont think that QGIS being in the base applications folder for all users should be problematic. That being said, I have no idea why I am getting this error when trying to open my saved QGIS projects. I have used QGIS quite a bit for school and personal interest. I never had a problem when opening saved QGIS projects on school computers other then them being so slow. I have been getting by with just saving layers which is annoying because it doesn't save the properties of layers either.
Please Help if you can, either with being able to save properties of layers on saved layers or with troubleshooting my issue of opening saved projects. Thanks!


texts: How to add text on top of symbols for vector layers in QGIS

texts: I read the QGIS 3.2 Changelog and saw a new feature to add text on top of symbols for vector layers. Does anyone knows how to do this?


texts: Connecting to PostgreSQL using QGIS using a VirtualMachine?

texts:
I installed Postgresql in Windows server 2012 Essential on Virtual Machine.

Everything works ok , but when  i try to use QGIS Desktop which is on desktop installed to access to Postgresql-data , the connection fails or (Error), known that i use the IP of virtual machine in Host. see image


texts: Reproject Equirectangular .cub file to Orthographic Project - using GDAL

texts:
I’m trying to plot an image of the Moon, with a good version of the topography (or even a nice picture). I then want to overlay the Apollo landing sites. This locations should be pretty accurate - although the final image need only be less than 8 cm x 8 cm. I’d like to use python Basemap and I’m reasonably familiar with its use. I probably want the final image to be an Orthographic projection.

If anyone has simple suggestions - such as a ready made netCDF in the right projection, please let me know!

I think the LOLA topography has what I want, but I’m struggling to transform file formats and the projection.

I have downloaded the file  Lunar_LRO_LOLA_Global_LDEM_118m_Mar2014.cub   from
https://astrogeology.usgs.gov/search/details/Moon/LRO/LOLA/Lunar_LRO_LOLA_Global_LDEM_118m_Mar2014/cub

The output of the command

  gdalinfo Lunar_LRO_LOLA_Global_LDEM_118m_Mar2014.cub  -mm is shown
  below. This is an Equirectangular projection in meters. The
  information is in Band 1, is Int16 and does not have a colour
  specified.

So I need to transform it to a netCDF file (this seems to work).
I also need to reproject it to an 'ortho' projection, and I need to downsample it. I am using gdal, which is new to me. Neither command is working properly - it seems that the resampling is messing with the range of the output (it should be from  approximately -5000 to +1000, but the range is much greater after resampling, and the colour becomes specified as grayscale.

downsample:

gdalwarp -tr 1024 1024 -r bilinear Lunar_LRO_LOLA_Global_LDEM_118m_Mar2014.cub a1_resampled.cub

reproject:

gdalwarp a1_resampled.cub b1.cub -t_srs "+proj=ortho +a=1737400. +units=m"

make a netCDF file:

gdal_translate -of netCDF b1.cub a_trans.nc

(gdal)bash-3.2$ gdalinfo Lunar_LRO_LOLA_Global_LDEM_118m_Mar2014.cub  -mm

Driver: ISIS3/USGS Astrogeology ISIS cube (Version 3)
Files: Lunar_LRO_LOLA_Global_LDEM_118m_Mar2014.cub
Size is 92160, 46080
Coordinate System is:
PROJCS["SimpleCylindrical Moon",
    GEOGCS["GCS_Moon",
        DATUM["D_Moon",
            SPHEROID["Moon",1737400,0]],
        PRIMEM["Reference_Meridian",0],
        UNIT["degree",0.0174532925199433]],
    PROJECTION["Equirectangular"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",0],
    PARAMETER["standard_parallel_1",0],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0]]
Origin = (-5458203.076608000323176,2729101.538304000161588)
Pixel Size = (118.450587600000006,-118.450587600000006)
Corner Coordinates:
Upper Left  (-5458203.077, 2729101.538) (180d 0' 0.00"E, 90d 0' 0.00"N)
Lower Left  (-5458203.077,-2729101.538) (180d 0' 0.00"E, 90d 0' 0.00"S)
Upper Right ( 5458203.077, 2729101.538) (180d 0' 0.00"W, 90d 0' 0.00"N)
Lower Right ( 5458203.077,-2729101.538) (180d 0' 0.00"W, 90d 0' 0.00"S)
Center      (   0.0000000,   0.0000000) (  0d 0' 0.01"E,  0d 0' 0.01"N)
Band 1 Block=1024x1024 Type=Int16, ColorInterp=Undefined


texts: When converting xls to DBF table, can I change the length of text fields?

texts:
I used the code:

arcpy.ExcelToTable_conversion(XLSfile, DBFfile)

The data comes across looking fine, but all text fields in the dbf have length of 254.  Is there a way, either in the excel file before I export it, or in the dbf after or during the export, to specify the length of that text field to something other than 254?

Things I've figured out so far:
ExcelToTable is pretty good at figuring out which field is which, I don't have any errors where the data came in as the wrong type.
The length of the actual data in excel doesn't seem to effect this (unless it is longer than 254 characters, which is another can of worms).


texts: Using multiple colours in the same symbol in arcmap 10.4

texts:
I am working on a map of the Rainy-Lake of the Woods watershed on the U.S.-Canada border.
I have created a shapefile of water quality monitoring stations in the basin using data from USGS and Environment and Climate Change Canada (ECCC). I have attached a screenshot of the Excel spreadsheet I'm basing it off of below.

Is there a way to split a symbol by multiple colours in arcmap to not only show whether the station is run by USGS or ECCC, but also whether they measure water quality, water flow, and water levels (using yes or no variables)? I have tried using the Multiple Attributes function but I'm limited to only three values.

In other words, is it even possible to do or am I stuck with Multiple Attributes?


texts: How to find upstream - lateral and downstream zones around urban areas?

texts:
I've a DTM and a shapefile of built-up areas. For a fire risk model, I've to find (within a defined buffer zone) two zones around built-up areas: one that is above/higher and lateral and the other below (downstream) the built-up area. Then I'll classify them according to model weights.

I'm using ArcGIS 10. I've made some attempts with spatial analyst but still haven't found a solution as it could be that in the same (dissolved) buffer zone an area is at the same time above a built-up area and below another built-up area I need to keep the "below" information and discard the "above" one.


texts: How do I find duplicates within a file geodatabase, and query within the groups of duplicate records?

texts:
Im currently attempting to query a file geodatabase, where I want to find all the duplicate PID (sub for an OBJECTID), and then within each one of these matching groups of records, find where the LAND_SIZE_ACRES has increased over time (each row of the same PID represents land data of sequential years).

I have been using "Select By Attributes," where I've used this formula, except it doesnt seem to work on a file geodatabase:

[FIELD_NAME] In (SELECT [FIELD_NAME] FROM [TABLE_NAME] GROUP BY [FIELD_NAME] HAVING Count(*)>1 )

Therefore I need to be able to query each group of matching records, and then find the variation within said groups of records. I just dont know where to start.

Example data:

Very few of the records actually change in land_size_acres overtime, hence the picture with no variation.


texts:
I'm working on a plug-in for QGIS and I'm a bit stuck on a problem with my dialog (realized on Qt Designer).

Indeed, I added a progression bar wich seems to work, but I can't really check this because, when I click on the "OK" button (wich launch
the execution of the plugin), my dialog just disappear.
Except this, the plug-in work well, so I was wondering if there's a way to keep my dialog - or, at least, a progression bar - displayed
during the processing, after the users have clicked on "OK".

I'd would like to show you some codes, but I'm not actually sure of what part is useful to show...

Another problem wich seems to me pretty close from the previous one : I've add to my script some error messages (QMessageBox) displayed
when the user click on "OK" when the input fields are not correctly filled.

The messages pops wells, but, when clicking "OK" on the message box, the message box AND the dialog disappear and the user have to
launch again the plug-in although I'd rather the user just back on the dialog to fill the input's field.

The code wich display the message box is this one :

if self.dlg.menu_mnt.currentText() == "" : #if combo box is empty
                    QMessageBox.warning(None, "Temperature Map", str("Veuillez indiquer un raster de MNT valide"))


texts: How can I snap 2 lines in QGIS so the st_distance between them is 0?

texts:
Sometimes when I snap the ends of 2 LineStrings together and I check the st_distance between those endpoints I get a value in the order of 1.1E-10 while at the same time when I check for equality of those endpoints I get the value 'true'.

How can I snap 2 lines so the st_distance between them is 0?

I am using QGIS 1.9.0-Master and a PostgreSQL database 9.1.7 with PostGIS 1.5.3-1.


texts: White PNG pulled from mapserver stored in PGSQL using leaflet

texts:
I have a polygon layer stored as a table within PostGIS/PostgreSQL using SRID 3857. I have created a MapServer layer pointing to this layer. However, when I try to load this layer into Leaflet as a WMS using the L.tileLayer.wms() function it only comes up as a white blank PNG over the default Leafet map area. I am pretty sure I am passing in a correct bounding box for the SRID 3857 coordinates. Here is a code snippit of what I am doing:

var map = L.map('map').setView([30.35792, -87.17677], 14);

L.tileLayer('http://localhost:8080/cgi-bin/mapserv.exe?map=C:/OSGeo4W/apps/gbp/gbp.map&SERVICE=WMS&VERSION=1.1.1&REQUEST=GetMap&LAYERS=parcels3857&STYLES=&SRS=EPSG:3857&BBOX=-31856626.987,11633877.801,-31821383.931,11657185.093&WIDTH=400&HEIGHT=300&FORMAT=image/png', {
    maxZoom: 18,
    layers: 'parcels',
    attribution: 'GBP Parcels'
}).addTo(map);

And my map file layer configuration:

  LAYER
    NAME "parcels3857"
    CONNECTIONTYPE POSTGIS
    CONNECTION "host=127.0.0.1 port=5432 dbname=gbp user=abc password=xyz"
    DATA "geom FROM parcels"

    PROJECTION
        "init=epsg:3857"
    END
    METADATA
        "wms_title" "parcels3857"
    END
    STATUS       on

    STATUS DEFAULT
    TYPE POLYGON
    CLASS
      NAME "parcels3857"
      STYLE
        COLOR 227 227 127
      END
    END
  END

Why does the layer only show up as a white PNG?


texts: Coordinates of cities and villages in The Gambia/Senegal

texts:
I am looking for sources of geo-referenced data on villages and cities (or all other things there are) in The Gambia (and Senegal if available). One issue with geocoding in rural West Afrika is the varying spelling of village names.

I own a (possibly incomplete) list of some 2000 Gambian village names, I was able to acquire some 1000 coordinates from HERE maps, Google maps, and MapQuest, through the respective APIs. But I estimate that this is missing roughly half of the villages in The Gambia. What is more. Some of these one thousand are rather bad/uncertain matches, so it would be good to compare with additional sources.

Are there any other sources, with a decent level of detail and fairly up-to-date?

Here is a list of villages that I was not able to find with any of the above sources. Since I know for sure where they are located I use those as a benchmark:

Kulukulel (district: Fulladu East): 13°17'47.1"N 14°08'29.3"W
Toubanding (district: Fulladu West): 13°28'03.7"N 14°41'09.6"W
Sami Kuta (district: Upper Niumi): 13°23'43.9"N 16°19'19.5"W

(I am aware that some of these  actually have labels with different spellings/entirely different names in google maps. However, this is so far of no help to me though, as these spellings mostly nowhere close to how the village names are usually spelled/pronounced)

Clarification in response to some comments: geonames for Gambia is in most cases rounded to minutes, which in the case of Gambia gives a grid in which places are as far as ~1.3kms, away from where they actually are. In the case of Gambian villages which are densely and rather uniformly distributed, this is not very helpful. If there are ways to obtain ungrounded geonames data I'd be happy.


texts: Change Symbology of layers from within ArcMap

texts:
In an active ArcMap session - I want to change the symbology of a layer which is repeated multiple times in the mxd. Can I do this?

The basic code I have for getting the layers is below. This is running within ArcMap (not as a .py script)

myMXD = arcpy.mapping.MapDocument("Current")
lstDataFrames=arcpy.mapping.ListDataFrames(myMXD)
lstLayers=arcpy.mapping.ListLayers(myMXD)

hazardLayer=arcpy.mapping.ListLayers(myMXD,'*Haz*')

allFrames=arcpy.mapping.ListDataFrames(myMXD)

for dataFrame in allFrames:
    myMXD.activeView=dataFrame
    for mapLayer in hazardLayer:
        mapLayer.visible=False

This just changes the visibility to off...can I change the color for example to some RGB value?


texts: Setting maximum size of HTML map tip window in QGIS

texts:
I have a shapefile with points, and for each point I have an associated image in a folder. I'm trying to display the image as map tips, and followed the simple procedure described by Nathan Woodrow in his blog.

The map tips work well, but the small window that pops up is relatively small by default, and it only shows the top left corner of the image.  It has horizontal and vertical scroll bars, but ideally I would like to show a larger pop-up window with each image automatically resized so it fits in.

Does anybody know how to do this?