texts:
I am trying to get tiles of raster contained within a bounding polygon in a shapefile. My output raster should have values for all pixels within the polygon feature, and NA or 0 everywhere outside the bounding polygon. I am trying to use gdal.warp for this. This is the code which I tried, and while it works to trim the original raster to the extent of the polygon, I am not able to make the pixel values zero outside the polygon.

InputImage = 'inputs.tif'
Shapefile = 'outcrops.shp'
RasterFormat = 'GTiff'
VectorFormat = 'ESRI Shapefile'

# Open datasets
Raster = gdal.Open(InputImage, gdal.GA_ReadOnly)
Projection = Raster.GetProjectionRef()

VectorDriver = ogr.GetDriverByName(VectorFormat)
VectorDataset = VectorDriver.Open(Shapefile, 0) # 0=Read-only, 1=Read-Write
layer = VectorDataset.GetLayer()

feature = layer.GetFeature(0)
geom = feature.GetGeometryRef()
minX, maxX, minY, maxY = geom.GetEnvelope() # Get bounding box of the shapefile feature

OutTileName = 'SomeTileName.tif'
OutTile = gdal.Warp(OutTileName, Raster, format=RasterFormat, outputBounds=[minX, minY, maxX, maxY], dstSRS=Projection, cutlineLayer=layer, cropToCutline=True)
OutTile = None # Close dataset

Raster = None
VectorDataset.Destroy()
print("Done.")

Can someone point me towards where I am going wrong?


texts:
I have a couple of .png images, they represent relief shadows, only black and white color

In ArcCatalog Preview or in ArcMap, I see that image in example, but I when I open it in Adobe Photoshop, I see only black image

Can you please help me, because I need to edit some black lines in Photoshop but I can't see picture.


texts:
I have an app built on top of Google Maps right now and everything works great.

I want to add one feature that seems exceedingly difficult. On my site, I want users to be able to type "Boston, MA" or something like that and have it find a boundary polygon for Boston and display it on the existing google map.

What is the simplest solution to this problem?

I was looking at installing PostGIS and OSM but that process seems like overkill for what I need.


texts: Displaying Dynamic Layout Title using ArcPy for ArcGIS Pro?

texts:
I have a large ArcPy script where 80 layouts (created from 1 template layout - "Layout_King") are created and exported into a PDF, but right now each layout has the same generic title so sometimes it is hard to know what data is being shown in the layout. This problem could be solved if the layout's title was Dynamic with the feature class that is being portrayed in the layout.

The following script is my first pass at achieving a layout title that would be  "Dynamic" with a feature class' name:

import arcpy
from arcpy import env
import sys
import os

elm_name = "My Dynamic Title"
p = arcpy.mp.ArcGISProject(r"C:\arcGIS_Shared\Python\CenterHeatMaps.aprx")
for lyt in p.listLayouts("Layout_King") [0]:
    for elm in lyt.listElements("TEXT_ELEMENT"):
        if elm.name==elm_name:
            text_elm=elm
            break

for m in p.listMaps():
    for lyr in m.listLayers("BCBS*"):
        new_title=lyr.name
        elm.text=new_title
        lyt.exportToPDF(r"C:\arcGIS_Shared\Exports" + "\\BCBS" + elm.name[13:] +
        ".pdf")
        print(elm.name[13:])

However, when I tested the code I got the following error:

Traceback (most recent call last):
  File "<string>", line 8, in <module>
TypeError: 'Layout' object is not iterable

What I'm doing wrong? Also, does it look like I'm at least on the right track with the code?


texts: Getting village map into correct place on OpenStreetMap raster background?

texts:
Can anyone tell me the steps to display my village map (as a tif/tfw pair of files) over an OSM raster layer (from the OpenLayers plugin) in QGIS? Because every time I try my village ends up in the middle of the ocean south of Lagos!

I know that questions like this have been asked before, and I have looked at several answers and tried to follow what they said, but so far without success. I have looked at the documentation. I understand from EPSG 3857 or 4326 for GoogleMaps, OpenStreetMap and Leaflet there are complications. One of the articles I have tried to follow (also unsuccessfully) is at http://maps.cga.harvard.edu/qgis/wkshop/basemap.php.

The tfw file is:

0.6351247107
0.0000000000
0.0000000000
-0.6350766811
532808.8175623553
273200.4324616595

I have tried combinations of:

loading the tif file with CRS EPSG:4326 and EPSG:3857
switching on and off 'on the fly' CRS transformation
changing the project CRS to EPSG:4326 and EPSG:3857
warp (reprojecting) my village file with source CRS EPSG:3857 and target CRS EPSG:4326, and then using the resulting tif file
switching between degrees and metres

I am usually good at following instructions, but this has me beaten. I thought I understood a little about about projections, but the  more I read and what I try does not work, the more I am confused.

What are the steps I should take?


texts:
I added to an original osm.map my own routes (nodes, ways, relations). I need a way to retrieve from my OSRM server these routes (for navigation).

I try to change the speed of my routes in order to give them a better priority (in lua profile) but the response from OSRM server is a mix of my route and other road.

Is there a way to retrieve (only) my preferred route from OSRM server ?


texts: I am new at geoserver. I connected to databse from admin panel and published layers. But I could not understand difference of WMS, WFS and Rest service. What is the difference?


texts: European road usage statistics per highway

texts:
I'm looking for an open data source which features the amount of cars per road per day or per year for Europe(!). I've found several national reports, which are only available as PDF and as such a not really processable. The ones by EuroStat are not really feasible. I would like to process this data into a raster to work with it in my database as part of my routing solution. Any suggestions?

P.S.: If there is already a rasterized solution out there, I would certainly prefer that.


texts: Create "grid" of polygons with the same 3D area based on an underlying raster variable (e.g. DTM height)

texts: This could apply to point or raster values but imagine a raster of a continuous variable such as height (DTM). Is there a way to create a complete coverage of polygons that will have an equal area, taking into account the roughness of the underlying terrain? I realise that it would depend on where you seed the polygons (i.e. no unique solution) but if we accepted a certain degree of error it should be possible to do this with a nested regular grid (i.e. squares of differing fixed sizes)?


texts: raster format batch conversion in R, error gdal_translate

texts:
I am running a batch conversion, .grd to .tif. This code worked fine with .rst to tif but now I get the error:

  /ndmi.LC81270592013272.grd' not recognized as a supported file format.

This is my code:

batch_gdal_translate(infiles=dir,outdir=output_folder,
                     outsuffix="_conv.tif",of="GTiff",pattern=".grd$", ot="Int16")


texts: Reduce image using two image inputs - one with zones the other with data to summarize

texts:
I would like to reduce an image where the inputs are two images and the output is a table. One image will have zones (integer values) and the other data to summarize. For each zone in an image, I would like to calculate the mean of values in the other image that overlap that zone. In ArcGIS one could use the tool Zonal Statistics as Table.
https://pro.arcgis.com/en/pro-app/tool-reference/spatial-analyst/zonal-statistics-as-table.htm
//1- zone image (DEM binned 100m intervals)
//2- data to summarize (mean minimum temperature)
//3- figure out how to summarize mean tmin per dem zone and export as a csv

//zones (dem binned x 100m)
var DEM_binned = image.divide(100).ceil().clip(geometry);
Map.addLayer(DEM_binned, {min:0, max:40, palette:['blue', 'purple', 'green', 'yellow', 'orange', 'red', 'white']}, 'dem');

//data to summarize (mean long-term annual tmin)
var tmin = imageCollection.select('tmin');
var tmin = tmin.mean().clip(geometry);
Map.addLayer(tmin, {min: -5, max: 20, palette:['white', 'blue', 'green', 'yellow', 'red']}, 'tmin')

//get mean tmin per zone as a table

GEE working example script: https://code.earthengine.google.com/0f50cbfbbb8ee8eabad07d9569b9a58a


texts: Let's say I want to move my 2D coordinate system by A and B meters in pyproj. How would I do it?


texts: How to improve OSM layer resolution in OpenLayers plugin while exporting?

texts:
I'm having trouble with the OpenLayers Plugin exporting a map to a png or pdf file:

The output resolution of the exported osm-layer is very, very poor. Zooming to 100% in the png file shows big pixels. The rest of the map - so the output from shapes and fonts - is in the usual high resolution e.g. 300 or 600 dpi.

Is there any way to improve the osm-background-resolution?


texts: Filling null attributes with mean attributes of neighboring polygons using ArcGIS Desktop?

texts: I have a layer of polygons, where in some polygons there is missing data (null). I would like to fill these missing attributes by the average attributes in neighboring polygons. Is there a simple way to achieve this?


texts: How to remove vehicles from an aerial photograph

texts: I have an aerial photograph of an area. The image includes visible traffic signs and vehicles. How can I remove the vehicles from the image? Is it possible to do so in software like Erdas or Envi using methods like Density Slicing? If possible, please tell me the procedure.


texts: Exporting symbology from ArcGIS Desktop to ArcGIS Online?

texts:
We're planning to put up a bunch of our maps into ArcGIS Online and we'd like to create a streamlined process to do so, with as much automation (Python) as possible. This is our first test with a map of a city park with lots of unique symbols. We figure we'd get a map up there, then figure out what can be scripted.

Anyhow, so, we have an .mxd of a park where pretty much each layer has fancy symbols for things like parkbenches, bathrooms and the like. Upon fiddling around with ArcGIS Online, it looks like it only takes .txt, .csv, and .gpx files. I managed to get layers into "My Content" but it doesn't let me add it into a map. When I go "Add>Search for Layers" and tell it to search in My Content, nothing shows up. I then shared the layer in My Content so it was public on ArcGIS Online, and tried seeing if the search func. would see the layer in ArcGIS Online. No luck.

How do I get it to take attribute data, but also the symbology data from my layers?


texts:
I am using OpenLayers, Geoserver 2.1 CQL language. The problem is that I do not know how to address the layer. I am trying to call a function like this:

   var createWmsLayer = function(name, url, params, options) {
        map.addLayer(new OpenLayers.Layer.WMS(name, url, params, options));
    };

I am creating the layer like this:

var ikke_listet = new  createWmsLayer("Ikke rødlistet",
     "http://kart.naturkart.no/geoserver/wms",
                {layers: 'gbif_ns:ikke_listet',format: 'image/png',
                transparent: 'true'}, {singleTile: true, ratio: 1,
                  isBaseLayer:false, visibility: false});

Then I am calling a function on the layer to redraw it with CQL params to GeoServer.

function CQLfilter(param){
        filter_list = []
        filter_list.append("collector LIKE '%H%' ")
        ikke_listet.mergeNewParams({'CQL_FILTER': "collector LIKE '%H%' "})

But I get this error in Firebug:

ikke_listet is not defined

So my question is: What is the name of my layer?

I am using the "name" part of my function to write out some styles for the layerlist (I use ExtJS) and put a padding and different icons in the layerlist to show what icons represent which layer. Is this a wrong way to do this styling?

Any suggestions as to what I use in the function to add CQL call to the layer and redraw it?

Aragon, I was able to get the layername from your suggestion of the layer array:

function CQLfilter(param){
var mLayers = map.layers;
mLayers[9].mergeNewParams({'CQL_FILTER': "collector LIKE '%H%' "})
// This gave me the names of the layers
// for(var a = 0; a < mLayers.length; a++ ){
//    alert(mLayers[a].name + '#' + a + '#')
//};
}

So I wonder if there is another way to style my layer list because now I get the layernames like this:
 "Ikke rødlistet". But of course I can still get to the layer by referencing the layer array position.

I managed to run the function now and it filters correctly. I only had to ignore the first base 4 layers and I have to find out how to collect the values from the search form. But manually this code now works and I have much progress :-)

function CQLfilter(param){
var mLayers = map.layers;
// have to get the search by document.getElementById or jQuery I guess. Testing manually now
var search = "kingdom LIKE '%Plantae%' ";
for(var a = 4; a < (mLayers.length); a++ ){
    alert(mLayers[a].name + '#' + a + '#')
    mLayers[a].mergeNewParams({'CQL_FILTER': search})
};
}

Hi again, I am able to get the layer now. Works fine with this function on the CQL filter and the reset filter function at the end.

function CQLfilter(sok,cat){
var mLayers = map.layers;
switch (cat)
{
case 'SN':
   param = "scientificname LIKE '%"+ sok + "%' ";
  break;
case 'NNG':
   param = "norsknavngruppe LIKE '%"+ sok + "%' ";
  break;
case 'NNA':
   param = "norsknavnart LIKE '%"+ sok + "%' ";
  break;
case 'PROJ':
   param2 = new Array()
   param2.push("prosjektnummer LIKE '%" + sok + "%'")
   param2.push("prosjektnavn LIKE '%" + sok + "%'")
   param = param2.join(" OR ");
  break;
default:
   param2 = new Array()
   param2.push("scientificname LIKE '%" + sok + "%'")
   param2.push("norsknavngruppe LIKE '%" + sok + "%'")
   param2.push("norsknavnart LIKE '%" + sok + "%'")
   param2.push("collector LIKE '%" + sok + "%'")
   param = param2.join(" OR ");
  break;
}
for(var a = 4; a < (mLayers.length); a++ ){
    mLayers[a].mergeNewParams({'CQL_FILTER': param});
    mLayers[a].setVisibility(true);
}}
function resetCQLfilters() {
    var mLayers = map.layers;
    for(var a = 4; a < (mLayers.length); a++ ){
    delete mLayers[a].params.CQL_FILTER;
    mLayers[a].redraw();
}}


texts: How can I process vector difference between 1 group layer and 1 polygon layer?

texts:
I have one layer group called 'left_buffer', containing 40 buffered polygon layers (name: 1,2,3,...,39,40). I want to use difference vector tool between 'left_buffer' and 'bldg' layer, but I'm getting error message.

for i in range(40):
    processing.run("qgis:difference", {'INPUT':
    'left_buffer'[i+1], 'OVERLAY':'small_bldg',
    'OUTPUT':'C:/Users/sky/Desktop/1/' + str(i+1)
    + '.shp'})

error message: Could not load source layer for INPUT: e not found

How can I fix this issue?


texts: QGIS - How to differently colorize labels from two separate fields using color expression script?

texts:
I have 2 fields, I concatenate them [concat(to_string(field1), ' ', to_string(field2))] to show their labels side-by-side.

Is it possible to colorize from "text color expression" (using some scripts) first label in "color1" and second label in "color2"?

Clarification

I want the final result to be like below ("field1" data to be shown in red, "field2" data to be shown in green)


texts:
I want to make an ee.Image from an ee.Dictionary of values.

I generated a dictionary using the following piece of code that takes a sample pixel's data from an image collection and generates a dictionary of values where each value of the pixel in the image collection(monthlyNDVI) is now an entry in the dictionary (created a time-series for a given pixel).

var reference_signal_list = monthlyNDVI.iterate(referenceTimeSeries,
ee.List([]));
print('Reference Signal', reference_signal_list);

var referenceTimeSeries = function(current, previous) {
var temp = current.reduceRegion(ee.Reducer.first(), reference_point, 30);
return ee.List(previous).add(ee.Number(temp));
};

I now want to create an image from this dictionary in a way that each entry in the dictionary becomes a band in the resultant image.


texts: Parsing multivalue parameter using ArcPy?

texts:
I am trying to get my head around how multivalues derived from the arcpy command GetParameterAsText gets looped through a script. Could some please provide a sample script that does the follow basic task so that I may study it for further insight?

Get a set of shapefiles from arcpy.GetParameterAsText(0) (the number of shapefiles can vary from to 10 say) and then loop through each shapefile, displaying the following with arcpy.AddMessage saying something to the effect:

  "This is shapefile " shapefile name

When looped through all of the shapefiles, print the message "No more shapefiles"


texts: Extract Certain Elements from Global OSM Database

texts:
The goal is to find all waterways in Open Street Maps globally.

This is the equivalent for a random area in China in overpass turbo. Instead of doing it for this bounded area, I want to query globally.

Overpass turbo gets overloaded by a global query but I found osmosis in a related related2 question which seem to be able to get the job done.

wget -O- https://planet.openstreetmap.org/planet/planet-latest.osm.bz2 | bzcat | osmosis --rx - <filter args here> --wx my-compressed-output.osm.gz

I struggle however in translating the overpass query to the osmosis query and got lost in osmosis documentation.

How would I filter the waterways in osmosis? I think I need --way-key (--wk) or --way-key-value (--wkv) but like to check before executing a massive query.(50GB download)

Alternative solutions (non osmosis) are welcome too.


texts: DEM produces pixelated hillshade when CRS changed to WGS 84 / UTM

texts:
UPDATE:
I now realize my problem with pixelation occurs when I change the CRS to WGS 84 / UTM. Anyone familiar with why this would happen or how to avoid it?
ISSUE:
I am having trouble with hillshade on a QGIS proejct I have been working on for some time. The hillshade I get using a DEM comes out with significant pixelation (See image below). I've used the same DEM in a separate project and I don't get the pixelation at all. I would start a new project except that I have very elaborate Map Themes set up in the original (problematic) project that I can't seem to transfer to a new project.
A stripped down version of the project QGIS files and layers can be found here:
https://drive.google.com/file/d/1Kt6TwXI2SBzjjT21kUr2MJMXcTVV3YEp/view?usp=sharing
This is the process I use to generate the hillshade in a new project where I don't get pixelation:

Open a new QGIS project
Add the DEM as a raster file
In symbology

Set render type to Hillshade (directly from the dem)
Set resampling section to 'cubic' in the Zoomed in field.
set transparency to  50%
Click  Apply

If I use the DEM with the above process I get smoothed hillshade from the DEM (no pixelation) but the same DEM produces a pixelated hillshade in my original project. I've tried eliminating all other layers to see if it's a layer conflict but the hillshade is still pixelated even when it's the only layer. It seems to me that there must be some setting or preference in the old project (linked above) that prevents me from smoothing the hillsahde using resampling.
Anyone have any idea how to fix my original project so resampling can be used to smooth the hillshade?
Here's a sample of the map with pixelated hilshade.


texts: Geocode quality: Nominatim vs. PostGIS Geocoder vs. Geocoder::US 2.0

texts:
I occasionally need to geocode >> 500,000 US street addresses (no features or intersections). I'd like to set up a geocoding service locally to handle these bulk geocodes.

I've used the Tiger-based Geocoder::US 2.0 locally in the past. Initial benchmarking gives reasonable accuracy against Google-geocoded addresses, especially when I don't rely on the Ruby parser and I pass the address elements to the geocoder explicitly.

I've also read a good deal about the PostGIS Geocoder. Since my data are already in PostGIS this sounds attractive, but I haven't set it up locally to try it out.

Also, someone recommended I try Nominatim locally. I don't know much about it except that it's not Tiger based (as far as I can tell), and I'm wondering if I'll get more reliable results.

My question to anyone with experience with the above solutions: How does the accuracy of the three compare? Specifically, which provides the most consistently-reliable geocodes for US street addresses? I'm not overly concerned with speed or set up complexity.


texts: Seeking top tier conference in GIScience?

texts:
In many fields, such as computer sciences, it seems to be quite clear which conferences are top tier and how a paper accepted to those conferences compares to journal publications in the same field.

Is there a general agreement in GIScience about which conferences are considered top tier? And how do publications at those conferences compare to papers in, e.g. Transactions in GIS or the International Journal of Geographical Information Science?


texts: Generate new geometric network from intersection between curent network and polygons in ArcGIS Desktop?

texts:
I have generated a stream network using the "Strea Reach and Watershed" tool from the TauDEM toolbox.

Within the stream reach dataset there are fields that give an ID to the stream segment and fields which state the upstream and downstream segments, i.e. a network.

I have a polygon shapefile with field boundaries. What I would like to do is create a point wherever the stream network crosses a field boundary and then join the points together with straight lines to represent links between the fields based on the stream network. I have managed to create points where intersections between the network and the field boundaries are, I am just unsure how to join these together with lines that follows the network.

I have illustrated what I mean below with crudely drawn red lines. Figure 1 is the basic form of what I would like to achieve i.e. connections between the field boundary intersections along the network.

Ideally I would like to achieve what is in figure 2 which is lines that follow the network stemming from the centroids of each field polygon, i.e. one line can represent all flow lines crossing into an adjacent field.

Figure 1 - Field links based on intersections with field boundaries

Figure 2 - Field links from centroid of polygon

As I stated I've used the TauDEM tool box in ArcGIS 10 so far, but am open to any other tools in ArcGIS Desktop.


texts: What kind of coordinate uses this WMS? Also how to convert to lat/lon

texts:
On this map http://infomapa.rosario.gov.ar/emapa/mapa.htm which I believe is MapServer, I get some coordinates after inspecting with Chrome like this:

    "geoJson" : "{\"type\":\"Point\",\"coordinates\":[5433059.165834456,6360608.810436543]}"
}, {
    "numero" : 3004,
    "linea" : "106",
    "geoJson" : "{\"type\":\"Point\",\"coordinates\":[5433151.161165902,6360413.851751617]}"
}, {

1 - How do I know which type of coordinate are these?

2 - How do I convert to lat, lon?

As an example I know:

   "geoJson" : "{\"type\":\"Point\",\"coordinates\":[5434921.2799996715,6359301.360077026]}"

should be almost equal to this

-32.907269, -60.695515


texts:
I am trying to use the aforementioned plugin with QGIS 1.7 and Postgres 9.0/PostGIS 1.5.

Added a PostGIS spatial layer w/out a problem and then clicked on the RT Sql layer button which displayed the main dialog. When I select Database-> , I got this error

BTW - the ...'s is the QGIS plugin path. Left it out for clarity):

An error has occured while executing Python code:

Traceback (most recent call last):
  File ".../python/plugins\rt_sql_layer\ManagerWindow.py", line 112, in dbConnectSlot
    self.dbConnect( self.sender() )
  File ".../python/plugins\rt_sql_layer\ManagerWindow.py", line 163, in dbConnect
    self.refreshTable()
  File ".../python/plugins\rt_sql_layer\ManagerWindow.py", line 272, in refreshTable
    model.loadFromDb(self.db)
  File ".../python/plugins\rt_sql_layer\DatabaseModel.py", line 250, in loadFromDb
    self.tree.constructTreeFromDb(self.db)
  File ".../python/plugins\rt_sql_layer\DatabaseModel.py", line 96, in constructTreeFromDb
    list_tables = db.list_geotables()
  File ".../python/plugins\rt_sql_layer\connectors\postgis.py", line 275, in list_geotables
    items[i] = geo_item
IndexError: list assignment index out of range

Python version:
2.5.2 (r252:60911, Feb 21 2008, 13:11:45) [MSC v.1310 32 bit (Intel)]

QGIS version:
1.7.0-Wroclaw Wroclaw, 63ecdd7

I never used python in the past but I had a look at the source code and it seems that the query the plugin is running to get the spatial tables into a tree view is not returning anything? Had a go to running the sql statements in PostGIS and they all seem to return data. Then again- not sure if I am looking at the right place!

Have you seen this before? Has it got anything to do with the s/w versions? Plugin was developed for an earlier version and something has changed in the latest?


texts: Google Earth Engine band did not match any bands errord

texts:
According to the dataset, there should be a band named: total_precipitation_surface. However, when I try to get data from this band I get:

> Collection.first: Error in map(ID=2019061000F000): Image.select:
> Pattern 'total_precipitation_surface' did not match any bands.

I works perfectly with other bands in this dataset, but for some reason it does not work for the band I want.

My code is:

var fc2 = ee.Geometry.Point([-159.5469188,21.9080489]);

Map.centerObject(fc2,10); //the last numbe is zoom level

print (fc2)

var data_collection = ee.ImageCollection('NOAA/GFS0P25');
var dataset = data_collection.filterBounds(fc2)
                    .filterDate('2019-06-10', '2019-06-15')
                    .select('total_precipitation_surface');

// Create a graph of the time-series.
var graph = ui.Chart.image.seriesByRegion({
  imageCollection: dataset,
  regions: fc2,
  reducer: ee.Reducer.mean(),
  //seriesProperty:'uniqe_name'
})
print(graph)


texts: QGIS SVG with editable color and project var value

texts:
in QGIS we can make SVG colors editable by adding params:

<svg width="100%" height="100%">
<rect fill="param(fill) #ff0000" stroke="param(outline) #00ff00" stroke-width="param(stroke-width) 10" width="100" height="100">
</rect>
</svg>

The optional parameter behind is the default.

Is it possible to have project vars as defaults? If so how is the correct syntax?

<svg width="100%" height="100%">
<rect fill="param(fill) @myfillcolor" stroke="param(outline) @myoutlinecolor" stroke-width="param(stroke-width) 10" width="100" height="100">
</rect>
</svg>


texts:
So I'm new to QGIS (2 days) but have spent HOURS on it everyday trying to learn it to process raw SRTM data make 3D relief maps on another software, so I can work my way around it but can't really do any of the coding. Therefore, please consider that if your answer will refer to coding. The files I'm using are here
I downloaded the SRTM images (98 tiles), merged them into one, reprojected them, and now I want to use the following raster calculation on them: (Pixel Value + 23) ÷ 1244 * 65535. 65535 is because it's a 16-bit TIFF and the other values have to do with the range of pixel values in the merged STRM map. This is to avoid terracing in the 3D relief map when created in the other software and to get rid of negative values.
When I try and run the calculation, it gives me this error:
"The following layers were not correctly generated.
• C:/Users/dell/OneDrive - UOM/Blender/Map Making/QGIS Files/UK/SRTM Map Scaled.tif
You can check the 'Log Messages Panel' in QGIS main window to find more information about the execution of the algorithm."
The calculation would work and create a raster if I'm using one of the downloaded SRTM files as the 'reference layer', but that obviously doesn't give the right extent and projection. It gives the error I copied when using the reprojected, clipped or merged layers. I've tried doing the calculation in a new project file with no luck.


texts: Merging attribute and geometric features in QGIS?

texts:
I have polygons shapefile in which some of the polygons have same attributes as seen below:-

In QGIS, how do I merge the polygons having the same attributes into one polygon?

In other words, am looking for a way to merger polygons/attribute of the same layer, that is similar way of merging vector layers found at: Vector >> Data Management Tools >> Merge Vector Layers


texts: Is there a way to use the over function in r to retrieve all polygons that do NOT have points in them?

texts:
I have used the over function to retrieve the polygons of a map that contain points. I'am wondering, for my own learning purposes, if there is a way to use the over function to retrieve the polygons that don't have points in them. Is this possible?

Here is some code to return polygons of a map WITH points in them:

coordinates(data) <- ~ long_xcord + lat_ycord

# Set the projection of the SpatialPointsDataFrame using the projection of the shapefile.
proj4string(data) <- proj4string(shapeData)

p <- over(data, shapeData, returnList = TRUE)

p <- do.call(rbind.data.frame, p)

can the above be modified to return the polygons that DO NOT have points in them?


texts:
I want to classify (and separate) houses from trees (high vegetation) in point cloud data using PDAL.

I use height filter for this purpose. Our high Vegetation lies in the range between 5 and 10 meters. So, I apply filter of range 5m to 10m.

The filter output gives me few houses and all of trees.

I think we can differentiate these remaining tree from houses using tree canopy width; however, is there any filter which can differentiate trees from houses?

Here is my filter details:

{
  "pipeline":[
    "./ground_non_ground.las",
    {
      "type":"filters.hag"
    },
    {
      "type":"filters.ferry",
      "dimensions":"HeightAboveGround = Z"
    },
    {
      "type":"filters.range",
      "limits":"Z[5:10]"
    },
    {
      "type":"writers.las",
      "filename":"./tree.las"
    }
  ]
}


texts: Issues creating ArcSDE sde_logfile for certain users (username begins with a numeral)

texts:
I am having issues when editing an sde feature class from certain users. The way my database is set up is that the authentication roles are updating constantly by comparing the current authentication roles with an LDAP server.

Everything seems to be working fine except for people with user names beginning with a number. I get an error that says "unable to create logfile system tables. User perhaps lacks permissions or resources to create tables [ERROR: syntax error at or near "6" Line 1: CREATE TABLE 6zn.sde_logfiles (logfile_name VARCHAR....)

If the syntax read CREATE TABLE "6zn".sde_logfiles... it would work just fine. Seems like a simple (but scary) solution would be to edit sde functions/triggers/etc. to add quotes for arguments that start with a number in the CREATE TABLE but I have no idea where to start...

Unfortunately, there is no way to get around user names that don't start with numbers.


texts: How to return zero's if ee.Reducer.count() returns empty image?

texts:
I have an algorithm that is computed using 2 ImageCollections over many areas and timeframes. When the algorithm is finished I export the results in an image, and I add 2 bands that tally the number of pixels from each ImageCollection the final result was based on. I do this using ee.Reducer.count() as follows:
var s2toa = ee.ImageCollection("COPERNICUS/S2")
                  .filterBounds(roi)
                  .filterDate('2015-01-01', '2016-01-01');

var nPixels = s2toa.select('B2')
                   .reduce(ee.Reducer.count())
                   .rename('nPixels');

print(nPixels);

However, it occasionally happens that no images are available at a given location within the timeframe of interest. In these instances the above approach returns an empty image:
var s2toa = ee.ImageCollection("COPERNICUS/S2")
                  .filterBounds(roi)
                  .filterDate('2013-01-01', '2014-01-01');

var nPixels = s2toa.select('B2')
                   .reduce(ee.Reducer.count())
                   .rename('nPixels');

print(nPixels);

In these instance, I would like to replace the empty image by an image of zeros, so that I can log the fact that zero pixels from that collection were used. Is there a way to do this programmatically? Perhaps there is something easy I have not thought of...


texts: Raster symbolization in open source web gis?

texts:
We would like to build a simple web map application that displays raster data (e.g. a digital surface model) and symbolizes values (e.g. height) using a color ramp.

Is it possible to have the user manipulate the symbolization, maybe even allow to perform a classification?

For example like this: http://maps.nrel.gov/biomass  - only that we would like to additionally change the classification on the client side.

I am aware of WCS and read about GDAL in this answer, but it did not help much.

Can anyone point me to some actual examples doing raster manipulation?


texts: Seeking Transit Feed (GTFS) Toolset for QGIS?

texts:
I can see that ArcGIS has a transit feed toolset that enables converting the .txt files from GTFS data directly into usable feature classes for GIS.

I am trying to develop catchment areas in the Greater Metropolitan Area of Montreal accessible within 20, 40 and 60 minutes.

Is there a toolset like the ArcGIS one (see link below) that I can use in QGIS?

https://pro.arcgis.com/en/pro-app/tool-reference/conversion/an-overview-of-the-transit-feed-gtfs-toolset.htm


texts: QGIS - Network join issue with cost allocation using OSM dataset

texts:
I've been trying to create a cost allocation map, using the v.net.alloc tool on QGIS with Grass.
Some of the results don't connect to the network fully, even though they should in real life. This has occurred for multiple points in different locations I've used the v.net tool to connect the nodes to the network lines. Any suggestions?


texts: Comparing field values of two feature classes using ArcPy?

texts:
I have an intermediate fc #1 that contains data that I would like to append to my production fc # 2(final output).  When working with the data, I would only like to append data from my fc # 1that is not in fc  # 2, to avoid duplicates, etc.  My field that I would like to check is SRNumber I would like to only append from fc # 1 to fc # 2 when SRNumber is not in fc # 2.  Can this be done via ArcPy?

Maybe get values of fields from fc # 1 SRNumber store in list and write anything not in list to fc # 2?


texts:
I am using ST_Intersects() to know whether two geometries intersect or not. This is my code:

$var_name = pg_query($conn,"SELECT ST_Intersects(p1.column_name,
p2.column_name)
FROM table_name AS p1, table_name AS p2
WHERE p1.id < p2.id");
if(pg_fetch_assoc($var_name) == 'false')
{
echo "accepted";
}
else
{
echo "rejected";
}

I know that the geometries are not intersecting and when I run the same query in Postgres the result is false as expected.There are only two rows in the table.

But here it is executing the else statement. Where am I wrong?


texts: I am looking for annual precipitation maps for Minnesota and South Dakota. Preferably in a raster or polygon format. Basically, something I can spatially tie to data in a large database. I have been on NOAA, USGS, and several other sites that typically have this type of data, but I'm not having much luck. Any ideas?


texts: What does Scale mean in the reference of Desktop GIS or WebGIS?

texts:
I'm having trouble understanding the meaning of the concept of scale when it comes to maps on screen, i.e. either in a Desktop GIS or in a WebMap.

When it comes to paper maps, the idea is very clear. The Scale is the ratio of a particular distance on the map, to the same distance on the ground.

But when it comes to maps on screen, the same concept cannot be applied. Firstly different screens do not have the same pixel size. It might be that a 14 inch monitor and a 24 inch monitor could have the same resolution (in terms of pixels). Then comes the issue of dpi. An iOS device with retina display has a dpi several times that of my work-station monitor. So 1 cm on my screen might show a different distance than 1 cm on your screen.

So what does the scale of a map in digital cartography mean? What does it mean, when I see the scale of a map is 1:180,000 in ArcMap?

Aside
This question came about because of several ArcMap Documents I received from a client who publishes paper maps.
The paper map is published at 1:180,000, but at that scale the map looks terrible. It's only when I zoom in, does the map look good. With trial and error, I found that the map looks best on screen at about 3 times the scale, i.e. at 1:60,000


texts: Update points geometry while snapping them to nearest line

texts:
I have two tables in the PostGIS database - points and lines. And I want to update points geometry to snap them to line if they lay within a certain distance.
By now I managed to create a new table with snapped points which has only geometries. This works:
CREATE TABLE new_table AS
SELECT ST_ClosestPoint(lines.geom, points.geom) as snapped_point
FROM lines
INNER JOIN points on ST_Dwithin(lines.geom, points.geom, 5)
ORDER BY points.id;

But what I really need is to update the existing table. So I try:
UPDATE points
SET geom =
(SELECT ST_ClosestPoint(lines.geom, points.geom) as snapped_point
FROM lines
INNER JOIN points on ST_Dwithin(lines.geom, points.geom, 5)
ORDER BY points.id);

And it did not work. How do I put these two pieces together?


texts: How to use Google Maps geocoding with Openlayers?

texts: I'm making a Openlayers map with some shapefile and I have Google Maps base layer. I want to create a text area and a button. When I input a address and click the button, the map should show me the google maps address in the map. Can you help me please ?


texts: Fetching GPS point data from a REST service endpoint using jQuery/AJAX and displaying them as tracks on an OpenLayers map

texts:
How do I structure this ajax call to receive GPS point data from a REST service point and display the points as tracks on an OpenLayers map. So far I have the following code;
dataProvider.requestJSON = function(trackId) {
    $.ajax({
        url: "https://services.arcgis.com/Sf0q24s0oDKgX14j/arcgis/rest/services/gpsData/FeatureServer/0/query?where=Source= " + trackId + " &outFields=*&f=geojson",
        method: 'GET',

        success: function(response) {
            visualizer.sendDataToMap(response);

        },
        error: function(error) {
        }
    });
};

The idea is to have the user input the trackId on the application and visualize the respective track. The trackId is called source on the GeoJSON. How can I fetch the GPS data from the service endpoint and display them on an OpenLayers map as tracks? I would like to implement the solution with jQuery/AJAX, AngularJS and OpenLayers.


texts: I'm trying to draw a pie chart on QGIS, over a polygon layer. This pie chart should include data from 6 fields in the attribute table. These fields contain the percentage of several categories of data and all of them sum up to 100% (I've checked it personally). However, the diagrams option in layer properties gives no result. Anyone knows what may be going on?


texts:
I use Google Map Tiles (lvl 16 to 25). But, I can't reach more than 22 maxZoom level when I use the maptype Roadplan and 20 when I use the maptype satellite..
I would like to increase the maxZoom to 25 but it doesn't work.

<!DOCTYPE html>

        <title>Orthomosaic</title>
        <meta charset="utf-8" />


texts: Capping the amount of a certain cost per route with Network Analyst

texts: With the Network Analyst extension in ArcGIS 10.2, is there a way to set the maximum amount of a certain cost that a route must not exceed? For example, can I require routes to be below a certain length? Or, if a route makes several stops incurring some cost (i.e. picking up some packages or containers) at each stop, can I route it such as not to exceed the number of items picked up for all stops it visits?


texts: Dynamically move map canvas and main window using PyQGIS

texts:
I am trying to create plugin that will open new window with some layer, Dynamically "pan" my window, when main window is pan. So boundingbox coordinates should be the same in both windows.
I need new window to move with main QGIS window. I've found out in documentation useful code for map canvas:

import sys, os
class MyWnd(QMainWindow):

    def __init__(self, layer):
        QMainWindow.__init__(self)

        self.canvas = QgsMapCanvas()
        self.canvas.setCanvasColor(QtCore.Qt.white)
        self.canvas.setExtent(layer.extent())
        self.canvas.setLayerSet([QgsMapCanvasLayer(layer)])
        self.setCentralWidget(self.canvas)

        actionZoomIn = QAction("Zoom in", self)
        actionZoomOut = QAction("Zoom out", self)
        actionPan = QAction("Pan", self)
        actionZoomIn.setCheckable(True)
        actionZoomOut.setCheckable(True)
        actionPan.setCheckable(True)

        self.connect(actionZoomIn, QtCore.SIGNAL("triggered()"), self.zoomIn)
        self.connect(actionZoomOut, QtCore.SIGNAL("triggered()"), self.zoomOut)
        self.connect(actionPan, QtCore.SIGNAL("triggered()"), self.pan)

        self.toolbar = self.addToolBar("Canvas actions")
        self.toolbar.addAction(actionZoomIn)
        self.toolbar.addAction(actionZoomOut)
        self.toolbar.addAction(actionPan)

        # create the map tools
        self.toolPan = QgsMapToolPan(self.canvas)
        self.toolPan.setAction(actionPan)
        self.toolZoomIn = QgsMapToolZoom(self.canvas, False) # false = in
        self.toolZoomIn.setAction(actionZoomIn)
        self.toolZoomOut = QgsMapToolZoom(self.canvas, True) # true = out
        self.toolZoomOut.setAction(actionZoomOut)
        self.pan()

    def zoomIn(self):
        self.canvas.setMapTool(self.toolZoomIn)

    def zoomOut(self):
        self.canvas.setMapTool(self.toolZoomOut)

    def pan(self):
        self.canvas.setMapTool(self.toolPan)

And using class below I get coordinates from main canvas and set coordinates to my canvas (part of QGIS plugin):

def run(self):
    """Run method"""
    import re

    wnd = MyWnd(self.iface.activeLayer())
    wnd.show()
    canv = qgis.utils.iface.mapCanvas()
    coords = canv.mapRenderer().extent().toString()
    coords = re.findall("([0-9]*\.[0-9]*)", coords)
    wnd.canvas.setExtent(QgsRectangle(float(coords[0]), float(coords[1]), float(coords[2]),float(coords[3])))
    wnd.canvas.refresh()

I think this part of code can be used in loop (if coordinates change - I will set new coordinates in own canvas). But I don't sure if it is a good way.
Is there any other way to solve this issue? Provide me with right way, please.


texts: CalculateField_management - updating fields with X and Y coordinates for start and end of line

texts:
I've searched all the help documents and can't see where my script is going wrong.

I've tried adding the fields as FLOAT, replacing "PYTHON_9_3" with just "PYTHON" but it just doesn't work. The field get created fine.


texts: Changing graduated color class values with ArcPy?

texts:
I am trying to use python to adjust my graduated color symbols' class breaks and labels. The layer is a query layer tied to an SQL DB and I want to be able to execute my python (after the database has been updated) to automatically author a map and have the layer/legend reflect the changes made in the database. I've tried this script but it doesn't work and I get an error I don't recognize.

import arcpy
mxd = arcpy.mapping.MapDocument("current")
df = arcpy.mapping.ListDataFrames(mxd, "Layers")[0]
lyr = arcpy.mapping.ListLayers(mxd, "Unsafe_Driving", df)[0]
if lyr.symbologyType == "GRADUATED_COLORS":
    lyr.symbology.valueField = "CSA.DBO.%Unsafe_Driving_1.VARIANCE"
    lyr.symbology.classBreakValues = [-.01, .01, .0866]
    lyr.symbology.classBreakLabels = ["-8.66% to -1.00%", "-.99% to 1.00%","1.01% to 8.66%"]

arcpy.RefreshActiveView()

  Runtime error
  Traceback (most recent call last):
    File "", line 7, in
    File "c:\program files (x86)\arcgis\desktop10.1\arcpy\arcpy\arcobjects_base.py", line 87, in _set
      return setattr(self._arc_object, attr_name, cval(val))
  RuntimeError

The goal here is to have 3 classes:

Lowest value to -1.00%
-.99% to 1.00%
1.01% to the highest value.

Since the lowest and highest values can change, do I need to use a query to dynamically insert those into the lyr.symbology.classBreakValues?


texts: PostGIS shapefile import manager error: syntax error

texts:
I'm trying to use PostGIS Shapefile Import/Export Manager to import a shapefile that was created in QGIS. The shapefile has all the blocks of the 2010 U.S. Census for Oregon and was exported as a polygon from a TIGER database that came from the U.S. Census website.

I was having one problem with the shapefile not importing and fixed that with the first answer from this question, i.e. I needed to create a postgis extension in the PostgreSQL database I wanted to work in.

Now, I get a new error when I click on "Import". I get the following error message:

Importing with configuration: or_census_popdens, public,  , D:\data\data\PostgreSQL\data\OR_census_popdens.shp, mode=c, dump=1, simple=0, geography=0, index=1, shape=1, srid=4326
Shapefile type: Polygon
PostGIS type: MULTIPOLYGON[2]
Failed in pgui_copy_start(): ERROR:  syntax error at or near ")"
LINE 1: ...","area_ha","popdens_ha","area_sqmi","popd_sqmi", ) FROM std...
                                                             ^

Shapefile import failed.

It seems there is some kind of end of file problem with the attribute table? In the attribute table itself, "popd_sqmi" is the final column, but the code for the error shows a "," trailing that column and just before the closing ")" in the command . I don't know exactly what is causing this error though. How should I fix this?


texts: Symbolizing individual highway network segments in ArcMap

texts:
I have a highway layer that I need to style in ArcMap in such a way to differentiate between individual highway sections, ideally using no more than 6-7 distinct colors.  Each highway is divided into one or more sections.  Highways have Road Number and Section Number attributes that I can use for styling.  The symbology will need to be exportable to ArcPad.  I can calculate/add new fields to the table, if that helps.

Here's an image of how a section of highway may look like:

Note that in the image, the colors aren't repeated, but they will definitely have to be over the entire highway network. All I'm looking to accomplish is show highway section breaks with a limited subset of colors.


texts: How to add/create features in postgis layer without geometry?

texts:
I have a Postgis layer loaded in QGIS 2.8.6. There is empty POINT geometry field that I have to fill with geographic information.

How do I add features to this layer without adding new geometries? My intention is to complete the point field geometry with coordinates information.

For exemple, for the row 0 I have all the information but I have to complete it with the geometry information (In this case, UTM information).


texts: Distance between two points is not proportional above or below the Equator

texts:
I am trying to develop a tool that would find a point (longitude, latitude) that would be located a chosen distance away to the left from given point using Python and pyproj. I am not very familiar with geodesic calculations and I have come up with an idea to use proportion to calculate the (long,lat) of the distant point using a given point and a reference point, that has the same latitude but different longitude:

from pyproj import Geod

def get_point_placed_distance_away_from_point(point, reference_point, distance):
    geod = Geod(ellps="WGS84")
    point_located_by_distance_away = {"type": "Point", "coordinates": None}

    point_coordinates = point["coordinates"]
    reference_point_coordinates = reference_point["coordinates"]
    lats = [point_coordinates[1], reference_point_coordinates[1]]
    lons = [point_coordinates[0], reference_point_coordinates[0]]
    reference_distance = geod.line_length(lons, lats)

    point_located_by_distance_away["coordinates"] = [point_coordinates[0] +
                                                     (reference_point_coordinates[0] -
                                                     point_coordinates[0]) * distance /
                                                     reference_distance, point_coordinates[1]]

    return point_located_by_distance_away

The method works for points located near the equator. However the further the points are located above or below the Equator, the bigger the error is. I do not understand why, since even if the ellipsis above equator is smaller, the reference distance should be enough to calculate the proportional distant point.


texts:
I'm currently working on QGIS 3.18 in the area of geo-referencing. I have two orthophotos (they are same) but the first one is untreated, the second was corrected by differential GPS. I want to measure the difference in geo-referencing between these images and particularly the shift.
My idea is to identify noticeable pixel in these images and measure the shift between these pixels, but I have no idea how to make that. How can I do this?


texts: Unable to access GeoServer from another computer in the same domain

texts:
I installed GeoServer 2.5.2, GeoNetwork 2.10.3 and Tomcat 7 on my Windows 7 computer {MyComputer} in our company's network. Now I would like to access GeoServer, GeoNetwork and Tomcat from other computers in the network.

However, I could only access via the web browser GeoNetwork and Tomcat but not GeoServer.

NO  -  http://{MyComputer}.{domain}.com:8070/geoserver/web/

YES -  http://{MyComputer}.{domain}.com:8040/geonetwork/srv/eng/main.home

YES -  http://{MyComputer}.{domain}.com:8080

Is there any configuration that I missed?


texts: I am struggling to convert .dwg to .shp, ideas?

texts:
Any ideas as to how can I convert dwg to shp?

I tried Kosmo and gvSIG. gvSIG kind of worked but it seems like I lost information..
in addition everything was converted into lines, no polygons. I can manually recreate the polygons, but there must be a better program to do this conversion.

Ideas?


texts: Setting CRS to ESRI:102001 using GeoPandas

texts:
I've been passed some GeoJSON shape data that I've been told is in the CRS ESRI:102001, although that property has not been written into the GeoJSON file. I've been able to confirm that is the correct CRS by converting it to EPSG:4326 using Python GDAL:

ds = gdal.VectorTranslate(input_file, srcDS=srcDS, srcSRS='ESRI:102001', dstSRS='EPSG:4326', format = 'GeoJSON', layerCreationOptions = ['RFC7946=YES', 'WRITE_BBOX=YES'])

I'd prefer to do the same in GeoPandas as I'm already using that for other tasks, while I'm only using GDAL for this one step.  I've had a few goes at setting the CRS with no luck so far.  This was my first effort:

inputGDF.crs = {'init' :'esri:102001'}
output_wgs84_GDF = inputGDF.to_crs("epsg:4326")

This fails on the 2nd line, with error:

pyproj.exceptions.CRSError: Invalid projection: +init=esri:102001 +type=crs: (Internal Proj Error: proj_create: cannot expand +init=esri:102001 +type=crs)

My second effort came from the proj4 info from:

https://spatialreference.org/ref/esri/102001/proj4/

inputGDF.crs = {'init' :'+proj=aea +lat_1=50 +lat_2=70 +lat_0=40 +lon_0=-96 +x_0=0 +y_0=0 +ellps=GRS80 +datum=NAD83 +units=m +no_defs'}
output_wgs84_GDF = inputGDF.to_crs('epsg:4326')

This throws a similar error:

pyproj.exceptions.CRSError: Invalid projection: +init=proj=aea +lat_1=50 +lat_2=70 +lat_0=40 +lon_0=-96 +x_0=0 +y_0=0 +ellps=GRS80 +datum=NAD83 +units=m +no_defs +type=crs: (Internal Proj Error: proj_create: cannot expand +init=proj=aea +lat_1=50 +lat_2=70 +lat_0=40 +lon_0=-96 +x_0=0 +y_0=0 +ellps=GRS80 +datum=NAD83 +units=m +no_defs +type=crs)

Am I missing something?


texts:
I have two datasets, the first is a polygon dataset of wellsite pads in NE British Columbia. The second is a point dataset representing well heads in the same area. The point dataset has an attribute describing the location of the well head, this is the attribute I want to use to label with.

The problem is, is that it is not a one to one relationship. One wellsite pad can have multiple wellheads on it and there is no attribute linking the two. What I want to do is label each well pad (polygon) with a list of stacked labels of well heads that are on that pad. See screenshot example of the kind of things I'm after.

I sketched up some examples in CAD to illustrate my problem. If I label the wellhead data I get something like the below sketch

How do I have the well head labels for each pad neatly ordered and stacked?

See screenshot below.


texts: Need to snap multiple polygons so no gaps

texts: Is there a tool to snap together multiple polygons so no gaps? Not polylines, but polygons. I have multiple polygons covering 8,100 sq. ft of territory for separate utility shape files; it took 3 days to create 15 polygons but it will take me a week to edit all the gaps out of area! Please tell me ArcMap 9.3.1 has a tool for this!


texts: Union of two Point Features which are close to each other, but have got an offset?

texts:
I have one Point Shapefile with the correct geographic position, lets call it A.

I have another Point Shapefile with the attributes I need, but with an offset of 0,7m or so (it varies a bit from point to point, but its very close to Point A).Lets call it B.

I would like to combine those two features automatically, so that I will have a Point Shapefile with the correct geographic position from A and the needed attributes from B. There's thousands of them so I cannot insert the attributes manually.
Is there a way to work with buffers or something?

(QGIS, etc.)


texts: I was wondering whether in QGIS you can break up a long text label attached to a polygon object into two lines and then move the label manually?  The label appears in small text as the polygon is small, so as well as this is there a way of increasing the font size of the new label?


texts: Adjacent DEM files have narrow triangular gaps when joined

texts:
Newbie here. I've downloaded several DEM files from GeoCommunity which correspond to adjacent regions. When I join them together to form a larger relief map there are triangular 'gaps' on all four sides of each 'tile'. It turns out this is because the elevations are all the same in these triangular regions. Visually, it looks like trimming the triangular regions would allow the remainder to fit properly.

Is this because the regions are not in fact rectangular? Something to do with UTM projection? If so how do you know which areas to trim?


texts: Rendering TopoJSON in Leaflet from GeoJSON PostGIS output?

texts:
I'm using the query found here SQL query to have a complete geojson feature from PostGIS to get a geojson that I can display on a leaflet map.

However, the amount of data is kinda big and the performance on Firefox is crappy (okay on Chrome), so I'm wondering if there is a way to convert the result of the query in Postgis into a topojson, in order to get a better performance


texts: Using ui.chart.image.series reducer with higher max pixel limit in Google Earth Engine?

texts:
Is it possible to use the ui.Chart.image.series function with an image collection, region, reducer, scale, AND a max pixel limit?

I'm dealing with 10m Sentinel-1 data, and even moderately sized regions hit the 1e7 default pixel count limit.

Anyway to override this in the time series plotting function?

I know how to do so using the "maxPixels:" argument in the image.ReduceRegion function.  Essentially I want something like this

print(ui.Chart.image.series(collection, polygon, ee.Reducer.count(maxPixels: 1e12), 10));


texts:
Using ArcGis 10.0. What I am looking for is a way to get point features that are placed in a polygon to have the same "Address" value.

This is the same as the "Spatial Join" feature, however, I am only trying to extract one field not all the fields in the polygon layer.

Is there a way to "Spatial Join" from a polygon to a point feature while selecting what fields from the polygon you wish to keep?


texts: Convert Pixel / array position to global coordinates Python

texts:
I am trying to convert positions in a raster representing a .tif to the corresponding global coordinates. Converting the whole array to a tif and load it to QGIS everything is referenced fine, but using the below calculation method for single points there is a slight offset (to east-north-east in the resulting coordinates....

raster.tif uses ETRS 89 UTM Zone 32N

Does anyone has an idea?

from osgeo import ogr, gdal, osr
import numpy as np

raster = gdal.Open("rasters/raster.tif")
raster_array = np.array(raster.ReadAsArray())

def pixel2coord(x, y):
     xoff, a, b, yoff, d, e = raster.GetGeoTransform()
     xp = a * x + b * y + xoff
     yp = d * x + e * y + yoff
     return(xp, yp)

print(pixel2coords(500,598))


texts:
I have an aplication, made with OpenLayers 3, GeoServer and ReactJS, showing some WMS. And for a strange reason I can't manage to refresh the image when I zoom in ONLY if the screen resolution is wider that 1314px. It's very strange, if the resolution is for example 1310, I can do normally zoom-in and zoom-out and the wms got refreshed after every zoom, BUT, if the resolution is wider than 1315px, ONLY zoom out works, but not zoom in.

Any one had experimented something like this before?

This is how I'm calling the wms (ReactJS):

    var map = this.state.map;
    var wmsLayer = new ol.layer.Image({
        source: new ol.source.ImageWMS({
            url: 'http://spif.XX.XX:8080/geoserver/cite/wms',
            params: {'LAYERS': 'cite:biodent_clientes'},
            serverType: 'geoserver'
        })
    });
    this.state.layersMaServer = wmsLayer;
    map.addLayer(wmsLayer);

There's no problem with GeoServer, the WMS are well loaded after every zoom-in or out, but I don't know why at zooming in isn't been refreshed on client side (broser). If I zoom-out just a little bit, then the image is refreshed with the correct wms, but why only under 1314px!

EDIT: The last logs of GeoServer:

https://paste.ofcode.org/32BP5Yvwi3JUsgugFnfPcdd

The first "event" (mouse scroll to zoom in) were made at 1314px wide, there's no problem on showing the wms. The last "event" was made at 1315px, there I can't make to OpenLayers to show the WMS, but they are correctly generated, I can see them on the Chrome Dev.


texts:
Can someone advise me where i might be going wrong when trying to open an Image in a Network folder using the Action function in QGIS
I have a point dataset with a column that has the image location on the Network. Network location is;
\gis\GIS\Environmental Health\Fly tipping\Images
A single record will therefore be \gis\GIS\Environmental Health\Fly tipping\Images\IMG_123.jpg
the action window is configured thus:

Using the Hotlink tool results in nothing opening
I've done this before with no problems and can also use a local drive location with no problems
I'm using QGIS 3.10 and Windows 10
thanks for any help in advance


texts:
Linux Mint 17.2.
After trying to reinstall QGIS through various processes referred to on SE, AskUbuntu & anywhere else*, I get various unmet dependency or held broken packages  errors.

Even with adjusting the repository details from saucy or xenial to trusty etc, I still can't get anywhere.

I have had this message a number of times:

The following packages have unmet dependencies.
 qgis : Depends: qgis-providers (= 2.14.3+dfsg-2~trusty1) but it is not going to be installed

I can't find out what this means to rectify it - anything I've followed up hasn't worked.

I am trying to install QGIS 2.18 with GRASS.

Can someone help?

*eg:
Broken packages and unmet dependency installation QGIS, PostgreSQL and postgis Ubuntu 14.04

QGIS install on Ubuntu 14.04 fails

https://askubuntu.com/questions/351899/qgis-quantum-gis-install-fails-unmet-dependencies


texts: How does i.fusion.brovey in Grass 6.4.3 work?

texts:
I'm using Grass 6.4.3 in Linuxmint 16.
I don't have a very big problem but I'm trying to understand what is Grass doing, if I'm applying the process correctly and if there is something wrong with my image files.

I'm trying to pansharpen a SPOT image. I have two files, one the panchromatic band and the other is the multiband file.

I import them to grass as, for example, 'spotimage.pan' (panchromatic) and the other file is automatic named by grass as 'spotimage.1', 'spotimage.2', 'spotimage.3' and 'spotimage.4'.

I think that in the original file, the bands order is inverted, I mean that the first band is NIR, the second one is Red and the third one is Green. This is because original image composite was done in order to obtain a false color image when is displayed as R,G,B = B1, B2, B3

So, if I do

d.rgb -o red=spotimage.1 green=spotimage.2 blue=spotimage.3

I obtain a false color view of the image.

Then, to pansharpen the image I use

i.fusion.brovey -s ms1=spotimage.3 ms2=spotimage.2 ms3=spotimage.1 pan=spotimage.pan outputprefix=spotimage.pansh

I inverted the order, where I should use the SPOT Band 1, I'm using the band 3 of my image and where I should use the SPOT Band 3, I'm using the band 2 of my image.

This because, in the dialog box of i.fusion.brovey, where it says "Name of input raster map (green: tm2 | qbird_green | spot1)" I suppose I have to use the name of the SPOT Band 1 and, where it says Name of input raster map (MIR; tm5 | qbird_red | spot3)" I suppose I have to use the name of the SPOT Band 3.

Is it correct?

The output are three raster files: pansh.red, pansh.green and pansh.blue.

But I cannot manage to obtain again a false color view with any combination of these pansharpened bands.

What am I doing wrong?

Are pansh.red the "pansharpened version" of spotimage.2, pansh.green the pansharpened version of spotimage.3 and pansh.blue the pansharpened version of spotimage.1? Remember thar the bands order is inverted.

Should I use i.landsat.rgb? Can I use it with a SPOT image?


texts: Stretched raster legend badly displayed ArcGIS Desktop?

texts:
I have a stretched raster in ArcGIS Desktop 10.3 with this symbology:

But when I put the layer in the legend, it's showed like this:

I've tried to change the styles, showing and hiding headers, changing size, text-size, etc. But it always seems broken. I've deleted the layer and put it again and nothing changes, I even have created the layer again.

What's going on?

Is not the first time I have this issue. Sometimes doing one of the mentioned things was solved, other times nothing worked.


texts: Using Clip on data from a GPKG results in empty result set - WARNING 000117

texts:
I'm trying to Clip the contents of a GPKG file (geodata.nationaalgeoregister.nl/brpgewaspercelen/extract/2020-concept/brpgewaspercelen.zip). I select the feature from a different layer as to reduce the output to the area of interest.
However every time I perform this opperation I'm met with:

WARNING 000117: Warning empty output generated.

I simply added the GPKG file from the add data button. The other layer has no problems, as I used that a milion times before when doing operations like this, but then on a normal layer file.


texts: How can I get the circle of a sphere (on Earth ellipsoid) around a point using the Latitude and Longitude?

texts: Specifically finding points along the circle. If possible without converting into Northing/Easting beforehand. I'm looking for the fastest possible way to calculate points in Latitude/Longitude in a circle around a point on the Earth ellipsoid. I've seen ways to figure out a circle of a sphere, however not while using Latitude and Longitude angles directly and not on an ellipsoid.


texts: Extracting certain area from DEM using ArcGIS Desktop?

texts: Using ArcGIS Desktop, how do I extract a certain area (uniform 1 km from the entire coastline) from a particular year's DEM image?


texts: Setting output format on MousePosition control

texts:
Having a bit of trouble getting the format to change on the MousePosition control. I want to change it to {y}N, {x}E but can't quite work it out.

var mousePositionControl = new ol.control.MousePosition({
    coordinateFormat: ol.coordinate.createStringXY(4),
    projection: 'EPSG:4326',
    // comment the following two lines to have the mouse position
    // be placed within the map.
    className: 'custom-mouse-position',
    target: document.getElementById('mouse-position'),
    undefinedHTML: ''
});


texts: Calculating average raster values around set of points using ArcGIS for Desktop (or R)?

texts:
I have a 30-m elevation raster as well as a set of points that overlay the raster.

For each point, I want to calculate the average elevation within a 5-km radius of the point.

One way I was thinking of doing this was by creating 5-km buffers around each point, then using those buffers to calculate the average values of the elevation raster.  This should work but I find it hard to believe that there is no easier way of doing this.  Does anybody know?

I would prefer a solution using ArcGIS 10.1 for Desktop but I could use R too.


texts:
I have a field called "Values" in ArcMap where I am just wanting to Display the 10% of the values on the map and the remainder of the values not be shown on the map.

Is there a way I can do this?
I was thinking perhaps I could do a selection? Select the top 10% values in the field?


texts: Reading geotiff as array gives Zero matrix or NoDataValue (-100000) in Gdal?

texts:
from osgeo import gdal
import numpy as np
red = gdal.Open("/home/ubuntu/Agriculture/data/sensefly/indices_sens_op/agri_sense_index_nir.tif")
red_array = red.ReadAsArray().astype(np.float32)
print (red_array)

out [[[0. 0. 0. ... 0. 0. 0.]
  [0. 0. 0. ... 0. 0. 0.]
  [0. 0. 0. ... 0. 0. 0.]
  ...
  [0. 0. 0. ... 0. 0. 0.]
  [0. 0. 0. ... 0. 0. 0.]
  [0. 0. 0. ... 0. 0. 0.]]

 [[0. 0. 0. ... 0. 0. 0.]
  [0. 0. 0. ... 0. 0. 0.]
  [0. 0. 0. ... 0. 0. 0.]
  ...
  [0. 0. 0. ... 0. 0. 0.]
  [0. 0. 0. ... 0. 0. 0.]
  [0. 0. 0. ... 0. 0. 0.]]]

the ortho is from Pix4d mapper and when called using numpy it gives only zero matrix or gives nodata value (-10000) in the matrix.


texts:
This query

SELECT
     a.id,
     b.id AS id2,
     ST_MakeLine(a.geom, b.geom)
 FROM riverscale.points a CROSS JOIN riverscale.points b
WHERE a.id != b.id;

generates this:

I am guessing the necessary information for creating a geomtery is in the column "st_makeline".

How can I convert this table into an object with a LINESTRING (or other) geometry that I can load in e.g. QGIS?

#

EDIT:
When I try to create a table with this query I get this and I can't load it in QGIS.

CREATE TABLE riverscale.test as
   SELECT
     a.id,
     b.id AS id2,
     ST_MakeLine(a.geom, b.geom)
 FROM riverscale.points a CROSS JOIN riverscale.points b
WHERE a.id != b.id;


texts:
I used Random Forest regression for my work, but now I have another problem because I have really big data. I am using a VHR and this is much raster information. My script in R is this:

rand.zero <- apply(rand.riq, 1, function(y) sum(length(which(y > 0))))
hist(rand.zero, breaks=50, col="grey", main="",
     ylab="Number of OTUs", xlab="Number of Non-Zero Values")

Remove zero data

remove_rare <- function( table , cutoff_pro ) {
  row2keep <- c()
  cutoff <- ceiling( cutoff_pro * ncol(table) )
  for ( i in 1:nrow(table) ) {
    row_nonzero <- length( which( table[ i , ]  > 0 ) )
    if ( row_nonzero > cutoff ) {
      row2keep <- c( row2keep , i)
    }
  }
  return( table [ row2keep , , drop=F ])
}

Remove non-zeros values

rand_table_rare_removed <- remove_rare(table=rand.riq, cutoff_pro=0.5)
dim(rand_table_rare_removed)

rand_table_rare_removed_norm <- sweep(rand_table_rare_removed, 2, colSums(rand_table_rare_removed) , '/')*100

Scale or data transformation

rand_table_scaled <- scale(rand_table_rare_removed_norm, center = TRUE, scale = TRUE)
dim(rand_table_scaled)
x = nearZeroVar(rand_table_rare_removed)

Regression

rand_table_scaled_IS <- data.frame(t(rand_table_scaled))
rand_table_scaled_IS$riq <- rand.riq[rownames(rand_table_scaled_IS), "riq"]
set.seed(45)
RF_IS_regress <- randomForest( x=rand_table_scaled_IS[,1:(ncol(rand_table_scaled_IS)-1)] , y=rand_table_scaled_IS[ , ncol(rand_table_scaled_IS)] , ntree=500, importance=TRUE, proximities=TRUE )
print(RF_IS_regress)

And finally the program shows me this error:

> Error in randomForest.default(x = rand_table_scaled_IS[,
> 1:(ncol(rand_table_scaled_IS) -  :    NA not permitted in predictors

But I removed zero data and don't exist Nas


texts: st_intersection only if specific conditions are satisfied

texts:
I'm new to sf package.
I have two sf files, one containing lat/long position of points (point_sf) and the other one containing sf polygon geometries (polygon_sf). I want to see which points belong to specific polygons (using sf_intersection this is easy) however certain conditions must be met (I believe this would also reduce multiple options).
For example, I want to "match" only points that are within specific altitude band, which is defined in one of the columns in polygon_sf, and that years (values exist in both files) are the same. I have tried the following but it gives NULL output.
point_new <- if(point_sf$year== point_sf$year&
                  point_sf$FL >= polygon_sf$from_FL &
                  point_sf$FL <= polygon_sf$to_FL)  {
  st_intersection(point_sf, st_difference(point_sf))}

Any suggestions?


texts: Sync-ing ArcGIS Online and collector for ArcGIS?

texts:
enter image description hereI am new with "collector for ArcGIS" and we have just recieved the ArcGIS Online account to be able to provide data for Collector.
Based on the General instruction, we have shared some data ( Shape, Layerpacket,...)in Arcmap and were able to see them in ArcGIS online.

We shared them through our organization in ArcGIS online. Though the data editable is, still we can not see it in Collector.

I collector we have two Option either my Maps or all Maps and still non of them is visible.

What could be the Problem?


texts: Buffer from Point or Polygon with Information from Raster

texts:
The intent of this  work is to define a buffer around a point or polygon based on information from a raster. For instance if a polygon representing a city and the raster defines how long it takes to pass through that cell (in hours for simplicities sake). I would like to construct a buffer around this city defining the maximum range that would sum up to 24 hours. Any cell can move to any other cell so long as they are moving towards the origin (The value of the origin cell is not included in the sum).

This will allow for the creation of non-uniform buffers based on environmental information rather than setting an arbitrary buffer around a polygon.

set.seed(8354)
r <-raster(matrix(ncol=5,nrow=5,abs(round(rnorm(25,5,5),0))))
e <- extent(c(0,5,0,5))
extent(r) <- e
plot(r)
points(2.5,2.5,cex=3,pch=19)

Below the images show the initial raster and the raster "buffered" out to outward values that equal 10. The intent to to move outward from the center cell in all 8 directions summing up values as progress is made. Upon achieving 10 the buffer is set to that extent, if the value exceeds 10, the current cell limits the extent.

One challenge I can already see is counting diagonally, in one instance the value may be less than 10 and in another greater, not really sure how to get around that as of yet.

I have been looking through the gdistance documentation as well as accCost as seen here. The transition function seems to be the hangup point. Much of that stems from my lack of understanding as to what this function is doing (I'll continue my reading). I've been using sum in the transitionFunction argument.

Values on the original raster

Outward summation defining buffer extent


texts: PyQGIS 3.0: cannot make it work under Windows

texts:
I've just tried to install QGIS 3.0 package in my Anaconda (Win7x64, Python 3.6 x64). I wanted to use QGIS widgets in my PyQt application. Alas, QGIS works incorrectly: at first, when I just tried to import qgis modules in Python, it said that it cannot find them. And when I manually set the PYTHONPATH, QGISPATH and PATH variables, it simply does not show my map (only the blank white screen).

Under Ubuntu 18.04 the script works correctly (there I just made apt-get install python3-qgis).

Can you tell me, how can I make it work under Windows too?

Here's my test script:

#!/usr/bin/env python3

import os

from PyQt5.QtWidgets import QWidget, QFrame, QHBoxLayout, QGridLayout

from qgis._core import QgsApplication, QgsRasterLayer
from qgis._gui import QgsMapCanvas

class TestWidget(QWidget):
    def __init__(self, parent=None):
        super(TestWidget, self).__init__(parent)
        self.main_layout = QHBoxLayout(self)
        self.canvas = QgsMapCanvas()
        self.main_layout.addWidget(self.canvas)

        base_layer_path = os.path.join(os.path.abspath(os.path.dirname(__file__)), "world.mbtiles")
        self.base_layer = QgsRasterLayer(base_layer_path, 'world')
        self.canvas.setLayers([self.base_layer])
        self.canvas.zoomToFullExtent()

if __name__ == '__main__':
    qgs = QgsApplication([], False)
    qgs.setPrefixPath(os.environ.get('QGISPATH', '/usr'), True)
    qgs.initQgis()

    w = TestWidget()
    w.show()

    qgs.exec_()
    qgs.exitQgis()

And here's my test map data (world.mbtiles): https://www.sendspace.com/filegroup/8OqffqxS3aWuRdwr4iolRQ

And here are my environment variables (set under Windows):

PYTHONPATH=C:/Anaconda3/Library/python
QGISPATH=C:/Anaconda3
PATH=C:/Anaconda3/share/qgis/python/qgis;C:/Anaconda3/Library/bin


texts: How to subset a SpatialPoints object to get the points located on each side of a SpatialLines object using R?


texts:
I am a programmer (not a cartographer!!) developing an application which needs to read and write GeoTIFF files. I have a number of questions about those GeoTIFF keys which are related to projections.

I would ideally like to be able to read/write files which originate from and/or will be read by a maximum number of other GeoTIFF-aware software packages. My primary compatibility target is ARC, however, so if there is anything especially peculiar about files written or read by ARC I would be particularly interested in those details. Also, I need to support only a relatively small number of projections, namely Mercator, UTM, polar stereographic and Lambert conic conformal in that order of importance. I am using a home-brewed projection package based closely on Snyder's USGS reference and code from the GMT package (as we need to maintain co-registration with GMT plots.) I have looked through the GeoTIFF spec at www.remotesensing.org and other related materials in some detail but am still confused on a number of points:

(1) The GeoTIFF spec declares various lon/lat keys, almost all of them with deliciously circular definitions, e.g.:

"ProjNatOriginLongGeoKey: Longitude of map-projection natural origin."

Can someone explain to me the difference between natural origin, false origin and center lon/lat? My limited understanding of the math behind Snyder's projections includes knowledge of the fact that virtually all projections require a 'projection longitude' parameter, for instance. With regard to Mercator, this is the longitude at which all projected x-values will be 0 (likewise for UTM, before the 500000 false easting is applied). So, is this longitude the ProjNatOriginLong or the ProjCenterLong in GeoTIFF jargon?

(2) If a GeoTIFF file contains no specified value for a key, e.g., ProjFalseEasting, is there an assumed default value for that key? If so are such values projection-dependent, e.g. default false easting might be 500000 for UTM but 0 for anything else? I have looked at a couple of GeoTIFF files generated by external software and there seems to be a wide variety of quantity of keys specified. The spec says nothing about default values but I suspect that various external readers/writers make many assumptions about these keys (which no doubt differ on a vendor-by-vendor basis).

(3) As I noted I am primarily interested in reading ARC files. A sample Mercator GeoTIFF from ARC had a defined ProjectedCSTypeGeoKey of 3395, not listed anywhere in the GeoTIFF spec. This threw me for a while until I started digging around in the ancillary files distributed with libgeotiff and found the line

3395,"WGS 84 / World Mercator",9001,4326,19883,9804,1,0, 4400,8801,0,9102, 8802,0,9102,8805,1,9201,8806,0,9001,8807,0,9001,,,,,,

buried in /usr/share/epsg_csv/pcs.csv. Is this file pretty much invariant, meaning that I can count on this entry never changing, or do I really need to write code which parses these interconnected CSV databases everytime I encounter an unrecognized ProjectedCSType in an imported file? Like I said, I only need to support a handful of the most common projections.

Thanks very much for any help in clarifying these points!

Roger Davis
Univ. of Hawaii


texts: Setting zoom level for labels of highway in JOSM?

texts:
I'm new in using JOSM and OpenStreetMap data . I import a shapefile of highways in JOSM and use it for styling the highways.Some highways are trunk ( in open street map standard ) and some of are primary . When i style highways, the name of highways  display in zoom level of  17 on map . I want to labeling the highways in zoom level of 12 .

How can i change the zoom level for name of the highways?


texts: Multiple import of xml metadata catalogues into GeoNode

texts:
In our GeoNode's database (Postgis), we have imported many shapefiles using the tool shp2pgsql and then the command manage.py updatelayers. After that, we created their metadata catalogues from the INSPIRE editor:
http://inspire-geoportal.ec.europa.eu/editor/

Is it possible to import the xml metadata catalogues in GeoNode db without importing again the data ?


texts:
I have created a model with ModelBuilder to interpolate values with the IDW tool, iterating through rows.

How can I restrict interpolated values to regions within a 20km buffer?

I tried setting the "search radius" to 20km, but that only prioritizes points within the radius - the radius increases until the minimum number of points is encompassed. Instead, I want anything outside of the search radius to be set as missing or no data.


texts: Loop and perform loop in order depending on attribute

texts:
I have a lot of points in one layer. All of which represent the same type of real world feature.

I want to perform a delete identical within XY of each point but in a specific order. (by splitting them into a layer each if there is no smarter way to do this within the layer??)

I want the order to start with the point closest to a specific line feature so I have added a distance attribute to each point with the near tool.

Thus the first time the operation runs it deletes identical points close to the point with the distance attribute closest to 0.

The second time it should then move on to the next point but it will skip those points deleted because they were too close to the first point. and repeat until all points have been dealt with.

How can I make such a loop either dealing with just the one layer or with having each point split into its own layer (might be easier?)

The pic below illustrates the line and blue points will be kept while red ones deleted, as a type of generalization. The red circle represents the XY tolerance.

I want to try and make this in FME or ArcMap ModelBuilder.


texts:
I used r.horizon to calculate horizon on a raster every 10° without any offset. I have got 36 files now.

Is the first file from 0° to 10° (then centered on azimuth=5°) or is it centered on 0°, in this case the first file is from 355° to 5° (knowing that GRASS start 0 in the East counterclockwise)?


texts:
I am working with OpenLayers (and GeoServer) and I am showing a WMS layer with countries in a map.
Each object (country) in the layer has an attribute called NAME (the name of the country).

I created a filter like the following in order to see only Italy (it's an example)

myFilter = new OpenLayers.Filter.Comparison({
    type: OpenLayers.Filter.Comparison.EQUAL_TO,
    matchCase: false,
    property: 'NAME',
    value: 'Italy'
});

This work perfectly and if I apply it to the layer and call redraw() method I can see only one country.

Now I would like to do the same thing with sytles, in order to highlight Italy but keep visible all the other countries in the map.
I created a style like the one below, but now I don't know how I can apply this to the layer.

var style = new OpenLayers.Style();
var ruleLow = new OpenLayers.Rule({
    filter: myFilter,
    symbolizer: {fillColor: "green"}
});
style.addRules([ruleLow]);

Can I apply the style to a WMS layer? Or only to vector layer?

Thank you


texts: Start Editing Return Unspecified error (Exception from HRESULT: 0x80004005 (E_FAIL))

texts:
I'm using IEngineEditor.StartEditing method to start an edit session on SDE Workspace. My workspace is having a right value because I have no problem to get a feature class using it. I'm not sure where I did wrong, or what causing this error. Please help.

 IPropertySet propertySet = new PropertySetClass();
 propertySet.SetProperty("SERVER", MyServer);
 propertySet.SetProperty("INSTANCE", MyInstance);
 propertySet.SetProperty("DATABASE", MyDatabase);
 propertySet.SetProperty("USER", MyUsername);
 propertySet.SetProperty("PASSWORD", MyPassword);
 propertySet.SetProperty("VERSION", MyVersion);

 Type factoryType = Type.GetTypeFromProgID("esriDataSourcesGDB.SdeWorkspaceFactory");
 IWorkspaceFactory workspaceFactory = (IWorkspaceFactory)Activator.CreateInstance(factoryType);

 IWorkspace workspace = workspaceFactory.Open(propertySet, 0);

 IFeatureWorkspace pFeatureWorkspace = workSpace as IFeatureWorkspace;
 IFeatureClass pFeatureClass = pFeatureWorkspace.OpenFeatureClass("MyFeatureClassName"); //I can get this feature class. This was for later use
 IMap pMap = mapControl.Map;

 engineEditor.StartEditing(workSpace, pMap);


texts: How to customize the QGIS GUI using Python?

texts:
Is it possible to hide/remove/customize the default UI of QGIS from a Python plugin?

I want to add my own menus and toolbars and remove some included by default in QGIS.

If it's possible, which reference should I use?


texts: Creating table with data from raster to samples

texts:
How do I obtain a table with the data of a raster for my single points in the map?
The image might be helpful.  What I want to do is to have a table with the value from the 'Elevation' raster to my '1001 genomes' samples.


texts: How to calculate the area of several classes of a raster in a polygon of a vector file with QGIS?

texts:
I have some trouble with QGIS

I have two layers of information:

The first one a vector that contains 500 objects, with attributes such as number of population, and...
A raster that has three classes: degrade (-1), stable (0) and improvement (1), plus -9999 which corresponds to non-values

What I would like to obtain is in the attribute table of my vector layer, the percentage (area) of the three classes of my raster

something like this, in the row one :
- the id, the population, the area of class 1, the area of class 2, the area of class 3.
and so on for my other 500 lines.

I tried thing with Zonal statistics but I couldn't obtain what I wish.


texts: How can I use Arcpy to set label expression parser to Python?

texts:
I looked and looked but I can't find anything that lets the user to change the parser default to python for label expressions. The default is vb so when I run my arcpy script it list my python expression in the vb which doesn't work because it's in a different language. Can someone help me?

import os
import arcpy # if outside ArcMap

# inside ArcMap, use "CURRENT" for document name
mxd = arcpy.mapping.MapDocument(r"C:\Users\aa2zz6\Desktop\test\Untitled.mxd")
df = arcpy.mapping.ListDataFrames(mxd, '')[0]

for lyr in arcpy.mapping.ListLayers(mxd, "P_Meters"):
    print lyr.name
    if lyr.supports("LABELCLASSES"):
        print lyr.name + " supports label classes"
        for lblClass in lyr.labelClasses:
            print lblClass.className
            lblClass.expression = lblClass.expression = '"{}" + [HOUSENUMBER] +  "{}"'.format("<CLR red = '255' green = '255' blue = '255'><FNT size = '10'>","</FNT></CLR>")
            if lblClass.showClassLabels:
                print "    Class Name:  " + lblClass.className
                print "    Expression:  " + lblClass.expression

lyr.showLabels = True
del mxd


texts:
How can I add data based on neighbouring values?

If have a shapefile, dividing the country into different areas.
I have population data for each area.
I would like to add new fields to my data, detailing the population of the neighbouring areas, for any area that shares a land border.

For example, in the case of area A, I would have the population of B and C in the adjacent fields (i.e. neighbour1 is B and neighbour2 is 3).


texts: How to compile QGIS from source with Grass?

texts:
I recently decided to reinstall QGIS and build against the latest version of GDAL - something I have been putting off for ages because I suspected it would be a nightmare. Lo and behold it has been, so after a full-day of trying to get this to compile following these: instructions I'm reaching out for some help.

I successfully built GDAL 2.0 this morning (at least I think it was successful gdalinfo --version is now showing GDAL 2.0.0, released 2015/06/14 and have been working my way through the instructions above.

I'm currently stuck at 47% with make failing here:

[ 47%] Built target qgis.d.rast7
Linking C executable ../../../../output/lib/qgis/grass/modules/qgis.g.info7
/usr/lib/grass70/lib/libgrass_gproj.so: undefined reference to `OSRGetLinearUnits@GDAL_1.8'
/usr/lib/grass70/lib/libgrass_gproj.so: undefined reference to `OSRExportToWkt@GDAL_1.8'
/usr/lib/grass70/lib/libgrass_gproj.so: undefined reference to `OSRGetAttrValue@GDAL_1.8'
/usr/lib/grass70/lib/libgrass_gproj.so: undefined reference to `SetCSVFilenameHook@GDAL_1.8'
/usr/lib/grass70/lib/libgrass_gproj.so: undefined reference to `OSRExportToProj4@GDAL_1.8'
/usr/lib/grass70/lib/libgrass_gproj.so: undefined reference to `OSRMorphFromESRI@GDAL_1.8'
/usr/lib/grass70/lib/libgrass_gproj.so: undefined reference to `OSRExportToPrettyWkt@GDAL_1.8'
/usr/lib/grass70/lib/libgrass_gproj.so: undefined reference to `OSRMorphToESRI@GDAL_1.8'
/usr/lib/grass70/lib/libgrass_gproj.so: undefined reference to `OSRImportFromProj4@GDAL_1.8'
/usr/lib/grass70/lib/libgrass_gproj.so: undefined reference to `OSRDestroySpatialReference@GDAL_1.8'
/usr/lib/grass70/lib/libgrass_gproj.so: undefined reference to `VSIFree@GDAL_1.8'
/usr/lib/grass70/lib/libgrass_gproj.so: undefined reference to `OSRIsProjected@GDAL_1.8'
/usr/lib/grass70/lib/libgrass_gproj.so: undefined reference to `OSRNewSpatialReference@GDAL_1.8'
/usr/lib/grass70/lib/libgrass_gproj.so: undefined reference to `OSRIsGeographic@GDAL_1.8'
/usr/lib/grass70/lib/libgrass_gproj.so: undefined reference to `OSRGetUTMZone@GDAL_1.8'
collect2: error: ld returned 1 exit status
make[2]: *** [output/lib/qgis/grass/modules/qgis.g.info7] Error 1
make[1]: *** [src/providers/grass/7/CMakeFiles/qgis.g.info7.dir/all] Error 2
make: *** [all] Error 2

I'm not quite sure what this means - though some searching suggests I may still have an old install of GDAL lurking, though not sure I can find it - as I removed a couple earlier when make failed much earlier on.

I thought it could be something to do with Grass 7 which I installed following these enter link description here instructions. But I've also tried to compile with GRASS_PREFIX7 set to OFF in ccmake. But make fails at the same place.

p.s. using:
Ubuntu 14.04
GRASS7
GDAL 2.0


texts: Are landsat 8 SR and TOA rasters calibrated/corrected data or the data to be used for calibration/correction?

texts:
I downloaded Landsat 8 Level-2 scenes from https://espa.cr.usgs.gov/ for land cover classification. Following are the product list in the package, cause I ticked Top of Atmosphere Reflectance (TOA), Surface Reflectance (SR) and Pixel QA to be processed.

LC08_L1TP_129026_20140712_20170421_01_T1 = *

*_b1 to b11
*_sr_band1 to band7
*_toa_band1 to band9
*_bqa
*_MTL
*_ANG
*_pixel_qa
*_radsat_qa
*_sr_aerosol
*_sensor_azimuth_band4
*_sensor_zenith_band4
*_solar_azimuth_band4
*_solar_zenith_band4

I am new to image calibration and correction things and I thought the image quality will be significantly increased after these processing. My understanding was that the calibration and corrections will be done on same raster band of same scene, For example all of the SR and TOA and all other angle corrections can be processed on same raster layer of Band 1 of LC08_L1TP_129026_20140712_20170421_01_T1 scene all together. But what I got separate files from the ESPA web, such as apart from original band data (b1-b11), I got surface reflectance (b1-b7), top of atmosphere reflectance (b1-b9) and several other azimuth and sensor related files.

Now I am confused, are they corrected/calibrated data of a certain band (which one is corrected and calibrated if they are) or are they data to be used for calibration/correction process of the original band data?


texts: Snapping multiple points to nearest edge of polygon using ArcGIS for Desktop?

texts:
I have a point grid (representing the centroid of a map tile) laid across a landbase; however, some of those points fall outside the boundaries (in the ocean).  I need to somehow snap those points back inside the polygon (landbase) so I can use an identity tool and return results for all of the points.  Is there a simple way to do this? (There are hundreds of points)

There may be a better way to do this too... I'm open to suggestions.

I'm trying to use identity to combine two datasets (map tiles and counties), but I want to make sure each point (map tile) only falls within one of the other polygons (counties).  I hope this is at least somewhat clear...


texts: How to restrict Mapshaper's `-dissolve` to merge only certain specified features on a layer

texts:
I'm try to merge a few specified features on a topojson map using Mapshaper. For this example, I'm trying to merge Somalia, Somaliland and Puntland to match a dataset that doesn't separate them.

Mapshaper has command features, -dissolve and -dissolve2, that seem to do what I want. However, they seem to apply to every feature in the selected layer, and I can't see any obvious way in the Mapshaper docs to select only certain features within a layer to apply to.

The options name and target apply to files and layers respectively, not features. Likewise all the usage examples I can find don't select specific features to merge.

There's a feature to limit it by field:

   Name of a data field to dissolve on.

... but I can't find any way to pipe in a list of field values. I can specify the name field, but can't find a way to specify that the name field should be "Somalia" or "Somaliland".

If possible, I'd prefer not to move the features onto a temporary dummy layer then move them back again.

Ideally I'd use something similar to how -each and -filter allow to match by an expression but that doesn't seem to exist for -dissolve. For example, this, which works perfectly for -filter, doesn't work:

 -dissolve '(NAME == "Somaliland" || NAME == "Puntland" || NAME == "Somalia")'

...because that parameter must be a field name and can't be an expression:

  Data table is missing field: (NAME == "Somaliland" || NAME == "Puntland"  || NAME == "Somalia")


texts: Dissolving a TIN using python/ogr binding

texts:
I have a TIN model where I remove some of the triangles (see fig. 1). Most of the triangles touch on one edge, some on one vertice. I would now like to generate a polygon shapefile containing the remaining triangles, dissolved into one polygon (of multiple rings, see fig. 2).

Fig 1. - TIN with some triangles removed. Vertices of touching triangles are exactly identical, no sliver polygons.

Fig. 2 - Note, that this is only one polygon with a single line in the attribute table

Pictures were made by using the QGIS "Vector-->Geoprocessing Tools-->Dissolve" tool.
Is this possible with the ogr python bindings?
Here's what I've tried so far:

from osgeo import ogr

tri = ogr.Geometry(ogr.wkbLinearRing)
tri.AddPoint(0,0)
tri.AddPoint(0,1)
tri.AddPoint(1,1)
tri.CloseRings()
poly = ogr.Geometry(ogr.wkbPolygon)
poly.AddGeometry(tri)

tri2 = ogr.Geometry(ogr.wkbLinearRing)
tri2.AddPoint(1,1)
tri2.AddPoint(0,1)
tri2.AddPoint(0,4)
poly2 = ogr.Geometry(ogr.wkbPolygon)
poly2.AddGeometry(tri2)

poly.Union(poly2) # returns None
poly.Intersect(poly2) # returns False
poly.Touches(poly2) # returns False (?)
poly.Disjoint(poly2) # returns False

I'm missing something like poly.Dissolve(), maybe used on a ogr.wkbGeometryCollection?


texts:
I used to use spyder as python IDE and was learning QGIS developing recently. But I cannot import qgis.core in spyder although the library can be imported in QGIS Desktop Console successfully.

Here is my setting environment:

My OS is Windows 7;

Spyder was integrated in Anaconda2-4.1.1-Windows-x86_64 and install together, in which default interpreter is Python 2.7.12(64bit)

QGIS 2.18 was installed by Standalone installer 64bit but python is 2.7.5 64bit.
Then I modified the environment variables by appending

“C:\Program Files\QGIS 2.18\apps\qgis\python” in PYTHONPATH and

“C:\Program Files\QGIS 2.18\apps\qgis\bin;

C:\Program Files\QGIS 2.18\apps\qgis\python;

C:\Program Files\QGIS 2.18\apps\qgis\python\qgis\core;” in Path.

However, I still cannot import it in Spyder’s IPython, error says”DLL load failed:”

Same error even I change the interpreter to QGIS’s 2.7.5 with Desktop open.

In term of OSGeo4W Shell, that is another problem.

Still failed to import it when finding no module at first opening, but it works only once if run Germán Carrillo’s BAT file.

Which means that I have to run this batch file every time I start OSGeo4W Shell.

So it is weird my addition path in PATH and PYTHONPATH have no effect at all.


texts: Showing WKT multipolygon in Mapserver results msGEOSError?

texts:
I have multipolygon WKT and I want to show it as a layer in Mapserver, I followed Mapserver as in this link http://mapserver.org/mapfile/feature.html but I can see it didn't mention MULTIPOLYGON as supported feature and if I try to add MULTIPOLYGON WKT layer, I am getting below error:

[MapServer Error]: msGEOSError():0'ƒ╖▀_tΦ≡^☺ÇΘ^☺∞µ^☺Φ≡^☺0v`♦⌠≡^☺↑τ^☺♦▀_t0v`♦Φ≡^☺

Is there any other way ?


texts:
I am using GDAL to clip my DEM by the selected polygon in the active layer with QGIS 3.4. I am facing a problem with the memory layer. I read the questions asked before related to this topic PyQGIS Processing: "memory: not found", Loading temporary layer in processing PyQgis 3.4
However, the solutions that are proposed are not working for me.
With the following code, I do not get any result as the memory layer is not created

res1 = processing.run("gdal:cliprasterbymasklayer", {
    'INPUT':demPath,
    'MASK':QgsProcessingFeatureSourceDefinition(areaPath, True),
    'NODATA':99999,
    'ALPHA_BAND':False,
    'CROP_TO_CUTLINE':True,
    'KEEP_RESOLUTION':True,
    'OPTIONS':'',
    'DATA_TYPE':6, # Float32
    'OUTPUT':'memory:test'})
QgsProject.instance().addMapLayer(res1['OUTPUT'])

So my workaround to solve the problem is the following.

##Clip DEM
res1 = processing.run("gdal:cliprasterbymasklayer", {
    'INPUT':demPath,
    'MASK':QgsProcessingFeatureSourceDefinition(areaPath, True),
    'NODATA':99999,
    'ALPHA_BAND':False,
    'CROP_TO_CUTLINE':True,
    'KEEP_RESOLUTION':True,
    'OPTIONS':'',
    'DATA_TYPE':6, # Float32
    'OUTPUT':QgsProcessingUtils.tempFolder()+'/OUTPUT' + str(random.getrandbits(128))+'.tif'})
rlyr_demClp = QgsRasterLayer (res1['OUTPUT'], "rlyr_demClp", "gdal")

However, I would like to ask whether there is a way to use the memory layer in QGIS 3.4. I am using the same processing tool in QGIS 3.10 without any problem as follows.

# Clip DEM in QGIS 3.10
res3 = processing.run("gdal:cliprasterbymasklayer", {
    'INPUT': rlyr_dem,
    'MASK': vlyr_ezgStats,
    'NODATA': 99999,
    'CROP_TO_CUTLINE': True,
    'DATA_TYPE': 6,  # Float32
    'OUTPUT': 'TEMPORARY_OUTPUT'})


texts: Drawing on jpg-map image with proj4 (perl)

texts:
From the point of view GIS professionals this question will be sure very primitive, but please keep in mind - i'm totally newbie in this area - and this is my 1st question here.

Have several GPS coordinates and need draw them as pins the into world-map image. On the internet i found one "public domain" image of the world. (attached a resized smaller version here - in reality I need 2000 pixel wide image).

I know perl - so, decided to use:

GD.pm (or Image::Magick) - for graphics operations
and installed the Geo::Proj4 perl module to cartographic conversions.

The manual says:

  $proj->forward(LATITUDE, LONGITUDE)

  Perform a forward projection from LATITUDE and LONGITUDE (LL) to the
  cartographic projection (XY) represented by the Geo::Proj4 instance

  ...

  Returned are two values, usually X and Y in meters, or whatever
  units are relevant to the given projection

If I understand right - that is enough for me, need only convert the meters to my pixels in my map image.

For using the above forward routine, it need an $proj instance what must define the map projection TYPE, like:

$proj = Geo::Proj4->new(proj => "merc", ellps => "clrk66", lon_0 => -96);
#or
my $proj = Geo::Proj4->new(init => "epsg:28992");

#and like (many-many-many possibilities).

As the manual says - It is not easy. :(

  The object defines the target projection, but that's easier said than
  done: projections have different parameter needs. The parameters which
  can (or need to) be used are listed with cs2cs -lP. The manual page of
  cs2cs explains how the configuration works.

Now the questions:

In what projection is the attached map? - or how to determine it?
How i can decide how to initialize the Proj4 instance for this image?
I'm totally lost how to continue - the cs2cs -lP. returns me approx 150 different projections :(
if this image is not usable - is here some better way to get an better defined map? Because need generate many-many-many bitmap images, can't use google API or like)

Thank you for a patience.