texts: Remove cloud cover and calculate EVI in Sentinel-2 2a by GEE

texts:
I am trying to calculate EVI by GEE for every 15 days for two years. I have put the code to remove the cloud mask and wondering if my code is right?
// Lists of Years and Months
var startDate = '2018-12-31';
var endDate = '2020-12-31';

var Start = ee.Date('2018-12-31');
var End = ee.Date('2020-12-31');
var Months = ee.List.sequence(1,12);
var Years = ee.List.sequence(ee.Number(Start.get('year')),ee.Number(End.get('year')));

var i;
// 49
for (i = 0; i < 49; i++) {
  // print(i);
  print("---------------------------------");
  var step = 15 * i;
  var newDate = ee.Date(startDate).advance(step, 'days');
  print(newDate);

var for_export = ee.ImageCollection('COPERNICUS/S2_SR')
                  .filterDate(newDate, newDate.advance(24,'hour'))
                  // Pre-filter to get less cloudy granules.
                  .filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE',30))
                  .filterBounds(AOI)
                  .map(maskS2clouds);
}

function maskS2clouds(image) {
  var qa = image.select('QA60');

  // Function to mask cloud from built-in quality band
// information on cloud
var maskcloud1 = function(image) {
  var QA60 = image.select(['QA60']);
  var clouds = QA60.bitwiseAnd(1<<10)
           .or(QA60.bitwiseAnd(1<<11));
  return image.updateMask(clouds.not());
};

  // Bits 10 and 11 are clouds and cirrus, respectively.
  var cloudBitMask = 1 << 10;
  var cirrusBitMask = 1 << 11;
  // Both flags should be set to zero, indicating clear conditions.
  var mask = qa.bitwiseAnd(cloudBitMask).eq(0)
      .and(qa.bitwiseAnd(cirrusBitMask).eq(0));

  return image.updateMask(mask).divide(10000);
}

var s2_QA = for_export.map(maskS2clouds)
print(s2_QA);

// Add EVI using an expression.
var addEVI=function(image){
var evi = image.expression(
      '2.5 * ((NIR - RED) / (NIR + 6 * RED - 7.5 * BLUE + 1))', {
      'NIR' : image.select('B8').divide(10000),
      'RED' : image.select('B4').divide(10000),
      'BLUE': image.select('B2').divide(10000)}).rename('EVI');
      return image.addBands(evi);
}

//Mask out VIs image values not in the range (0,1)
function VIsMask(image) {

  return image.updateMask(image.gte(0).and(image.lt(1)));
}

// Need to change for Landsat5 or Landsat8
var EVI_Collection = s2_QA.map(addEVI).map(VIsMask);

print(EVI_Collection)

print(for_export);

Map.addLayer(for_export.select('EVI').median().clip(AOI), {min:-1,max:1,palette:['red', 'green']}, 'EVI for_export');

var batch = require('users/fitoprincipe/geetools:batch');

// for batch export
batch.Download.ImageCollection.toDrive(for_export, 'Folder',
                  {scale: 10,
                  region: AOI.getInfo()[AOI],
                  type: 'float'});

dqf


texts: Measuring distances when crossing UTM zones?

texts:
I like to calculate distance between two UTM coordinates which are in different zones. How I can do this?
I found that we can map these coordinates in to a new projection for solving the problem. Which is the best map projection where Cartesian coordinate related calculations are applicable like UTM?
If some points me to extended UTM zones and how it helps, that will also be very helpful.


texts: XYZ Layer in OpenLayers has missing tiles

texts:
In OpenLayers I'm having a missing tile problem where it generates a missing link icon and a pink background. I found several threads including this one saying to add this code
    OpenLayers.Util.onImageLoadError = function() {this.src = '/img/blank.png';};

Sorry if this is dumb, but should my blank.png just be a 1px transparent png24? And also where should the img folder go?

Sorry, I've tried 1px, 256px, png24, png8, gif and putting the img folder and/or file in every conceivable place and I'm having no luck. Maybe you'll see something in these screenshots.

I've also tried an image with color to see if I can notice any actual change, so that's why the blank.png file looks red in the Finder thumbnail.


texts: OGR2OGR unable to initialize connection to the server for [ConnectionString]

texts:
Im trying to import a shape file to SQL Server 2017 using ogr2ogr
ogr2ogr -overwrite -f "MSSQLSpatial" "MSSQL:server=.\[SQL instance];driver=SQL Server Native Client
RDA 11.0;database=test;" "C:\Program Files\QGIS 3.18\bin\shapefile.shp" -a_srs "ESPG:32616"

But im getting the following error:
ERROR 1: Unable to initialize connection to the server for MSSQL:server=.\[SQL instance];driver=SQL
Server Native Client RDA 11.0;database=test;,

Try specifying the driver in the connection string from the list of available drivers:
SQL Server
SQL Server Native Client 11.0
SQL Server Native Client RDA 11.0
ODBC Driver 13 for SQL Server
ODBC Driver 17 for SQL Server

ERROR 1: MSSQL Spatial driver doesnt currently support database creation.
Please create database with the Microsoft SQL Server Client Tools.
ERROR 1: MSSQLSpatial driver failed to create MSSQL:server=.\[SQL instance];driver=SQL Server Native Client RDA 11.0;database=test;

NOTE: The database is already created, GDAL version: 3.1.4, shapefile spatial reference system: EPSG:32616


texts: Merging different Landsat scenes in Grass GIS?

texts:
I am new to Grass GIS.

I have 17 scenes from Landsat8 covering an area. I have created color infrared images of bands 5,4,3 and now I want to merge all scenes into a mosaic, for further analysis. I can't seem to manage to do this, neither with r.patch or with i.imagine.mosaic.

Does anyone know another way for merging?

I tried also in QGIS, I exported from Grass as tiffs, loaded them into QGIS and used the merge command, but the images become grayscale and with some white corners overlaying.

I am attaching a screenshot of my result in QGIS, at the suggestion of @Joseph.


texts: Duplicate values in combobox in QT Designer form

texts:
I have a custom form created with QT Designer on which I have several comboboxes among other stuff. All the fields have the correct names corresponding to the database.

But something weird and frustrating happens when I'm trying to use the form in QGIS. All the comboboxes except one, when I click on them, they display duplicate values of the original values I have designed them with. E.g if the combobox originally has the values 1,2 and 3, when I click on it, it displays 1,2,3,1,2,3 or 1,2,3,2,3 or 1,2,3,3,1 or some other combination of the existing values.

If I choose from the first tuple of the values (the original ones), nothing is stored in the database. If I choose a value from the duplicates it is correctly stored. The one combobox that displays only the correct values (no duplicates) doesn't store anything in the database. What's up with that?


texts: How to uncheck a layer in LayerSwitcher on map load?

texts: I have four layers in the LayerSwitcher. When the map is loaded all the layers are checked by default. I do not want to show all Layers as checked.How do I uncheck a layer when the map is loaded?


texts:
When creating a feature, I want that, if I change the value of field A that field B is also updated with the same value, like an onChange event.
For Example: Imagine you have a basket of apples and field A contains the max amount of apples that fit in the basket. Field B is the current amount of apples in the basket. So when I create the feature "basket" I want that the basket is always full. This means, that the max amount is equal to the current amount. So if I enter the max amount it should also set the same number to the current amount.
In the end, I just want to save time to tip the number twice.
Any idea to achieve this?


texts:
Since I have no way, other than using BBOX strategy with my vector layer + wfs protocol
Is it required to use BBOX strategy with WFS protocol in OpenLayers 2?

I'm wondering if there is a way to force BBOX to load all layer features, not just display data that is within the map's extent.


texts:
I'm activating the built in ArcGIS measure tool using the following code:

ControlsMapMeasureTool measureTool = new ControlsMapMeasureToolClass();
measureTool.OnCreate(_axMapControl.Object);

// set current tool to be measure tool
_axMapControl.CurrentTool = measureControl as ITool;

// deactivate measure tool
_axMapControl.CurrentTool = null;

Whenever I deactivate the tool the window opened by the tool is not closed.  How can I cause it to close?


texts: Converting LatLng coordinates to point shapefile

texts:
I have more than 50,000 LatLng coordinates defined under SWREFF99 1800 projection. I want to know, is it possible to convert these coordinates into a point shapefile ?
if yes, then how can I do this ?


texts: Trouble registering datasets with ArcGIS Server 10.1 before publishing an MXD

texts:
I'm having trouble registering a dataset while attempting to publish an MXD to ArcGIS Server in 10.1. ArcMap and ArcGIS Server are 10.1, and the MXD was originally created in a previous version of the software (9.3.1 or 10.0, I'm not sure).

The MXD includes a joined layer. The feature data comes SQL Server 2005 over a direct connection. Many of the other layers in the MXD also come from SQL Server 2005 over direct connect so I don't think this is part of the problem. The joined data comes from a separate SQL Server schema in the same instance.

When analysing the MXD before publish I am told a few things about this joined feature class. 1) it doesn't have a spatial index and 2) - more importantly - its data source is not registered with the server so its data will be copied to the map server machine. This is obviously not what I want.

The feature class's direct connection is successfully registered, and other non-joined feature classes that use the same connection do not complain about dataset registration so I know that's not the problem. The issue is coming from the separate dataset that is used for the join. Its connection parameters (SQL Server instance, db name, user and pass) have all been registered with the server but it still complains. If I right-click the issue in the analyse results and say 'register with server' (or whatever it says) a dialog shows what looks like the exactly same connection details I already registered (though ArcMap won't show me the full text of the connection string). When I click OK it tells me that there is a parameter missing from the connection properties. This error is repeated twice.

I previously saw a similar error message when trying to first register direct connections. It contained the same 'missing property' message but then also said the version of the database server was not compatible. I first interpreted this to mean that ArcGIS Server 10.1 needed SQL Server 2008, but I later discovered this actually meant that the SQL Server native client had not been properly installed on the ArcGIS Server machine. After a manual install I could successfully register direct connections.

Is there some issue with registering joined datasets like this with ArcGIS Server?
Am I just using incompatible software versions in my stack here? Please don't just guess at the answer to this second question, I'd like to have some real and genuinely useful information before changing things.
Is there some issue with publishing MXDs created in previous versions of ArcMap?

Any help very much appreciated.


texts: Where to find a list of geographic 'objects'?

texts: Is there an extensive list of geographic "objects", "entities" or "features" that can be discovered in the geographic environment, relatively static and therefore can be geo-located and represented on a map, for example building, restaurant and also park and lake would be part of this hypothetical list, as opposed to, for example apple, cat, and fun.


texts:
Is there a way in Shapely to calculate intersections of geometries with some tolerance?

Imagine we have lines that are the same, but their coordinates are written with different precision.

LINESTRING(83 55, 84 56)
LINESTRING(83 55, 84.499999998 56.4999999998)

Or imagine a line had decimal fractions, and was rounded in another program, and returned back. I saw that rounding caused unexpected geometry behaviour (like lines that should intersect, do not because of 10th digit precision), so am trying to work around this.

The problem I'm solving is to find if a route in OSRM has passed through a certain segment in OSM highways, that are imported in Python separately from OSRM.


texts: What are the main steps when building a balloon callout using ArcObjects in arcengine 10 in VS2010 with C#

texts:
I'm looking for a good BalloonCallout example / tutorial / sample.
edit:
Using C# with ArcEngine 10
edit #2:

What are the main steps when building a balloon callout using ArcObjects in arcengine 10?

Here's what I've tried so far:

This happens in a dynamic data layer

DynamicGlyphFactory dynamicGlyphFactory = (DynamicGlyphFactory)dynamicDisplay;
IDynamicSymbolProperties2 dynamicSymbolProps = dynamicDisplay as IDynamicSymbolProperties2;

IFillSymbol fillSymbol = new SimpleFillSymbolClass();
TextSymbolClass textSymbol = new TextSymbolClass();
IBalloonCallout balloonCallout = new BalloonCalloutClass();

simpleFillSymbol.Color = rgbColor;

dynamicSymbolProps.set_DynamicGlyph(esriDynamicSymbolType.esriDSymbolFill, fillSymbol as IDynamicGlyph);
dynamicSymbolProps.SetColor(esriDynamicSymbolType.esriDSymbolFill, 20f, 20f, 20f, 255);

IDynamicGlyph fillGlyph = dynamicGlyphFactory.CreateDynamicGlyph((ISymbol)fillSymbol);

balloonCallout.AnchorPoint = anchorPoint;
balloonCallout.Style = esriBalloonCalloutStyle.esriBCSOval;
balloonCallout.Symbol = fillSymbol;
balloonCallout.LeaderTolerance = 10;

textSymbol.Background = (ITextBackground)balloonCallout;
dynamicSymbolProps.TextBoxUseDynamicFillSymbol = true;
IDynamicGlyph glyph = dynamicGlyphFactory.CreateDynamicGlyph((ISymbol)textSymbol);
dynamicSymbolProps.set_DynamicGlyph(esriDynamicSymbolType.esriDSymbolText, glyph);

iDynamicDisplay.DrawText(point, someString);


texts: ST-Links convert Feature Class to DB Table error

texts:
I'm using ArcMap with the ST-Links SpatialKit Conversion tool to convert a Feature Class to a DB Table. But I'm getting an error: ERROR: 42601: type modifier is not allowed for type "geometry".

I can't seem to find anything on that error.

Feature Class (shapefile)

Geometry Type : Polygon
Projected Coordinate System : WGS_1984_UTM_Zone_21S
Projection : Traverse_Mercator
Geographic Coordinate System : GCS_WGS_1984

DB

PostgreSQL 8.3
Postgis 1.5

Any ideas? I have other tables in the DB that were generated with Quantum GIS. I also can't open the other spacial tables on ArcMap correctly. In Quantum GIS I have no problems. But I need to work with ArcGIS.

Thanks,


texts: Seeking Indian city points and village boundaries?

texts:
I am looking for Indian city/village points and village/Taluk boundaries for India.

A perfect example what I am looking for is here for the State of Haryan, produce by
Haryana Space Applications Centre.

Till now I am little lost, can some body point me in right direction and share some sources from where I can get this data for commercial purpose.

Just for update for all who are following

I have found here that Survey of India has "Village Bdy.Database (#42)". This is exactly what i am looking but till today no response from them.

Any body have some idea how to purchase this layer from them.


texts: zooming map in QGIS map canvas without changing the scale

texts:
I use QGIS to create maps mainly for printing in scale 1:100.000 to 1:25.000. It would be great, if the map in the canvas could be zoomed without changing the scale of the canvas.
The only way to get "into" the map is to change the scale. But I want to enlarge my view like with a magnifying glass. (I know, I could use magnifying glass from operating system, but I'm looking for a QGIS-Canvas solution)

Is there already a way? Is it worth to place a feature-request? Is anyone else interested in that?


texts: Not able to extract values from raster TIFF in EPSG:26912, but can in EPSG:4326?

texts:
I have an elevation file (.tif) in EPSG:4326 and point data for wildfire occurrence in the same projection. When I use this projection, I can use the following code to extract the elevation value, for the corresponding point and add that into the dataframe:
wfdf=gpd.read_file('wildfire_all_1980_2016_gdb/wf_all_1980_2016_clip.shp')
elev=rasterio.open('elevation/elevhillshade.tif')

fig, ax=plt.subplots(figsize=(12,12))
wfdf.plot(ax=ax, color='orangered', alpha=0.5)
show(elev, ax=ax)

coords = [(x,y) for x, y in zip(wfdf.DLONGITUDE, wfdf.DLATITUDE)]
wfdf['Elevation'] = [x[0] for x in elev.sample(coords)]
wfdf

However, I need this in EPSG:26912 to be more accurately represented. I have tried to reproject both layers into EPSG:26912, but then the elevation values just result in my default NoData value -9999 for every point. So I then created a whole new raster for elevation, already in the correct projection, using mosaic tool in ArcGIS. I checked in the software, and there were no NoData values, the raster is all good. However, same thing happens - I read it in, I reproject the wfdf fire occurrence point data to EPSG:26912 as so:
wfdf=wfdf.to_crs(epsg=26912)

I run the same code, and again, all values for Elevation are -9999. It only seems to read data values in the EPSG:4326 and I cannot find a solution online.
If I run the Extract Values to Points in ArcGIS for the correct coordinate system, it works also, so it's something happening in Python.


texts: Accessing one feature at time to do netwoprk analysis from ArcPy?

texts:
I am a bit stuck with my code right now.
I have a two feature classes, one for schools, one for residential buildings. So far I have built a model with the modelbuilder, whichs works perfectly as well as the exported python code. The aim is to create service areas around the school and then calculate the distance from the building, which are located in the service area to the school.
My problem now is, that I want to do this one by one. At the moment it calculates all the service areas at once.
I know I can access single features with the arcpy.da.SearchCursor but then I am stuck.

So far I got this code for my problem:

searchField = "FID"
with arcpy.da.SearchCursor(schoolPoints,searchField) as cursor:
for row in cursor:

    singleSchoolCounter += 1
    where = '"FID" = ' +"%s" %singleSchoolCounter
    #arcpy.AddMessage("Shapefile %s has %s rows " %(schoolPoints, singleSchoolCounter))
    for schools in range (0, singleSchoolCounter):
        ServiceAreaSchool = "ServiceArea" +"%s" %singleSchoolCounter
        arcpy.Select_analysis(schoolPoints,oneSchool,where)
        arcpy.AddMessage("Start calculating Service Area for School #%s" %(singleSchoolCounter))

        #Calculate the service Areas for the schools
        #create a new Service Area Layer
        arcpy.na.MakeServiceAreaLayer(networkDataset, "ServiceArea"+"%s" %singleSchoolCounter, "Meters", "TRAVEL_TO", "800, 1000", "SIMPLE_POLYS", "NO_MERGE", "RINGS", "NO_LINES", "OVERLAP", "NO_SPLIT", "", "", "ALLOW_UTURNS")
        #Add the schools
        arcpy.na.AddLocations(ServiceAreaSchool, "Facilities", oneSchool, "", "")
        #solve the Service area layer
        arcpy.na.Solve(ServiceAreaSchool)
        serviceAreaCounter += 1
        arcpy.SaveToLayerFile_management("serviceArea" +"%s" %singleSchoolCounter,outAreaFile,"RELATIVE")
        #Save solved Area to Disk
        arcpy.AddMessage("%s successfully built" %(ServiceAreaSchool))

It works until solving. There I get this error:

  ExecuteError: ERROR 030024: Solve returned a failure. Invalid
  argument. Failed to execute (Solve).

Maybe someone can help?


texts: I am trying to generate tiles on the fly in windows using mod_tile. I have followed the tutorial http://switch2osm.org/serving-tiles/manually-building-a-tile-server-14-04/ on ubuntu and successfully got it running but I want to do the same in windows. Could anyone help me by sketching up a tutorial for the above.


texts:
i have to make a software which is track a vehicle on Google map using GPS. I have almost done this project but i am stuck at a point, the client require a system that track the bus in a route. I am done route but how can i track vehicle on the route or out of route, i look few example like post GIS or spatial Index on MYSQL but i don't know where to start

In MYSQL i am create a linestring my database structure is

Field, Type, Null, Key, Default, Extra
id,int(11),NO,PRI,NULL,auto_increment
route_name,varchar(45),YES,,NULL,
from_id,varchar(45),YES,,NULL,
to_id,varchar(45),YES,,NULL,
line_string,geometry,NO,MUL,NULL,

where a line_string column is geometric column where i am store my all route coordinates now i have a lat long which is locate on this line or not have to calculate this


texts:
How can I change the field names (not the alias) in ArcGIS 10.2?
I already found this page on ArcGIS support: "Renaming shapefile fields (Defense Mapping)" but there is is not "Defense Mapping" on Customized Mode Categories In my ArcCatalog!
Without access to Defense Mapping, how can I change the field names?


texts: Create Point - measure distance to point layer- accounting for terrain

texts:
How can I - select a point / location- they pass us a location, that will be a point from one of a series of preexisting layers (which may include lat long or mailing address, but not generally. Generally it will be a preexisting point from a point layer.

We need to measure a straight line from that - selected point - to a bunch of points in 1 other layer.

That's it. - looking for an output of a table of distances and point information.

I have 10.1 advanced (with spatial analyst etc. and enterprise server.


texts:
I am currently working with HEC-GeoRAS

I am trying to import data from a shapefile into the feature class River. I would like to insert a type of data DOUBLE in the attribute ARCLENGTH of the table but I do not have the option. ARCLENGTH only allows float data.

How can I solve this problem?


texts: Is there a maximum directory or char length for file paths in ArcMap 10.0? I spent hours trying to trouble-shoot an untyped crash in both ArcMap and Catalog when trying to reclassify and subsequently delete a large raster. Moving the files to a shorter directory seems to have fixed the issue, but now I am curious if there is a certain filepath length you cannot exceed.


texts:
I created grids (A1, A2....,B1,B2,...etc.) in print layout (main properties-grids). I also use atlas. I display many points and I need to calculate location of every point within a page of atlas and a grid. It is simple to calculate location of a point within a page atlas, because atlas is driven by a grid layer. I don´t know how to simply calculate location of a point within a grid.
On a picture below black line is atlas grid frame (overlapping), red line are grids created at print layout
How can be it done?


texts:
Is it possible to change the default CRS by code?

I want to change between "Prompt for CRS" and "Use a default CRS" option by code in QGIS 2.18.


texts: working with geojson and vincent on python?

texts:
I want to import a geojson file into python so I can map it with a visualization package vincent and merge with other data in a pandas data frame.

To be specific, the said geojson file is: https://gist.githubusercontent.com/justgrimes/a63e8ec28868d8ce337e/raw/13cd7d8867e579d8b4f6b44f002a99940fba0504/DCNeighborhoodBoundariesWaPo.geojson. It's a map of DC with neighborhoods, put together by Justin Grimes.

Right now, I'm just trying to visualize this map on notebook. Here's my code:

import vincent
map=r'https://gist.githubusercontent.com/justgrimes/a63e8ec28868d8ce337e/raw/13cd7d8867e579d8b4f6b44f002a99940fba0504/DCNeighborhoodBoundariesWaPo.geojson'

geo_data = [{'name': 'countries',
             'url': map,
             'feature': "features"}]
vis = vincent.Map(geo_data=geo_data, scale=5000)
vis

but I have nothing in the output. What am I doing wrong here?

I don't yet know much about GIS and Python.


texts: QGS Project file with multi-geometry shapefiles conversion to a single KML file with styles

texts: How can we convert an entire QGS Project file with multi-geometry shapefiles into a single KML file with styles from QGIS?


texts: ArcGIS10: Publishing GPtools etc. in AddIns

texts:
I'm creating an ArcGIS 10 add-in in .NET C# which contains some tool commands and an extension class.  This extension is being launched on application start-up.

If I create any custom GP Tools, custom symbology renderers, custom toolboxes etc. in code do these still need to be registered during an install or can I launch them (as in have a custom symbology class in a layer properties dialog) dynamically from the extension class during run time?

Assume that the user of the tools does not have admin rights to their machine.

Hopefully I've made myself clear if not just ask!

Thanks
Rob


texts:
I am trying to change a title based on a feature class named Subject. They are parcel addresses, and the StrUnit field is the unit number. The issue is, not all of the parcels are apartments, so when it displays the title, it will show the address then the word "none". I tried to create a cursor to search the row where the unit number would be, if it's null it would delete the row, if not, keep it. Problem is, I'm not getting any error messages and it's still displaying None. Any ideas?

import arcpy
relateFC = r"Subject"
relateFieldsList = ["StrNum", "Street", "StrSuf", "StrUnit"]

with arcpy.da.UpdateCursor(relateFC, relateFieldsList) as cursor:
    for row in cursor:
        if row[3] == " " :
            cursor.deleteRow()

with arcpy.da.SearchCursor(relateFC, relateFieldsList) as cursor:
    for row in cursor:
        print("COMPARABLE SALES FOR ID # " + \
                        "Subject" + '\n' + '{} {} {} {}'.format(row[0], row[1], row[2], row[3]))


texts:
I want to get SQL command from ogr2ogr. My og2ogr command is:

ogr2ogr -append -f "PostgreSQL" PG:"host=localhost user=sa dbname=name password=secpsw" commm.dbf -nln "commision"

and i want it as like this:

MyOgr2ogr command > output.sql

ogr2ogr result is SQL commnand, which insert Dbf file to PostGIS. But I get not this SQL command. SQL command is hidden and I wanto to show it

For example shp2pgsql command produce SQL command and i can use shp2psql > insert.sql


texts: Remove a image from a ImageCollection if a band doesn't exist

texts:
I'm using the gee python api. The Hansen & Al. glad alert dataset is currently broken and the 20 last images are missing the 'conf19' band. In my process I need to use it anyway as a mosaic so I need to create a new Image Collection that contains all the image that have the 'conf19' band.
all_alerts  = ee.ImageCollection('projects/glad/alert/UpdResult')

#need to filter the images that miss the 'conf19' band here

var alerts_2019 = all_alerts.select('conf19').mosaic().clip(country);

How can I drop all the images that doesn't contains a specific band ?


texts: Is it possible to add vector tile layer published by GeoServer layer using ArcGIS JS API？

texts: GeoServer can publish vector tile layer in format of “application/x-protobuf;type=mapbox-vector”, and I believe the ArcGIS JavaScript API supports this format. So is it possible to add vector tile layer published by GeoServer layer using ArcGIS JSAPI?Or can the ArcGIS JS API only add vector tile layer published by ArcGIS Server?


texts:
I've imported my data to database by osm2pgrouting. This tool created network topology which contains vertices and edges between them. Problem is when I need get route to specific point which is not corresponding with any vertices of my database. For example I added screen of my route where you can see that my route is longer as should be. Red dots are vertices generated by osm2pgrouting, and blue point is my required final point by coordinates:

I select route with dijsktra algorithm:

SELECT * FROM pgr_dijkstra('
    SELECT gid AS id,
         source,
         target,
         cost_s AS cost,
         reverse_cost_s AS reverse_cost
        FROM ways',
    11545,
    19355,
    directed := true)e
    left join ways w on w.gid = e.edge
    order by e.seq

What is the best and simpliest way to cut it in same level as my required point. Geometry of my route is linestring which is composite of many points.

Is it possible to select some of the points which are nearest to my required point and cut line there?


texts: Calling gdalogr:gridinvdist in python console in QGIS

texts:
I'm trying to apply grid IDW interpolation to a point layer using python in QGIS, but I keep running into errors about one or more parameter values being incorrect... Should be a straightforward process but... What could be wrong?

import processing
processing.runalg('gdalogr:gridinvdist',
                        {"Input layer":"C:/Users/user/Desktop/test/shp/test.shp",
                        "Z field":"z",
                        "Power":2,
                        "smothing":0,
                        "Radius 1":4,
                        "Radius 2":5,
                        "Max points":6,
                        "Min points":1,
                        "angle": 0,
                        "nodata": 0,
                        "rtype": 5,
                        "Output file": "C:/Users/user/Desktop/test/tif/test2.tif"})

>>> execfile(u'c:/users/user/appdata/local/temp/tmpskr43u.py'.encode('mbcs'))
Error: Wrong parameter value 6 for parameter Max points.

I've tried playing with the values but no luck, it would return same error with the new number I've input. I've also tried commenting out one or more parameters but still won't do the trick and would create other parameter errors.


texts: How to create a point table from a remote csv file?

texts:
The current approach that I am using to load a csv file into a PostGIS enabled database on a server is the following:

-- create table
CREATE TABLE points(
 id int,
 lat double precision,
 lng double precision
);

-- copy data from csv with headers
COPY points (id, lat, lng)
FROM '/home/mydata.csv'
WITH
    DELIMITER AS ','
    CSV HEADER ;

-- Add point geometry column to table
SELECT AddGeometryColumn('public', 'points', 'geom', 900913, 'POINT', 2);

-- Populate column with point geometries
update points set geom = ST_SetSRID(ST_MakePoint(lng, lat), 900913);

-- Create a spatial index on points
CREATE INDEX idx_points ON points USING GIST (geom);

Although this creates a table with point geometries, there appear to be some restrictions when using the COPY command as you need to be a superuser. Also, I can't figure out how to upload a csv file from my laptop to the server using this approach (I need to transfer the file using ftp and then run the script).

Is there another strategy that I should be considering for generating a table with point geometries from a csv? I want to perform this entire process in PostGIS/SQL.


texts: Selecting the most recent entries from a table in ArcGIS

texts: I have a point shapefile that has many event occurrences stored for ever point.  I am trying to select the most recent event for every point.  Is there a way to do this through the "Select by Attributes" in ArcGIS?


texts: Create table in one schema, update geometry_columns View in Public schema using GeoAlchemy2

texts:
Using the SQLAlchemy ORM with GeoAlchemy2, I can create the table in the public schema and the geometry_columns view gets updated as expected and the created table also includes a gist index on the geom column automatically:

from sqlalchemy import create_engine
from sqlalchemy.ext.declarative import declarative_base
from geoalchemy2 import Geometry
from sqlalchemy import Column, Integer

engine = create_engine('postgresql://username:pwd@localhost/my_db')

Base = declarative_base()

class Aa19(Base):
    __tablename__ = 'aa19'
    id = Column(Integer, primary_key=True)
    geom = Column(Geometry('POINT', 4326))

Aa19.__table__.create(engine)

But what I want to accomplish is to create a table in one schema, while ensuring the geometry_columns view in the public schema gets updated accordingly. I can't figure this out. I thought by including a variable that sets the search_path maybe I could do this, and then the create_engine() call used above gets replaced like this:

db_schema = 'target_schema,public'
engine = create_engine('postgresql://username:pwd@localhost/my_db',
                       echo=True,
                       connect_args = {'options':'-csearch_path={}'.format(db_schema)})

But this doesn't work and I get the error message:

ProgrammingError: (psycopg2.ProgrammingError) relation "public.aa19" does not exist

That's fairly obvious and not surprising. I've also tried to set the search path to just my target_schema, hoping that it would automatically update the geometry_columns view in the public schema. But this doesn't work out either and I get the message:

ProgrammingError: (psycopg2.ProgrammingError) type "geometry" does not exist
LINE 4:  geom geometry(POINT,4326)

My question is: How can I create a table with geom column (and gist index) in one schema, AND populate the geometry_columns metadata in the public schema?

I'm new to this technique of interacting with an RDBMS using python. I've gone through the whole ORM tutorial and most of the Core tutorial, as well as having perused the docs. I'm sure the solution is in there somewhere, but I haven't come across it.


texts:
First I have a set of points in UTM 48S / WGS84 (EPSG 32748) as test.csv:

X,Y,Z
9231824.04,787648.04,711.1442

Second I defined a test.vrt for ogr2ogr to read the data as follows

<OGRVRTDataSource>
    <OGRVRTLayer name="test">
        <SrcDataSource>test.csv</SrcDataSource>
        <GeometryType>wkbPoint25D</GeometryType>
        <LayerSRS>EPSG:32748+4979</LayerSRS>
        <GeometryField encoding="PointFromColumns" x="X" y="Y" z="Z"/>
    </OGRVRTLayer>
</OGRVRTDataSource>

Thirdly I then convert to UTM 48S using EGM2008 datum as follows:

ogr2ogr -f CSV output.csv test.vrt -t_srs EPSG:32748+3855

The resulting file output.csv equals though the source:

X,Y,Z
9231824.04,787648.04,711.1442

I would expect for the Z-value to be difference. My OGR version is as follow:

ogr2ogr --version
GDAL 2.2.3, released 2017/11/20

What am I doing wrong?


texts: Selecting multiple attributes, greater and less than

texts:
I am trying to select by attributes from a polygon shapefile. I need to select areas with a temperature greater than 10 oC and lower than 15oC using gridcode but I can't seem to work out the correct formula.

Can anyone help?

The process so far has been Geoprocessing, selection by attributes,

"GRIDCODE"=  >10, <15

to which an error window comes up stating invalid SQL statement.

I am using ArcMap 10.3


texts: How to remove identical lines from layer in QGIS?

texts:
I converted polygon layer to line layer using "Polygons to Lines" tool.
Then I used "Split" function of the Networks plugin.
What I need to do next is to remove identical lines that were shared boundary of the touching polygons.
I've tried using MMQGIS plugin and v.clean (break,snap,rmline,rmdupl) but nothing works. You can see more details in attached image.

Edit: you can download shapefile here: https://drive.google.com/file/d/0B25RJIaeA41SVmNOUndWWGNlWU0/view?usp=sharing


texts: I have image .tiff (format geotiff). How I can convert the image to 3d model (terrain)? Blender or world machine can do it?


texts: Create Polygon with Length, Width, and a center coordinate

texts:
Using Shapely or like library, How can I create a rectangular polygon with known Length, Width (in meters)(Drawing over EPSG: 3857), and center coordinates (EPSG: 4326)?  The output needs to be in GeoJSON format.

import shapely.geometry
from shapely import affinity
from shapely.geometry import shape, LineString
import geopandas as gp
import geojson
import json
from ast import literal_eval
from functools import partial
import pyproj
from shapely.ops import transform

def main():
    # create geodataframe with some variables
    gf = gp.GeoDataFrame({'lat':38.6159031816694, 'lon':-90.2106388, 'width':100, 'height':60}, index=[0])

    # set default crs for dataframe
    gf.crs = {'init' :'epsg:4326'}

    # create center as a shapely geometry point type and set geometry of dataframe to this
    gf['center'] = gf.apply(lambda x: shapely.geometry.Point(x['lon'], x['lat']), axis=1)
    gf = gf.set_geometry('center')

    # change crs of dataframe to projected crs to enable use of distance for width/height
    gf = gf.to_crs(epsg = 3857)

    # create polygon using width and height
    gf['center'] = gf['center'].buffer(1)
    gf['polygon'] = gf.apply(lambda x: shapely.affinity.scale(x['center'], x['width'], x['height']), axis=1)
    gf = gf.set_geometry('polygon')
    # gf = affinity.rotate(gf, 11.60, origin='centroid')

    # return to inital input crs
    # gf = gf.to_crs(epsg = 4326)

    geopoly = gf['polygon'].to_json()
    a = geojson.loads(geopoly)
    # print(a.bbox)
    startbbox(a.bbox)

def bbox2wkt(minx, miny, maxx, maxy):
# first validate bbox values
    assert isinstance(minx,float) or isinstance(minx,int)
    assert isinstance(miny,float) or isinstance(miny,int)
    assert isinstance(maxx,float) or isinstance(maxx,int)
    assert isinstance(maxy,float) or isinstance(maxy,int)
    assert (minx < maxx), 'failed: %s is not < %s' % (minx,maxx)
    assert (miny < maxy), 'failed: %s is not < %s' % (miny,maxy)
    return "[%(minx)s, %(miny)s], [%(minx)s, %(maxy)s], [%(maxx)s, %(maxy)s], " \
           "[%(maxx)s, %(miny)s], [%(minx)s, %(miny)s]" % locals()

def startbbox(str1):
    # bbox = '-90.21153711528392,38.61548204271258,-90.20974048471568,38.6163243181534'
    bbox = map(str, str1)
    bbox = ','.join(bbox)
    rtn = bbox2wkt(*map(float, bbox.split(',')))
    data = list(literal_eval(rtn))
    gd_info = dict(type="Polygon", coordinates=[data])
    s = json.dumps(gd_info)
    g1 = geojson.loads(s)
    g2 = shape(g1)
    gf = affinity.rotate(g2, 90, origin='centroid')
    project = partial(
        pyproj.transform,
        pyproj.Proj(init='epsg:3857'),  # source coordinate system
        pyproj.Proj(init='epsg:4326'))  # destination coordinate system
    g2 = transform(project, gf)
    wow3 = geojson.dumps(g2)
    wow4 = json.loads(wow3)
    with open('data222.json', 'w') as f:
        geojson.dump(wow4, f)

if __name__ == "__main__":
    main()

Now I just need to clean it up.


texts: Java Script API V3.6 displays base maps in a chequered pattern

texts: I am currently running ARC GIS for Server 10.1 on an intranet. I have installed the Java Script API V 3.6. I am currently trying to view a cached image service using the Java Script API. When the service loads the tiles do not load in their proper places. Does any one know a fix for this?


texts:
I have 7 bands of Landsat and I would like to combine all layers into one layer.

Can ArcMap 10 do this like Erdas software can?


texts: PostGIS function testing geometries, inside and boundary

texts:
Is there a st_????(A geometry, B geometry) in PostGIS that returns true when B is inside A and possibly shares a boundary ?
I would say no after reading the doc but something may have escaped me.


texts: Web AppBuilder's accessibility changes break TAB in my custom in-panel 2D widget

texts:
I don't know if anyone else has encountered this, but I recently discovered that ever since v2.11 of Web AppBuilder, the TAB key navigation has broken in my in-panel 2D widget. (Modal and 3D are fine.)
I discovered there's a new file called client/stemapp/jimu.js/accessibleUtils.js which attempts to keep focus within the widget. The idea seems to be that if you press TAB on the last focusable node within the widget, it will "wrap around" to the first focusable node in the widget. Similarly, pressing SHIFT-TAB on the first focusable node would move the focus to the last focusable node.
In my case, the problem seems to be that when WAB tries to discover the first and last focusable nodes, my widget does not yet contain any. This results in the widget container becoming focusable and being the first and last focusable nodes. The result is when I press TAB inside a text box that I added later, the focus shifts to the frame of the widget and no amount of TAB or SHIFT-TAB will make a difference.
Note that there was no mention of this breaking change in the release notes.
Also, it appears that widgets now need to keep track of their first/last focusable nodes in order to continue to support this behaviour. I see many widgets call initFirstFocusNode() and initLastFocusNode() in order to keep up with the dynamic adding/removing of controls.
Is it expected that custom widget developers need to do this same housekeeping? This could be very complicated in our case and I'm reluctant to add so much effort in order to keep TAB working.
Is there some kind of way to opt out of this? It seems the default behaviour should be to add no classes/event handlers rather than making the widget's dom node focusable and setting its role to document.
In the short term, I have added a workaround to eliminate this behaviour from my widget. If you are similarly affected, you could do this too.
// Your `define` statement will need to include
// "dojo/_base/html" and "dojo/keys", which I refer to
// as `html` and `keys` respectively.

startup: function() {
    this.inherited(arguments);

    // (Your widget's startup code here.)

    if (this.inPanel) {
        var _this = this;
        var tabEventListener = function(event) {
            if (event.keyCode === keys.TAB) {
                if (_this.domNode.firstNodeEvent) {
                    _this.domNode.firstNodeEvent.remove();
                    delete _this.domNode.firstNodeEvent;
                    html.removeClass(_this.domNode, "firstFocusNode");
                }
                if (_this.domNode.lastNodeEvent) {
                    _this.domNode.lastNodeEvent.remove();
                    delete _this.domNode.lastNodeEvent;
                    html.removeClass(_this.domNode, "lastFocusNode");
                }
                _this.domNode.removeAttribute("tabindex");
                _this.domNode.removeAttribute("role");
                // The listener now removes itself. Our work here is done.
                _this.domNode.removeEventListener("keydown", tabEventListener);
            }
        };
        this.domNode.addEventListener("keydown", tabEventListener);
    }
}

The above code adds a keydown event handler on the widget's domNode before WAB's event handler is added. This means mine will be executed first.
When TAB is pressed (regardless of SHIFT) we do the following:

Remove the keydown event handlers from the widget's DOM node for it being the first and last focusable node.
Remove the related css classes.
Remove the tabindex attribute from the widget's DOM node to prevent being able to tab onto it.
Remove the role attribute from the widget's DOM node. (It does have controls now, so it doesn't need that role.)
My event handler then removes itself as its work is done.

The end result is I can now TAB and SHIFT-TAB through my widget again. I don't get the wrap-around focus but then I also don't get focus stuck on a div, so I consider it a win.


texts: Is there a Python script to identify which database features are in Map/Feature Services?

texts: I would like to identify feature classes that are present in my map services.  Is there a way to automate this with python?


texts: Import point data from amazon s3 into QGIS

texts:
I'm working with spatial/temporal data, stored in s3 (AWS). The data is in the format:

ID,    Date,                 Lat,    Lon
23424  2017-02-03T23:15:15Z  23.434  -5.2343
23424  2017-02-03T23:20:26Z  23.4356 -5.2637
etc...

What I want to do is import this data into QGIS. It is very large, several gb.

My current method is to use an Athena query to download the relevant data, as a csv. I then import the csv into a PostgreSQL database on my local machine, change the lat/long columns into a point geometry column, and import into QGIS as a POSTGIS layer.

This is time consuming, and also uses up space on my local machine.

Is there a way of connecting to the data stored externally on s3 instead?

I am pretty new to python, and I haven't used the python console much.


texts: Shapefile, clip or project first when both are necessary

texts:
This is a question mostly out of curiosity, but I do it enough if there were a better choice I would stick with that.

Using ArcMap 10.6, when a shapefile needs both clipped and projected is there an order in which one workflow is faster? Or would I just need to experiment to find out?


texts: "Singleparts to multipart" vs. "Dissolve"

texts: It seems to me that the vector tools "Singleparts to multipart" and "Dissolve" do the same things. Can someone confirm this, please? Or do I overlook something?


texts: ArcCatalog crashes when publishing a composite locator service

texts:
I created two address locators in a shared file gdb which I then published as separate geocode services. After performing several test searches in ArcMap against the gdb locators and the services, everything seemed to be working correctly.

I then created a composite address locator in the file gdb, composed of the two separate geocode services and then published the composite locator on the server. When I get to the Service Editor and Analyze, the only warning I get is that the composite locator will be copied to the server, even though the gdb is accessible to the server. This is a minor issue though, as covered in the answer to this question.

The problem occurs after clicking Publish. ArcCatalog hangs for a few seconds, then crashes. I rebooted the machine and tried again, and the same thing happened. I then tried from a different admin account, and it crashed again. A post on the ArcGIS forum indicates that this might be an issue in ArcGIS itself when trying to include online services in a composite locator service.

I'd like to know if this is possibly a bug, or if there is some workaround or settings that must be changed. I left all the default settings unchanged. I'm using ArcGIS for Desktop 10.1 and ArcGIS for Server 10.1, with SP1.


texts:
Could someone shed some light into a method to link breaking outline borders without much distortion from the original like e.g. convex hull does?
Convex hull alternative ideas are of course welcome to be tried even if they have similar results.

Below two screenshot of what I try to remove and achieve a more continuous result:


texts: Can ArcGIS only merge two shapefiles each time?

texts: I need to merge multiple point shapefiles (in the same directory) into one point shapefile. After inputting these shapefiles and merging it, I found the final output only contained the data from two shapefiles. How to merge multiple shapefile in ArcGIS at once instead of merging only two at each time?


texts:
Is ServerVector broken in OpenLayers 3.5?

Just running the following gives an error:

var wfsTest = new ol.source.ServerVector({
    format: new ol.format.GeoJSON()
});

The error is: "Object doesn't support this action"

Am I doing something wrong??  It used to work fine in previous OpenLayers 3.x versions.


texts: How to use inbuilt tools in QGIS Modeler?

texts: Is it possible to incorporate inbuilt tools like, for eg: vector grid of Research tools in QGIS model builder ?


texts: Improve performance of QgsVectorLayerDirector.makeGraph() in PyQGIS 3

texts:
Porting a plugin that makes excessive use of PyQGIS networking capabilities I had to realize, that building the graph (makeGraph() method of class QgsVectorLayerDirector) in PyQGIS 3 (https://qgis.org/api/qgsvectorlayerdirector_8cpp_source.html#l00169) is about 4 times slower than in PyQGIS 2 (https://qgis.org/api/2.18/qgslinevectorlayerdirector_8cpp_source.html#l00128) (Vector layer with ~18k line features, graph with ~38k vertices and ~73k edges). I tried to find the reason for this from the source codes referred to, but failed because of my poor C++ skills.

Since the documentation seems not to be updated for the new PyQGIS 3 classes (https://docs.qgis.org/testing/en/docs/pyqgis_developer_cookbook/network_analysis.html), are there any recommendations to clean up input vector or sth. like that to make makeGraph() significantly faster? Or is this an issue tracking case?


texts: Search and zoom to street stretch: By NAME, FROM, TO

texts:
I have a streets feature class that has name, from, and to fields. Using this feature class, I want to zoom to locations such as these:

Main St, from 1st Ave to 2nd Ave (single feature)
Main St, from 2nd Ave to 5th Ave (multiple features)

Note: My street names aren't actually numerical in reality. I just used street names like 2nd Ave and 5th Ave because it's easy to imagine what comes in between.

I can do a couple of similar things using the Find tool:

Search for a street by name
Search for a specific address

However, these options aren't the same thing as searching for a street stretch, using the from and to in conjunction with the name.

How can I search & zoom to a street stretch using the name, from, and to fields?

I don't know if it's relevant or not, but there is a fourth field, sort (integer) that orders the features from west to east, south to north.


texts: Why were Series variables of ModelBuilder removed at ArcGIS 10.0?

texts:
I've read a number of posts regarding the sequential batch processing of a model in ArcGIS.  In particular this one, because I'm having the exact same issue.

Recently, I found an old 9.3 help page that shows that ArcGIS at one time supported a type of value called a series which is missing from the current 10.2 dialogue in model builder.

It appears that what I need to achieve my sequential processing of multiple input parameters for my model is what was once called a "series value" in 9.3.

I'm just wondering if anyone knows what happened to this and if there is some way to achieve the same thing in 10.2?

It's hard to believe that Esri would take functionality out of more recent versions.


texts: Setting all the values in a attribute field to the same value with PyQGIS

texts:
When I use  "Update existing fields" in "Field Calculator" in QGIS. It takes only a second to change all the values in a attribute field to the same value. However when I try to do that with PyQGIS, QGIS crashes.
I have a layer with line geometry with 82800 features inside. I want to change all feature's value inside 'visited' column to '0'. Here is how I try to do that,
layer = iface.activeLayer()

lines = [feat for feat in layer.getFeatures()]

layer.startEditing()
for line in lines:
    line["visited"] = 0
    layer.updateFeature(line)
layer.commitChanges()

Can I do that like I do that with field calculator via using PyQGIS (with expressions or something)?
QGIS version
3.16.0-Hannover


texts: Replacing values with non-English characters in attribute table field using ArcGIS field calculator and python parser?

texts:
I'm new to python and trying to write a python script using the field calculator in ArcMap for a specific field in an attribute table that will calculate and place the values in the other field. My Pre-logic Script Code:

def calcVal(inVal):
    values = {
"BABBBBBBBB":1,
"CCCCCCCCCC":2,
"DDDDDDDDDD":3,
“EEEEEEEEEEE”:4,
etc...}
    if inVal in values.keys():
        return values[inVal]
    else:
        pass

Character “A” in 1st condition ("BABBBBBBBB":1) is non english character. Field calculator does not recognize this text value. Coding is UTF8.

Is a special symbols in python with who I can substitute this non-english character, like % or *?

Or, it is possible to configure this code to recognize UTF8 coding?


texts: How to use Field Calculator to calculate area as percentage?

texts:
I am searching for python code to use in Field Calculator to calculate area as a percentage in relation to another field?

The area percentage needs be calculated based on the total area for each "Basin_ID" separately.


texts: pgr_drivingDistance How to split edge if max distance is reached halfway through the edge (PGROUTING)

texts:
I would like to generate a catchment within a network, up to a max cost, using the pgr_drivingDistance function. This returns all the nodes within this distance - but does not provide for partial edges to be split and terminated when the max cost is reached part of the way through the edge.

This is related to the answer in Return nodes along portions of edges with pgrouting but I am repeating it here, as I wonder whether the situation has progressed in the last three years. A duplicate is here PgRouting - How to clip links when reaching max costs? but with unsatisfactory answers that are fairly old.

A visualisation of a similar process through interpolation (post-hoc)is here: https://anitagraser.com/2017/09/11/drive-time-isochrones-from-a-single-shapefile-using-qgis-postgis-and-pgrouting/ - yet, it does not retrieve the part-way points.


texts: OpenLayers forEachLayerAtPixel and MapServer .MAP file COLOR

texts:
I'm trying to use map.forEachLayerAtPixel in OpenLayers to see what layers are clicked, but I have a layer in the background that is a grid served by Mapserver with a grid.map file that is just a bunch of outlines.

The grid.map file has an OUTLINECOLOR but no COLOR to fill it so it's transparent inside the polygon.

However, this code does not show me as having clicked inside the grid.map polygon when I do:

map.forEachLayerAtPixel(evt.pixel, function(layer) {
  var title = layer.get('title');
  console.log(title);
}

If I put a COLOR in the grid.map file it works great but I'd like to keep it transparent.

EDIT: I've tried 2 different STYLE blocks as noted in MapServer - Show Polygon Layer with non-transparent border and transparent background and still have the same issue where it works if I set OPACITY to 1 but not as soon as I set it to 0.

What am I missing to keep a polygon transparent but still be recognized by map.forEachLayerAtPixel?


texts: Is there a way to show roads connected but keep them on separate layers (with separate markings, line styles etc.)?

texts: In QGIS 2.0.1 Dufour, is there a way to connect different layer roads? Say i have 2 kinds of roads on 2 layers. One is main roads, the other is smaller ones. At the moment they go over each other depending on which of the layers is over another. Is there a way to somehow connect them but keep them on separate layers (with separate markings, line styles etc.)?


texts:
According to the page on Google Code (https://code.google.com/p/libkml/downloads/list), the last update was Feb 2010.

Is this library still being updated?


texts: SRID in parametric SQL View in Geoserver

texts:
I want to create parametric SQL View, which provide feature collection of my route, which bbox is in depending on entered parameters. I add some SQL Script which is returned some columns, and one of them is geometry column of type LineString with SRID 4326.

When I tried refresh atributes and check Guess geometry type and srid geometry was not set up as you can see on screen:

Type is still Geometry not Linestring, and SRID is -1. So my questions are:

Should I set this manualy to LineString? And if yes, should I set SRID of my column in database (4326) or that SRID what I need (900913).
Geoserver layer display this area which is covered in bbox, so if I understad, I have to set bbox of whole area, which is my route could contains? So if my data are whole planet so I need set bbox of whole planet or what ?


texts: Click on a point and perform a geoprocessing tool?

texts:
In ArcMap, I have created a Python Add-In tool that clicks on a point anywhere on the map, then uses the coordinates as an input to a geoprocessing tool that runs a select by location and buffer of various other parameters. It works very well and also will add and remove data from the map.

I am looking to recreate this tool in ArcGIS pro. Is there a way that I can do something like this without having to make some sort of AddIn? I have heard FeatureSet may do something like this. I can change up the words in the Geoprocessing tool to match the ArcGIS pro language, however, am stumped with being able to click on the map to get an input without having to do some sort of Visual Studio SDK ArcGIS pro workaround.


texts: Finding volume and lowest point in subbasins using GRASS?

texts:
I have created 156 subbasins from initial raster.
Now for any subbasins I need to have the volume of raster and the lowest point and insert in second moment in a table. I use Grass Gis software.

What is the better way for it?


texts: Properties of Google Earth Engine GeoTIFF export

texts:
I am adding a couple of properties to an export, named 'elev' and 'threshold'. My objective is to have those properties embebed in the exported GeoTIFF, however, when I open the exported GeoTIFF and try to find those properties using ArcMap or ArcCatalog I cannot find them. Anyone knows how to properly export/display those properties on exported TIFFs?

// Export the result.
  var nameexp1 = 'Mask_' + start + '_' + end;
  Export.image.toDrive({
  image: result.int().unmask(-9999).set('elev',demval,'threshold',mndwival),
  description: nameexp1,
  scale: 30, maxPixels: 228612704,
  region: ROI
});


texts: Modify dataframe extent using data driven pages and arcpy.mapping

texts:
I'm customizing an automated mapping method using arcpy.mapping and data driven pages.

Essentially I have a script that loops through each page, and modifies the location of the inset map, legend, and other map or graphic elements.

I've got everything working, except for some maps, it is necessary to change the extent of the main dataframe so that the inset map, or legend does not overlap with the region of interest.
As an example, the data frame may need to shift 500 metres North and East, so that there is no overlap. This is what I can't get to work.

I'm using a for loop, with conditional statements to make the necessary adjustments.  An example condition is:

if (pgIndex == 4):

    # Set position of the Legend
    elLeg = arcpy.mapping.ListLayoutElements(tempMap,"LEGEND_ELEMENT", "")[0]
    elLeg.elementPositionX = 0.0186
    elLeg.elementPositionY = 0.4124
    # Move scalebar
    scaleBar = arcpy.mapping.ListLayoutElements(tempMap, "MAPSURROUND_ELEMENT", "Alternating Scale Bar")[0]
    scaleBar.elementPositionX = 0.1
    scaleBar.elementPositionY = 0.1
    #Move North Arrow
    NorthArrow = arcpy.mapping.ListLayoutElements(tempMap, "MAPSURROUND_ELEMENT", "North Arrow")[0]
    NorthArrow.elementPositionX = 0.75
    NorthArrow.elementPositionY = 0.1
    #Move Legend NeatLine
    Neatline = arcpy.mapping.ListLayoutElements(tempMap, "GRAPHIC_ELEMENT", "")[0]
    Neatline.elementPositionX = 0
    Neatline.elementPositionY = 0
    #Move inset map position
    inset = arcpy.mapping.ListDataFrames(tempMap, "LocatorDF")[0]
    inset.ElementPositionX = 0
    inset.ElementPositionY = 3.25
    #Set extent for main data frame (THIS IS WHAT IS NOT WORKING)
    mainMap = arcpy.mapping.ListDataFrames(tempMap, "MainDF")[0]
    extent_4 = arcpy.Extent(423017.450685167,4991568.14747271,438257.481165228,4976328.11699265)
    mainMap.extent = extent_4


texts: Creating ArcMap point feature with ascending number symbol within it?

texts:
I made a map of a planned route, at each of the routes stops (total of 15) I would like to have a point feature with a number in the middle of the point (think pool ball with the number on the ball) to represent that particular route stop.

Does anyone know the most straight forward/simple way I can go about making something like this?


texts: I am stuck in a seemingly simple but strange thing. I am finding sinks in an area by LiDAR data. after finding the sinks, I want to determine their elongation and for this purpose I use minimum bounding rectangle tool. I have the area and the length of the rectangles, and I divide the area by length to obtain width. after having the width, I divide the length by width. But the problem is that I get a very high number for ELG and it does not look like it. I will attach the pictures .sink


texts: Is there a "Identity (Analysis)" tool for QGIS?

texts:
In ArcGIS, there is a tool called "Identity." Say you have a bunch of points. Some points fall in polygons of interest, some don't. Each polygon has a unique ID. This tool will add a field to the points and assign the unique ID of the polygon that contains it. Is there a tool that does this in QGIS?

This question does not ask the same question I am.


texts: Howto setup and use a remote GIS workstation (QGIS,GRASS,PostGIS)?

texts:
I'm using QGIS (and sometimes GRASS and PostGIS) on my 2009 MacBookPro for my thesis project.
I was wondering (during an endless batch generalizing in QGIS) if I can set up a server so I can work from the library/university using the power of a faster home desktop workstation.
So far I understood:

I can access a remote PostGIS server via QGIS
I can run GRASS via SSH, even with GUI
GDAL has PostGIS driver, but it's only for querying and inserting data, the rest of the work is done by the machine that issues the command (or am I wrong?)
Basic editing of a PostGIS can be done via QGIS, and I can issue SQL commands to make some more advanced editing

What I would like to figure out:

Can I use QGIS plugins remotely? is VNC the only way to accomplish that?
is there some documentation for PostGIS+QGIS other than creating and visualizing maps (I found only these types of tutorial)? which tools/plugins can I use, and which of them can be run remotely?

UPDATE:
you all confirmed my thoughts about QGIS, if I want to use the processing power of the server I have to control it via VNC/RDP/NX.
as for my second question, let me rephrase it:
If I use postGIS+QGIS and I want all the processing be done on the server, do I have to (learn and) use only SQL queries with postGIS functions?


texts: How to Use an ArcGIS Tile Package with Desktop

texts:
I did an ArcGIS Online training course a while back and learned a little bit more on Tile Packages. I followed the instructions, how-to-create-a-tile-package
and created a tile package for some aerial photography. The process completed and I now have a nice little 68GB tile package.

I want to use this package locally in ArcGIS Desktop, but whenever I add AirPhoto.tpk to my ArcMap it will freeze and stay like that. I have a pretty powerful PC and the resources aren't maxing out when I try to add it, it just seems to freeze the ArcGIS desktop. I've tried waiting about 12 hours for the process to complete, but it still stays frozen.

I've been able share a different package and store it on ArcGIS online, then bring it back to desktop. That option doesn't really work for the aerial photos as I would max out my credits storing it online. I've got ArcGIS for Server here and think I will try adding it there and sharing it.

What am I doing wrong here? Do I need to unpack it or something and if so what is the point of a tile package? I am currently running 10.3 Advanced, does an upgrade to 10.4 solve this?


texts: What is the latest update to the positions of the tectonic plates? Does anyone have a most up-to-date shapefile or vector of the tectonic plates?


texts: Weird result in pgRouting with Dijkstra's algorithm

texts:
I'm using Dijkstra's algorithm to perform pgRouting with MapServer. I've stored the following function in PostGIS:

CREATE OR REPLACE FUNCTION public.routing_function(
IN x1 numeric,
IN y1 numeric,
IN x2 numeric,
IN y2 numeric)
RETURNS TABLE(seq integer, path_seq integer, node bigint, edge bigint, cost double precision, agg_cost double precision, street text, geom text) AS
$BODY$
WITH dijkstra AS(
SELECT * FROM pgr_dijkstra('SELECT id, source, target, time AS cost FROM road_network',
--source
(SELECT source FROM road_neetwork ORDER BY geom <-> ST_SetSrid(ST_MakePoint($1, $2), 2100) ASC LIMIT 1),
--target
(SELECT target from road_network ORDER BY geom <-> ST_SetSrid(ST_MakePoint($3, $4), 2100) ASC LIMIT 1),
-- undirected
false)
),
the_geom AS(
SELECT dijkstra.*, b.street, ST_AsText(b.geom) FROM dijkstra INNER JOIN road_network AS b on dijkstra.edge = b.id
)
  SELECT * FROM the_geom WHERE CASE WHEN 'source' != 'target' THEN true ELSE false END;
$BODY$
LANGUAGE sql VOLATILE

It's working and I'm getting a result. But I'm wondering if the result is reliable, because instead of "deciding" to travel around the city (yellow line) which is faster, as the "cost" is time, the algorithm takes me THROUGH the city, as shown in the picture (the arrow indicates the street with the desirable route and probably the fastest). Is it normal ?

EDIT: One more example. The arrow indicates the desired route

And another one. Even more weird...!!!!


texts: Sorting by row already pre-selected polygons by location by line?

texts: I have two layers. One with small Polygons and one with a line. I want to select all Polygons which intersect with the line. Additionally the selected Polygons should be sorted by the order of the line, as a line has always a start and a end. The first polygon that intersects at the start of the line should be on top of the attribute table. The second one on second place and so on.


texts: How to import a postgis database in ubuntu (osgeo live 7)

texts: How does one import a postgis database in the ubuntu/osgeo live environment? In windows I would copy the backed up database file into the postgres bin folder and then use cmd to restore the databse. As I am new to the ubuntu/osge live environment I am not sure how to go about restoring a postgis database.


texts: QGIS - CRS setting for Slovenia national agency data

texts: This is just a humble information to my (newbee) fellow slovenian users of QGIS regarding correct setting of CRS for our country. For data provided by slovene national agency (GURS ...), use CRS marked as MGI 1901 / Slovene National Grid EPSG:3912 and data should project, export etc fine. I was quite lost some time and almost give up on QGIS as I had so many problems using wrong EPSG settings provided internally by QGIS supposed to be "Slovene national grid D40, SNG D96, EPSG:4765, EPSG:3794 etc". Hope it will work for you, too.


texts:
I'm currently using a div icon to show a small graph as a marker on a leaflet map. The graphs change dynamically based on the day of week using data from a database, which changes the height of the graphs in the div icon HTML.

As the graphs change height I need to change the value of the icon anchor so the graphs remain attached to the correct location at their bottom left hand corner, using the heights of the graph to set this. Does anyone know of the best way to access and modify the anchors of div icons/markers once they've been added to the map? (I would prefer not to have to remove the makers and re-add them if possible)


texts: Converting VMAP to shapefile using GDAL?

texts:
This is my first time trying to convert a VMAP to a shapefile and the examples online aren't as straightforward to me.

Any guidance?


texts: Orientation of segments of a polygon in ArcGIS

texts:
Using the longest side of a polygon, I would like to get the orientation of each side of the polygon with the longest side orientation as the reference.

I have attempted to get the general orientation of each segment of a polygon using:

def minBG(inFC, outFC):
     arcpy.MinimumBoundingGeometry_management(inFC, outFC, geometry_type="RECTANGLE_BY_AREA", group_option="NONE", mbg_fields_option=True)

The result gives MBG_Length and MBG_Orient in the attribute table of each polygon but unable to proceed from here.


texts: Tiles not displaying in ArcGIS Online at furthest zoom level

texts:
In an application I am "building" in ArcGIS Online (not server), our starting view is a global level. The issue is that Zoom Level #3 (example) cuts off part of our dataset in the extremes--alaska and the pacific. Zoom level #2 causes the edges of the map to turn null-colored (pic 1), even though the tiles appear to be loading in the cache (pic 2).

Grey bars on side

Tiles appear to load but are not seen

Is there a way to address this null space in ArcGIS Online?

I made a posting on geonet.esri.com. Hopefully it gets some attention.


texts: Universal Datum Transformations between Rasters in ArcPy

texts:
Problem: Translating a set of various rasters with different datum into the same datum using ArcPy

My Solution: Run rasters through a for loop, check their datum with an 'if' statement, run rasters with an unwanted datum through CreateCustomGeoTransformation(), transform these datums to desired datum.

I do not seek advice on how to do this with vector files.  These are raster files and they behave differently.

Question:
How do I make my fourth parameter, customGeoTransfm,  in the tool CreateCustomGeoTransformation() to accept any input datums and run a transformation to a single output datum, if I have to select set translations on the x,y,and z axes?

Script

GCS_Datum All rasters will be transformed into GCS_datum, string input.

geoTransfmName represents the name of the transformation(s)

inGCS represents the datum of each raster in the list in turn

outGCS represents the datum each raster will be transformed into

customGeoTransfm outlines transformation <-- this confuses me

inport arcpy

# start user inputs ----------------
GCS_Datum = "D_North_American_1983"
inFolder = "pathway to workspace"
# end user inputs ------------------

arcpy.env.workspace = inFolder
rasterList = arcpy.ListRasters("*", "ALL")

for raster in rasterList:
    geoTransfmName = "thisisatest"
    inGCS = arcpy.SpatialReference(raster)
    outGCS = arcpy.SpatialReference(GCS_Datum)
        #how to make customGeoTransfm universal?
    customGeoTransfm = "GEOGTRAN[METHOD['Geocentric_Translation'],PARAMETER['X_Axis_Translation',''],PARAMETER['Y_Axis_Translation',''],PARAMETER['Z_Axis_Translation','']]"
        #how to make customGeoTransfm universal?

    sr = arcpy.Describe(raster).spatialReference

    if sr.GCS.datumName != GCS_Datum:
    arcpy.CreateCustomGeoTransformation_management(geoTransfmName, inGCS, outGCS, customGeoTransfm)


texts:
I'm new to ArcObjects, not a GIS professional so please be patient with my lack of GIS knowledge.

I have an .Net Web application with ESRI.ARCGIS.SOAP references, looking to update them with rest.

How can I use REST in below code?

Update:
I tried to use arcgisruntime, but that was not supported in Asp.Net Web application

ESRI.ArcGIS.SOAP.GeometryServerProxy serviceURL= new ESRI.ArcGIS.SOAP.GeometryServerProxy();
serviceURL.URL="Some Service URL";

ESRI.ArcGIS.SOAP.SpatialReference inputSpatialReference = new ESRI.ArcGIS.SOAP.ProjectedCoordinateSystem();
 inputSpatialReference.WKID = 3857;
 inputSpatialReference.WKIDSpecified = true;

ESRI.ArcGIS.SOAP.SpatialReference outputSpatialReference = new ESRI.ArcGIS.SOAP.GeographicCoordinateSystem();
outputSpatialReference.WKID = 4326;
outputSpatialReference.WKIDSpecified = true;

PointN pnt = new PointN();
 pnt.X = X;
 pnt.Y = Y;

Geometry[] inputGeometry = new Geometry[] { pnt };

ESRI.ArcGIS.SOAP.GeoTransformation geoTransformer = new ESRI.ArcGIS.SOAP.GeoTransformation();
geoTransformer.WKID = 2236;

ESRI.ArcGIS.SOAP.EnvelopeN envelope = new EnvelopeN();
envelope = null;

Geometry[] outputGeometry = serviceURL.Project(inputSpatialReference, outputSpatialReference, false, geoTransformer, envelope, inputGeometry);

foreach (PointN item in outputGeometry)
 {
xyCoord= item.X.ToString() + "," + item.Y.ToString();
 }


texts: Data dependent override on colour and linewidth of line pattern fill not working after upgrade to QGIS 2.12

texts:
I have recently upgraded to QGIS 2.12 which has resulted in a problem with one of my projects, whereby the data dependent override of the line width and colour of a line pattern fill and SVG marker line pattern has ceased functioning.

I have looked in the QGIS 2.12 changelog and didn't find anything that should have affected this.

What should happen is that the colour and linewidth of the pattern fill should vary according to the topographic curvature and the average slope within a vector segment.

a previous version with working colours on my blog

I use the following expressions:

for colour: case when RAT_LngtCrv_Mean > 0 then color_rgb(min(255,255*5*RAT_LngtCrv_Mean),0,0) else color_rgb(0,0,min(255,-255*5*RAT_LngtCrv_Mean)) end

for linewidth (using map units - metres): case when RAT_Slope_deg_Mean > 2.25 then max(20*((RAT_Slope_deg_Mean - 2.25)/ 90),0.4) else 0.0001 end

for size of the SVG arrow: case when "RAT_Slope_deg_Mean" > 2.25 then max(75*sqrt((RAT_Slope_deg_Mean - 2.25)/ 90),5) else 0.0001 end

for angle (working): 90 - RAT_Aspect_deg_mean

What actually happens is this, with neither the linewidth or colour taking the data dependent values:

I have also raised this as an issue.

A fairly minimal example qgis project file reproducing this issue can be found here on my Dropbox


texts: QGIS - Can't select Inputs for Select by Location algorithm in Graphical Modeler?

texts: All the other questions I've seen so far are a step ahead of me. I've defined an Input for a vector layer, and am trying to use the Select by Location Algorithm. I typed the name of my Input, but the modeler doesn't seem interested in connecting it. Any ideas?


texts: Adding value to newly added field in Routing model?

texts:
I have a model, that I created a routing service from.  I use the following tools in the model:
Make Route Layer
Add Locations
Solve
Add Field to Analysis Layer
Select Data

I need a unique identifier each time the model is run in async mode, to be associated with that particular job (or execution).  I added a field called RouteId to the Network Analyst Layer (or resulting "Route" layer). I made this "name" a model parameter.
When I execute the service from my C# code, I can assign a value to the "name" of the field, but I can't find a way in the code to place a "value" in that new field.
Here is a snippet of the code that I am using:

//Create a new GP String with parameter name of RouteId
 GPString gpIdString = new GPString("FieldName", "RouteId");
 parameters.Add(gpIdString);
 ExecuteAsyncMethod(gp, parameters);

Once it goes into the gp.GetResultDataCompleted method, I can get at the field "RouteId" but the value is null.

How can I put a value into the RouteId field?


texts:
I need to find a free satellite WMS which is precise enough to reach a zoom level of 17 minimum. Precision is an important asset in my application.

I am working with Leaflet and I am aware of the existence of this tool. But I am actually looking for users stories about different tools to help me make a choice!

Big companies are not an option. They have changed their billing politics and screwed up my app.


texts: Create single VRT from images with the same name

texts:
I have a list of rasters some with the same date that I have to turn into VRTs with the same date too. The problem is that when I create a loop in the Gdal list it overwrites the VRTs that have the same name. While I want it to add them in the same VRT. E.g:
root = 'F:...\\MSAVI'
out_vrt = 'F:...\\out_vrt\\'

raster = [os.path.join(root, file) for root, dirs, files in os.walk(root)
               for file in files if os.path.splitext(file)[1] == '.tif' ]
raster
['.../33TVE6/S2_20190712_079_33TVE6_B.tif',
 '.../32TQQ3/S2_20190712_079_33TQQ3_B.tif',
  ...]

dates = []
for file in raster:
    dates.append(file.split('\\')[6].split('_')[1])

dates
['20190817',
 '20190901',
 ...]

for day in dates:
    for file in raster:
        data_path = file.split('\\')[6].split('_')[1]
        name = 'S2x2A_' + day + '_MSAVI.vrt'
        tile = file.split('\\')[6].split('_')[3]

        if day in data_path:

            my_vrt = gdal.BuildVRT(out_vrt+ name, file)
            my_vrt = None

The output is a VRT with only one raster:
<VRTDataset rasterXSize="4001" rasterYSize="4001">
  <SRS dataAxisToSRSAxisMapping="1,2">PROJCS["WGS 84 / UTM zone 33N",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",15],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","32633"]]</SRS>
  <GeoTransform>  4.3500000000000000e+05,  1.0000000000000000e+01,  0.0000000000000000e+00,  4.4983300000000000e+06,  0.0000000000000000e+00, -1.0000000000000000e+01</GeoTransform>
  <VRTRasterBand dataType="Int16" band="1">
    <NoDataValue>-32768</NoDataValue>
    <ColorInterp>Gray</ColorInterp>
    <ComplexSource>
      <SourceFilename relativeToVRT="0">F:...\33TVE6\S2_20190712_079_33TVE6_B_MSAVI.tif</SourceFilename>
      <SourceBand>1</SourceBand>
      <SourceProperties RasterXSize="4001" RasterYSize="4001" DataType="Int16" BlockXSize="4001" BlockYSize="1" />
      <SrcRect xOff="0" yOff="0" xSize="4001" ySize="4001" />
      <DstRect xOff="0" yOff="0" xSize="4001" ySize="4001" />
      <NODATA>-32768</NODATA>
    </ComplexSource>
  </VRTRasterBand>
</VRTDataset>

How can I add rasters with the same date to get a VRT with more rasters without overwriting the previous VRT?


texts: Accessing information about ArcSDE historical markers via arcpy or SQL?

texts:
I would like to be able to access the list of historical markers that is presented in the "Change Version" dialog in ArcMap - https://desktop.arcgis.com/en/arcmap/latest/manage-data/geodatabases/changing-versions-in-arcmap.htm

I have searched the ArcSDE system tables, but I am unable to find this information - http://desktop.arcgis.com/en/arcmap/10.3/manage-data/gdbs-in-oracle/system-tables-oracle.htm

Ideally, I would like to be able to access the name of the historical marker and other pertinent information like the creation date/time - much like the information that is stored in the GDB "VERSIONS" table for transactional versions.


texts:
I have two tables: points(gid,id,datum,x,y,geom) and hexa_grid(gid,geom).

I must create new table hexa_grid_pip (gid,2008_01,2008_02 ... 2011_12, geom)
In column 2008_01 ... 2011_12 must be number points in polygon.
Columns 2008_01 ... 2011_12 are shortcuts of months.

Is it possible to do this by using script in PostGIS, or do I have to count month by month in QGIS?


texts: Clip inside of GeoTiff using shapefile in QGIS

texts:
I'm looking to distract the red figures from the background imagery in the attached image. Normally i would use the clipper command in QGIS, but this only works if i only want to have the red figures part of the background imagery.

Can someone tell me how to get the red figures cut out of the background?

The background image is a GeoTiff and the red figures are Shapefiles


texts:
in QGIS I´d like to label geometries with it´s layer name.

Example: Name of Polygon Layer "testlayer" -> each Polygon is labeled with "testlayer"

Is there some way to access the current layer name in rules like:

layer_property(_self, 'name')


texts: Seeking global (worldwide) list of street abbreviations?

texts:
I'm looking for a list of street and road names abbreviations (e.g Street -> St, Avenue -> Ave. and so on). All I've seen are country-specific lists. Is there a global one, or at least a dataset to extract it from?

I've found this list: http://wiki.openstreetmap.org/wiki/Name_finder:Abbreviations, and it's quite good, but obviously lacks countries/languages and even items for some of sections.

the OSM "short_name" tag is of no use cause it's absent for practically every road.


texts: Is PostGIS / PostgreSQL UK Public Sector Network (PSN) compliant?

texts:
Is PostGIS / PostgreSQL UK Public Sector Network PSN compliant?

... specifically so it can be used in a UK local government council?


texts: How can I select many rows from a table?

texts:
I want to create an action that select many rows from PostGIS table to show many lines in QGIS.

How can I select many rows from table with Python to add these rows in a layer?


texts: Applying QGIS Raster Calculator upon two different Rasters?

texts:
I want to Apply some calculations over rasters ,For which using Following Code.

  import time
  from qgis.analysis import QgsRasterCalculator, QgsRasterCalculatorEntry
  rasterpath='inputfile1'
  bohLayer = QgsRasterLayer(rasterpath,"temp")
  QgsMapLayerRegistry.instance().addMapLayer(bohLayer )
  entries = []
  # Define band1
  boh1 = QgsRasterCalculatorEntry()
  boh1.ref = 'boh@1'
  boh1.raster = bohLayer
  boh1.bandNumber = 1
  entries.append( boh1 )

  rasterpath='inputfile2'
  bohLayer1 = QgsRasterLayer(rasterpath,"temp1")
  QgsMapLayerRegistry.instance().addMapLayer(bohLayer1)

  # Define band2
  boh2 = QgsRasterCalculatorEntry()
  boh2.ref = 'boh1@1'
  boh2.raster = bohLayer1
  boh2.bandNumber = 1
  entries.append( boh2 )

  # Process calculation with input extent and resolution
  calc =QgsRasterCalculator('Float(boh1@1-boh@1)/Float(boh1@1+boh@1)','outputfile','GTiff', bohLayer.extent(), bohLayer.width(), bohLayer.height(), entries )
  calc.processCalculation()

But it Gives out 4 as parse Error instead of 0(as success)
what is the error in code ??


texts: How do I add my own data to an existing shapefile?

texts:
I want to add event data to a geographical map. I have district shapefiles for Zimbabwe that I downloaded from GIS.

However I now want to add my own data on to it. How do I do this?

I have added datafile and have a table. When I join the data to the table for the shapefile it does not show the results correctly.

However if I join the shapefile to the data table it does input correctly - but this means that I can not do any data tests as the correct table is not on the shapefile.

Sorry this is unclear I am new to GIS


texts: Separating multiple overlapping polygons using ArcPy?

texts:
How do I separate multiple overlays in a polygon using arcpy, according to image?

Each overlap must be a single polygon


texts: Path to add often used CRS to "recently used crs"

texts:
When I start a completly new QGIS I need to manually add my often CRS to have them in the list of "recently used crs". I would like to add all of them in a folder somehow and than add the path to the folder and automatically add all crs in the folder to "recently used crs".

Or maybe something like that can be done with a Python script? I'm bad with tasks related to CRS. You think something like that is possible?

I'm using QGIS 3.4.13


texts:
QTiles for QGIS v2 generated what I think are called slippy tiles. Apparently follow ons ran into complications. Tiles XYZ was "a stop-gap for QGIS 3.4 LTR users." and they said "for QGIS 3.8 and above, you do not need this plugin. The functionality is available from core processing toolbox."

Where is that funtionality? Thank you. I've searched the manual released Mar 06, 2020 and couldn't find anything.

I found it Menu:
Processing > Processing Toolbox > Raster tools > Generate YYZ Tiles


texts: Differences between tiling LAS datasets using file size versus tile dimensions?

texts:
The question refers specifically to the Tile LAS (3D Analyst) tool in ArcMap 10.7. I've tried both approaches to tiling my LAS dataset (specifying the value of the largest file vs. the dimensions of each tile).
Other than the immediately apparent differences (such as tile shape), what differences does each approach produce?


texts: QGIS WMS Layer - Returned legend image is flawed

texts:
I succeed in connecting to the following WMS in QGIS (2.14.15 LTR): https://secure.erlangen.de/arcgiser/services/Luftbild2016/MapServer/WMSServer?SERVICE=WMS&REQUEST=GetCapabilities

... with the following configuration:

The image is displayed, but in the Layer tree there is an error message Returned legend image is flawed:

It is not flawed, because legend image url http://secure.erlangen.de/arcgiser/services/Luftbild2016/MapServer/WmsServer?request=GetLegendGraphic%26version=1.3.0%26format=image/png%26layer=0 seems to work:

Since error messages like this might be confusing to simple minded people, what might be wrong here? Do I miss any settings?

UPDATE:

This seems only to occur in case of some configuration issues in our ArcGIS for Server / ArcGIS WebAdaptor environment and is not a general issue in QGIS not displaying WMS legends, as indicated by the screenshot below (WMS is http://sg.geodatenzentrum.de/wms_topplus_web_open)


texts: How to restrict the number of symbols in feature name in CartoCSS?

texts:
I make a tourist map with CartoCSS in TileMill. There are several zooms and I would like to show markers of the objects on 15 zoom and their labels (names) on 16 or 17 zoom. If the number of symbols in name is too much (for ex., more than 10) I want to label them on 17 zoom, if there are less than 10 symbols, I want to label them on 16 zoom.

Could any one help how to set the restriction on the number of symbols either in CartoCSS or sql?


texts:
My title is more than understandable,but the explanation part is that when i try to take a part that is larger than a big city i get an explanation that is to big for "trim".I tried to get whole europe in my area but I always get a info that area which I selected is to big to cut from OSM2GIS.

Is it possible to take that kind of map/area from OSM2GIS and where can I get WMS base for that area?

I hope the rest of explanation is understandable.


texts: How do I make a SHP file from waypoints?

texts:
I have a Garmin etrex 20 and I want to convert my waypoints into a shapefile for QGIS.

The waypoints need their name and background information with them in the shapefile.


texts: Simple python query: properly setting up arcpy.env.workspace and defining datasets

texts:
I'm creating a python script (for the first time in quite a long time). And the basics are confusing me a little. With the below code I keep getting the error:

"ERROR 000732: Input Table: Dataset /MARCHNET does not exist or is not supported
Failed to execute (AddField)."

In the "OK" folder there is a shapefile named MARCHNET and it is also currently open in an instance in ArcMap.


texts: How to get property of each image in a collection in google earth engine?

texts:
I am trying to convert ASTER image from DN values to reflectance. For this I require gain coefficient of each image in the collection I have gathered.

//Defining geometry
var geometry = ee.Geometry.Polygon([[72.8,29.8], [74.7,28.27], [77.6,30.1], [75.6,31.6], [72.8,29.8]]);

//Creating collection
var collection = ee.ImageCollection('ASTER/AST_L1T_003')
    .filterBounds(geometry)
    .filterMetadata('CLOUDCOVER', 'less_than', 2)
    //.filterDate('2004-01-01', '2004-03-28')
    ;
print('Collection: ', collection);
//Filtering data for specific day of year
var collection_121_134 = collection.filter(ee.Filter.dayOfYear(121,134));

Now this collection_121_134 contains 59 elements (i.e. 59 image tiles). For each tile I want to access its gain coefficient value which is found in:
Collection>>features( total no of tiles)>>any tile selected>>properties>>'GAIN_COEFFICIENT_B01'.

Please suggest a method through which I can access this information, which will help me in converting the images to their reflectance values.


texts: How To Calculate Height of Point B from Point A with X, Y, Height, Angle (Dip), Azimuth, and Depth (Distance)?

texts:
I have a survey point data which has coordinates X, Y, Height, Angle (Dip), Azimuth, and Depth (Distance).

For Example, point A:

  Easting: 290694

  Northing: 715927

  Elevation: 1060

  Angle: 65°

  Azimuth:45°

  Distance:150

Can you let me know how I can calculate the end point (End of trace) Height, please?


texts: Prevent labeling outside of a an inverted styled polygon

texts:
I have this situation:

All layers are filtered to the region except the layer "bahnstrecken" who is missing that attribute. Instead I use the polygon "bundeslaender" (inverted polygon style) to overwrite features of the layer "bahnstrecken" (black lines) outside of the shown region. This works perfect. But as soon I want to label the "bahnstrecken" (which I did with orange numbers), the labels show up outside the region. I only want to show the labels WITHIN the region.
I know the problem can be solved by just clipping the "bahnstrecken"-layer to the region. But I am just curious: Is there any other way how I can solve the problem, without manipulating the data?


texts:
When compiling mapfish I can't help but notice the quickstart instructions on the home site are not very good and have led to broken results so far. Any one able to clarify the steps? for instance when running the patch script should I save to any particular directory? the follow on command lines, should they be run in the same terminal window? as root? what are paths that need to be set on a debian system?

That and broken links and the use of $ and # in the scripts without background explanation of use in a terminal emulator and system specific intructions makes me wonder if mapfish is ready for prime time.

Thanks for any help

Here are the instructions

I can only post 2 hyperlinks

Usually if I can't build from instructions on a website it's because the documentation is not so hot, I'm not a CS major but I have a decent amount of linux experience.

http://trac.mapfish.org/trac/mapfish/wiki/HowToInstallMapFishWithPip

http://mapfish.org/doc/2.2/installation.html

Thanks Lewis

The problem is I've used both shells and I've pasted in the commands they supply (which include the # or $ symbol) I didn't use the symbol, just let the shell supply the character and the commands didn't correctly execute.

The python environment from running the download script does not have a clean way to exit the shell after the compile runs, has a bunch of > which looks like pico or some other text editor running off the shell. Makes sense as the download is text.

the (env) $ case is not returned on the next line as per the instruction set.

I've tried the virtual env as well and the problem always ends up being at the

$ source env/bin/activate

"install MapFish"

(env) $ pip install mapfish""

the first line when run does not lead to a (env) $ case being returned at the beginning of the line.


texts:
I have 17 rasters, some of which overlap one another. I wish to combine them such that where the rasters overlap the sum of their values is taken. I have tried combining them using raster calculator, simply using

raster1 + raster2 + ...

However this results in the error

"ERROR 010422: Error in porcessing the environment mask."

I have looked into this; http://pro.arcgis.com/en/pro-app/tool-reference/tool-errors-and-warnings/010001-020000/tool-errors-and-warnings-10401-10425-010422.htm

and

http://pro.arcgis.com/en/pro-app/tool-reference/environment-settings/mask.htm

All my data are in rasters and such should already have a mask enivronment as I understand it, with "NoData" being output where present.

What is the error in my approach?


texts: AttributeError: 'unicode' object has no attribute 'save'

texts:
I'm trying to multiply my raster to 100 from my ListRasters. I've adapted the code (especially the for loop portion) also from this site.
Now, I'm encountering an error at the last line, when saving the raster.

How can I successfully run this code?

    #import the module
    import arcpy
    from arcpy.sa import *
    from arcpy import env
    arcpy.CheckOutExtension("Spatial")
    env.overwriteOutput = True

    #set the workspace
    arcpy.env.workspace = r"C:\Users\Windows\Documents\JO_GIS_Analyst"

    #Create a list of raster files
    rasterList = arcpy.ListRasters("reprojected*", "TIF")

    #Loop through the list of ListRaster
    for k in rasterList:
        print k #check the presence of rasters
        OutRaster = (k)*100) #multiplying the raster by 100
        OutRaster_Name = "nofpt" + k  #specifying the output name
        OutRaster.save(OutRaster_Name) #saving the raster


texts: Show MapBasic Window equivalent in QGIS, e.g. Python Window

texts:
Is there a Show/Hide MapBasic Window (or equivalent) in QGIS such as in MapInfo?

I would like to learn SQL expressions parsing, or a Python Window viewer?


texts: Selecting features on PostGIS table from QGIS?

texts:
I use QGIS 2.14.11-Essen and PostGIS 2.1.8. This is an issue I had not with previous QGIS versions. I work with LTR version for my work so upgrade to current version (2.18.3) is not possible for now.

Selecting features tools work perfectly with shapefiles but not with tables loaded from a PostGIS database via DB Manager. Clicking on a feature (or englobing it with a selection rectangle) does not select it either visually or within attributes table. The only workaround I found is to select a feature from the attributes table (select a row) and therefore, the feature is coloring on the map canvas...

I checked this post which is not exactly the same issue (shapefiles) but the solutions didn't help me.

Eventually, I can't see any errors when I load the table or within DB Manager UI:


texts: Underlying distribution of points in an ellipse

texts:
I'm attempting to generate both a uniform distribution and a normal distribution of points that fall within an ellipse shape, made from a centre point, semi major and minor and rotation in degrees. I have code to do this (translated from MATLAB) and it plots as expected, but I was just looking for some assurance that it's doing what it is supposed to by others more knowledgeable than me. Would there be any need to account for projection issues or plotting this geographically that I'm not accounting for here?

from pyproj import transform, Proj
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns

def trans(lat, lon):
    return(transform(Proj(init='epsg:4326'), Proj(init='epsg:3857'), lon, lat))

def ellp(lat, lon, smaj, smin, rot):
    return(trans(lon, lat), smaj, smin, rot)

def dists(mean, sma, smi, a, d, norm=True):
    #https://www.quora.com/Is-it-possible-to-generate-random-points-within-a-given-rotated-ellipse-without-using-rejection-sampling
    theta = np.radians(a)

    if norm == True:
        rp = np.random.rand(1, d)
    else:
        rp = np.sqrt(np.random.rand(1, d))

    phi = np.random.rand(1,d)*2*np.pi

    x_c = rp*np.cos(phi)
    y_c = rp*np.sin(phi)

    x_e = x_c * sma
    y_e = y_c * smi

    x = np.array(x_e*np.cos(theta) - y_e*np.sin(theta)) + mean[1]
    y = np.array(x_e*np.sin(theta) + y_e*np.cos(theta)) + mean[0]

    return(y,x)

j = ellp(30, 40, 2000, 500 , 40)
f = dists(j[0], j[1], j[2], j[3], 10000, norm=False)
g = dists(j[0], j[1], j[2], j[3], 10000, norm=True)

plt.figure(figsize=(10,10))
ax = sns.scatterplot(x=f[0][0], y=f[1][0])
plt.axis('equal')
plt.show()

plt.figure(figsize=(10,10))
ax = sns.scatterplot(x=g[0][0], y=g[1][0])
plt.axis('equal')
plt.show()


texts: Make Query Layer only returns first row of each table

texts:
I have been trying to conduct an arcpy.MakeQueryLayer_management and it is only returning the first two values from each table in the SELECT.  When I run the same SELECT statement in Oracle SQL Developer it returns 39 rows in total.  Below is the sample code:

arcpy.MakeQueryLayer_management(DB, "Layer_Union", """SELECT FEAT.EVENT_ID,
SP.MEASURE AS ABS_START,
SPL.MEASURE AS ABS_END,
FEAT.TYPE_CL AS FEATURE_TYPE,
FROM TBL1 FEAT
INNER JOIN EVENT_RANGE ER ON FEAT.EVENT_ID = ER.EVENT_ID
INNER JOIN STATION_POINT SP ON ER.STATION_ID_BEGIN = SP.STATION_ID
INNER JOIN STATION_POINT SPL ON ER.STATION_ID_END = SPL.STATION_ID
INNER JOIN ROUTE RT ON SP.ROUTE_ID = RT.ROUTE_ID
INNER JOIN LINE LN ON RT.LINE_ID = LN.LINE_ID
WHERE
RT.NAME = '15'
AND
LN.DESIGNATOR IN (1,2)
UNION
SELECT FEAT1.EVENT_ID,
SP1.MEASURE AS ABS_START,
SPL1.MEASURE AS ABS_END,
FEAT1.TYPE_CL AS FEATURE_TYPE,
FROM TBL2 FEAT1
INNER JOIN EVENT_RANGE ER1 ON FEAT1.EVENT_ID = ER1.EVENT_ID
INNER JOIN STATION_POINT SP1 ON ER1.STATION_ID_BEGIN = SP1.STATION_ID
INNER JOIN STATION_POINT SPL1 ON ER1.STATION_ID_END = SPL1.STATION_ID
INNER JOIN ROUTE RT1 ON SP1.ROUTE_ID = RT1.ROUTE_ID
INNER JOIN LINE LN1 ON RT1.LINE_ID = LN1.LINE_ID
WHERE
RT1.NAME = '15'
AND
LN1.DESIGNATOR IN (1,2)
ORDER BY ABS_START""", "OBJECTID")

How do I get this query to return the same values as the SQL Developer.

Software: ArcGIS 10.4.1,
Oracle DB


texts: How do I group near points with GPS positions?

texts:
I'm an IT person so I don't know too much about projections and so on, I hope you can help me.

I have made an application for Android that collects a GPS position, so I have the latitude and the longitude at a given time. I want to save together those elements which are near to each other, in groups of terrain areas of the same phisical size; the problem is that I don't know the points beforehand and they can come from any position in the world.

My first idea (to explain a little bit more the problem) was to use the decimals of the latitude and the longitude for grouping; for instance, one group will be the positions with latitude between 35.123 and 35.124, and the longitude between 60.101 and 60.102. So if I get a position like lat=35.1235647 and lon=60.1012254598, this point will go to that group.

This solution would be OK for a cartesian 2D representation, as I would have areas of 0.001 units wide and heigh; however, as the size of 1degree of longitude is different at different latitudes, I cannot use this approach.

Any idea?


texts: How to get correct format of coordinates for export to geojson from QGIS?

texts:
I've created a custom map in QGIS with a few shapefiles. I'm using leaflet.js to render the map and the shapefile points and everything lines up correctly. The CRS is set to 4326.

When I export the shapefiles as geojson from QGIS, the coordinates are huge and look something like this: 9949941.925330460071564, -12405330.17722519300878

Leaflet is expecting coordinates that look like this: 82.9687,63.2336

Anyone have any insight into resolving this problem?


texts: General guidelines for topographic contour to DEM conversion?

texts:
It seems that many different DEM resolutions (10m, 30m, 100m) could be created from a contour map (with, for instance 50m elevation changes marked), but of course some of the resulting DEMs will be wrong in one of two ways: 1. overly 'accurate', or 2. highly inaccurate.

Are there any general guidelines when trying to produce the most accurate DEM from a contour map?


texts: Creating raster layer from data frame in R?

texts:
I have created a data frame which I need to convert to an image. This particular data frame is measured 174, 209, 36366 (nrow, ncol, ncell) respectively. This data frame is created and stored in R and I dont know how to export it from R or create a raster layer (of the same dimensions) for imaging from it.
It is effectively a table of pixel values. Can anybody help?