texts:
This Question is about Natural Earth Project. I tried to do some research about  creating my own WMS and came across this project. The web page says :

  Natural Earth is a public domain map dataset available at 1:10m,
  1:50m, and 1:110 million scales. Featuring tightly integrated vector
  and raster data, with Natural Earth you can make a variety of visually
  pleasing, well-crafted maps with cartography or GIS software.

How can I harvest the data for my WMS So that it can be used with
    Mapinfo or Other software?


texts:
Similar to Neil's question of March 8th 2013:
 QGIS not rendering WMS layers

I can down load the BGS 1:50K bedrock and superficial deposits maps OK, but they don't render except at scales greater than about 1:50,000 - this seems to vary; sometimes "only" 1:25000 renders correctly.

These scales are within the range quoted by the supplier (British Geological Survey).

Neil asked a very similar question, but there does not seem to be a definite explanation yet. I'm using QGIS 2.2.0-Valmeira.


texts: Deleting polygon defined area from raster layer in ArcGIS Desktop?

texts:
Which tool do I use in ArcGIS 10.1 for Desktop to exclude areas covered by vector polygon layers from a raster layer?

I want to keep portions of the raster layers that exclude the portions covered by the vector layer.


texts: How to change attributes with QGIS Python?

texts:
I have this piece of code which I would like to use for inserting values into a column. It does not throw an error, the changes are not written to the source (shapefile).

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

layers = iface.legendInterface().layers()

for layer in layers:
    name = layer.name()
    if name.endswith('x'):
        provider = layer.dataProvider()
        features = provider.getFeatures()
        for feature in features:
            feature.setAttribute('attr', 'a')

I understood there are two ways of editing data sources: directly with dataProvider or with edit buffer. If I got it right, it shouldn't be necessary to call commit action when using dataProvider.

I'm using QGIS 2.2 Valmiera.


texts: Using folder name of output file inside SAGA GIS Algorithm with PyQGIS?

texts:
I want to execute an algorithm using SAGA GIS. I am using the following Python code.

outputwinddata_asc = "W:/04 Tools/07 MapMaker/workingfiles/Thailand_Rectangle2_DTU50m.asc"
outputwinddata_asc1 = "W:/04 Tools/07 MapMaker/workingfiles/Thailand_Rectangle2_DTU50m_SAGA.asc"
processing.runalg("saga:reclassifygridvalues",outputwinddata_asc,0,-999.0,0,0,0,1,2,0,"0,0,0,0,0,0,0,0,0",0,True,0,False,0,outputwinddata_asc1)

Its showing (as given below) that the output file is generated successfully.

{'RESULT': 'W:/04 Tools/07 MapMaker/workingfiles/SAGA.asc'}

Unfortunately, I am not able to see any file with file name SAGA.asc in the folder W:/04 Tools/07 MapMaker/workingfiles.

Then, I tried by changing the folder of output file as given below.

outputwinddata_asc1 = "C:/Users/Sreeraj/Desktop/SAGA.asc"

Then the output file SAGA.asc is generated successfully in the folder C:/Users/Sreeraj/Desktop ; and I can see the output.

{'RESULT': 'W:/04 Tools/07 MapMaker/workingfiles/SAGA.asc'}

So, I found out that the difference is the spacing of folder names. C:/Users/Sreeraj/Desktop does not have any spaces in between folder names. But, for W:/04 Tools/07 MapMaker/workingfiles, there are spaces for 04 Tools and 07 MapMaker ( I mean 07 'space' MapMaker).

Unfortunately, I am not allowed to remove these spaces of folder names and I must keep the folder name W:/04 Tools/07 MapMaker/workingfiles as it is.

How can I solve this issue (without changing the name of output folder W:/04 Tools/07 MapMaker/workingfiles) ?


texts: Keeping OID field when exporting to CSV?

texts:
I imported a csv file and converted it to gdb format, and now I want to export it back to csv.

The file will be the exact same as before I imported it, except I want to keep the ArcMap generated OID.

However, when I export the file to csv, the OID field is lost, and contains null values only.

I understand why this is the case, because the OID is a temporary variable I guess.

Is there anything I can do to keep it though? Perhaps duplicate the field, but I am not sure how to do that in ArcMap.


texts: Running regression in each sample separately?

texts:
I'd like to smooth an NDVI image by regression analysis; response of time series NDVI image to Julian date to represent the seasonal changes in NDVI as a function of Julian day. The time series NDVI image has 26 bands (8-day images during April to October) and 12618/4144 col/row. I have chosen 30 point samples in each land cover (30 x 5 LC=150) from all band of time series NDVI image using stratified random sampling method. I need to discuss my further steps with you all.

Do I need to run the regression for each sample point (150 samples) and calculate regression equation extracted from step, 1 in each point sample (one by one) again to reproduce predicted NDVI image?

Is there any method to simplify these processes to save time?


texts: MapServer Pixmap Symbol GAP for Park Polygon Layer

texts:
I have a

Polygon layer (National Park)
PNG Pixmap symbol (from Maki icon set)
MapServer 5.6

If I use the PNG Pixmap symbol for my Polygon Layer, Mapserver fills the polygon layer with the png symbols.

I tried using

GAP 300

However it doesn't seem to work with pixmap (works fine with vector symbol)

Is there another way to use GAP for pixmap ?


texts: How to retrieve my position relative to a another on Google Maps?

texts:
I would like to SMS a URL with a call to the Google Maps API to visualize users' location relative to a location I choose, e.g.:

chosen location:
https://www.google.com/maps/place/Boston,+MA/@42.3142647,-71.110369,11z/data=!3m1!4b1!4m5!3m4!1s0x89e3652d0d3d311b:0x787cbf240162e8a0!8m2!3d42.3600825!4d-71.0588801

I am quite sure I did it once but I forgot how that works...


texts: Display connections between points in QGIS interactively

texts:
I have a layer (pointLayer) with point features:
Fields: ID, toID.
All points have a value in ID, some have a value in toID.
If a point "A" has a value in toID, I want to display a line between point "A" and the point with the ID = point "A"'s toID.
The line doesn't have to be stored permanently, but it would be good if it automatically appeared/disappeared upon a change of toID value.
How can I best achieve that?
My approach so far has been:

Creating a new GPKG layer of geometry type = Lines (lineLayer) and
adding one field: fromID.
Copying all points, changing the wkt from "PointZ" to "LineStringZ",
giving them some fake geometry and pasting them to the lineLayer.
Then updating the geometry like:

geom_from_wkt('LineStringZ(' ||
x(geometry(get_feature('pointLayer','ID',fromID))) || ' ' ||
y(geometry(get_feature('pointLayer','ID',fromID))) || ' ' ||
z(geometry(get_feature('pointLayer','ID',fromID))) || ', ' ||
x(geometry(get_feature('pointLayer','ID',attribute(get_feature('pointLayer','ID',fromID),'toID')))) || ' ' ||
y(geometry(get_feature('pointLayer','ID',attribute(get_feature('pointLayer','ID',fromID),'toID')))) || ' ' ||
z(geometry(get_feature('pointLayer','ID',attribute(get_feature('pointLayer','ID',fromID),'toID')))) ||
')'
)

But that means if I change something in the pointLayer, I have to update the geometry of the lineLayer to see the changes.
Looks like I can achieve similar result using PointConnector, but that involves generating/updating the "from - to" table and the input layer each time (because the first field must be "the joining attribute found in the from-to list").
Another idea I got was to somehow use a virtual layer for this, i.e.:
SELECT "ID", "toID", geometry
FROM pointLayer
WHERE "toID" <> ''
....
and somehow update, or create the geometry.

But my SQL knowledge is minimal.


texts: Change percentages to Real numbers in QGIS?

texts: I'm trying to convert percentages to real decimal numbers so they'll display properly. I tried  to_real() but maybe I did something wrong. Any suggestions?


texts: combing two layers dissolve and simplify

texts: I am trying to combine a dissolve with a simplify polygon layer I've created. When doing conversion of layer to kml it says I have too many vertices hence the need for simplification. How can this be done?


texts: Styling substring of label with different style without creating two separate label rules?

texts:
This question has been asked in the past but those questions were either not answered or suggested solutions that required multiple labels separated by newlines and placement rules to simulate a contiguous label.

Does the label engine in QGIS currently allow styling substrings based on rules? I'm writing a quick tutorial for a GIS class on labelling, and I'm seeing if a label like this:

"name_en" ||  '\n'  ||  'is a ' || lower("category") || ' site which is ' || case when "danger" = 1 then 'under threat.' else 'not under threat.' end

Can be modified so that the contents of the field "name_en" is displayed in bold, whereas the rest of the string is displayed using regular font. For reference, the data in question is UNESCO's list of World Heritage Sites, available here.

I can get there by splitting them up in two and devising a rule for "name_en" separately, but that feels like a wonky solution. I've tried several expressions in the 'data defined override' button, but the result always comes back as either all bold or all regular. It expects a boolean argument, so that kind of assumes the output (i.e. any given label) will always be either all bold or all regular.

I assume the answer will involve something like using the length of the name_en string to apply the bold style to the first n characters, or:

left(len(name_en))

But I'm not finding a way to use that kind of logic to style a substring of a label.


texts: Google Earth Engine script: runs forever without errors but never downloads the .tif file

texts:
I am trying to download a base map of California from Google Earth Engine using this code:

// Define boundary box
var region = ee.Geometry.Rectangle(-124.5,31.3,-110.5,39.8); // minX, minY, maxX, maxY

// using landstat data LE7_L1T
var collection = ee.ImageCollection('LE7_L1T') // dIFFERENT satellites 'LC8_L1T', 'MODIS', AST_L1T
  // collect data from these dates
  // the more data, the less cloud coverage
 .filterDate('2016-01-01', '2017-01-01')
 .filterBounds(region);

var composite = ee.Algorithms.Landsat.simpleComposite(collection, 75, 3);

// export it to drive
Export.image.toDrive({
 image: composite,
 description: 'cali-landsat-7-2016-2017',
 scale: 300, // 5000 * 4000 px aproximately, it is pretty big
 maxPixels: 270000000,
 crs: 'EPSG:32723', // projection
 region: region
});

I am trying to save the .tif file in Google Drive. Although the code has worked for me before with other areas, I do not know why is not working now.

I do not receive any error, the script runs for hours but I never get the file. Any solutions? Do you know if there are any limitations?


texts: Paste style to virtual layer from another layer

texts:
My virtual layer disappears when I paste a style from another layer. I've done this process before and it worked fine, so I'm not sure what the problem is now. Here are my steps:

create virtual layer with one field
copy style from layer

paste it to the virtual layer. You can see the symbology style is there. But, the map is blank?

I haven't changed anything else on the virtual layer after creating it. The only thing different this time was that a message popped up stating the virtual layer already exists (because I copied the style a  few times while troubleshooting), then asks if I want to overwrite it. I click yes. I'm not sure why I can't see the virtual layer (and yes, I 'zoomed to layer' and checked the CRS).
EDIT: Here is that warning message. Not sure if it has anything to do with what I'm experiencing.


texts:
I am experiencing a small screen space offset when using a WMS layer (origin is GeoServer 2.1.4):

The image shows a red square which is the top left pixel of a tile (GeoTIFF). An exact copy of the same tile is among the tiles in the background provided by my GeoServer. The four images represents the same view (the red square is the same pixel in all four images) just panned a bit inside my QGIS window. The yellow line outlines the corresponding pixel in my WMS layer - added to illustrate the misalignment. I would expect the red square to always align to the yellow line.

Can anyone provide me with some knowledge to why I'm having this misalignment depending on the extent shown on my screen?

I have posted a similar question here proposing that the misalignment could be from a difference in the aspect ratio between the width/height and the bbox in the WMS request. But I'm not sure this is the reason and thus the reason for posting two different questions on the same topic.


texts: I am building an application that requires me to plot UTM coordinates on a chart using CSV data. I am given two UTM coordinates for each record, although I would like the application to handle zone crossing with little loss of accuracy. I am not given a UTM zone, only the coordinates. Is it possible to handle UTM zone crossing without the information of which zone you are in given zone? Only given coordinates such as 539594.4,9071398.0 (probably not outside zone)?


texts: ArcGIS 10.2 Modelbuilder Spatial Join is not generating an output

texts: I have a model that runs spatial joins in a lines to points relationship and then adds and calculates fields to that join. The parameters are one to one and intersecting at 1 foot.  The search distance is so low because the points are already snapped to the lines. During my initial test run, the spatial joins worked fine in the model, but I realize that the difference here is that I was working in a local environment. This time, I am required to run the model using inputs from a network drive. Could this be my issue? Or am I missing something really obvious here? The spatial join is on a dataset of at least 20,000 lines, so it should be taking around 20 minutes to complete.  IN this new model, the operation is done in mere seconds, so I know its not doing what its suppose to be doing.


texts: Progress bar percentage is increasing and decreasing, never finishing

texts:
I am trying to clip a raster by shapefile perimeter using the Clip tool. The percentage has been jumping between 0% and 100% for over an hour, back and forth, climbing to 100% then dropping down to 35% or so all over again. Does anyone know how to fix this? I'm desperately trying to cut down this raster file so I can continue doing my work, since it has to re-draw the entire raster every time I add a point to a point shapefile.

I'm using ArcGIS10 with SP4 installed. This occurred in ArcMap while running the Clip tool under Data Management Tools > Raster > Raster Processing > Clip.


texts:
I want to remove some nodes from a shapefile using QGIS 2.10.1

I enable editing
I click on the Node Tool
I click on a node I want to select (it turns blue)

However, once I select this node I am unable to delete it.
The Delete Selected button is grayed out, as is cut, and pressing delete on my keyboard does nothing. I assume I am missing something very basic but I can't figure out what it is.


texts:
In a project using gdal I need to crop a geotiff file with a bounding box in pixel or coordinates epsg:4326. The cropping of large images works fine.

The next step is to manage cases where the image is smaller than the bounding box. The idea here is to avoid getting blank data.

When I crop my geotiff using the gdal_translate tools with -srcwin or -projwin an error occurs when the size of the bounding box is bigger than the geotiff size.

When I crop my geotiff using the gdalwarp tools with -te xmin ymin xmax ymax, my result get the size of the bounding box even if this one is bigger.

What could be the "average" of those two behaviours? Just to fit my result to the interception of the bbox and the geotiff file?

My process makes part of a .NET Project so I could do it into my code using .NET but is there a way or an option in gdal to do it?

I found documentation about the ogr2ogr command -clipsrc but it doesn't manage geotiff file.


texts: Is the ArcGIS Server JavaScript API free to use?

texts:
The 2008 blog post announcing the ArcGIS Server JS API says:

  Everyone can use the ArcGIS JavaScript API. There is no fee for using the ArcGIS JavaScript API or deploying an application built with it.

But Esri's Web Mapping Pricing page says "Contact us for pricing".

Which is correct?


texts:
I want to shorten a feature (a polyline representing a road) in QGIS 2.4

In QGIS 2.0 there was "Split Features" tool in the "Advanced Digitizing" toolbar.

Here the icon:

I'm not able to do this again.
Thanks.


texts: Where to get city boundaries without sea area as a geometry polygon?

texts:
I am using OpenStreetMap's Nominatim to find out boundary geometry of single cities: Helsinki
That polygon however contains the sea area and I would like to exclude that.
So where could I get city boundaries containing only land area or coastline geometry so I could exclude sea area from the city boundary polygon.


texts: I have OpenGeo Suite installed on a server. I want to have the PostGIS data stored on a different drive (D) on the server instead of on the C drive under my username. How do I do this?


texts: Align the border of a polygon to overlap the border of another using QGis

texts: I have two polygon layers. As seen in picture, the border of two layers do not fit. I want to align the border of layer 1 (blue) to overlap the border of layers 2 (pink). How to do it with QGIS?


texts: Find maximum and minimum value in a list of rasters

texts:
I trying to write a standardisation script that standardises each raster dataset in a list of rasters by the highest and lowest values found in the list of rasters.

The standardisation code is as follows:

for raster in rasters:
  output = (Raster(raster) - min_value) *100 / (max_value - min_value) + 0

Where the max_value and min_value are the maximum and minimum values in the raster list.

I initially though I would be able to use the GetRasterProperties_management function on the raster list but this doesn't work because the raster list is just a list!

The other idea I had was to loop through each raster in the list getting the max and min values and then outputting them as a table of values. Then finding the max and min value in each table to use in the standardisation equation. But I'm struggling to figure out how to output max and min values as a table.

So I would like help if anyone has got a better idea for the whole procedure or knows how to output max and min values from a looped GetRasterProperties_management into a table.

I'm using ArcGIS 10.


texts:
I have the following map:
https://aressk.cartodb.com/viz/408dd3ea-deea-11e4-9d17-0e9d821ea90d/embed_map

It has two layers.

On the Waypoints layer, I have Hover turned on for the Info Window.
On the SPOT Tracker layer, I have Click turned on for the Info Window.

The hover for the Waypoints works, but the problem is that if I hover over points from the SPOT Tracker layer I get a "No data available" message that I don't want.

If I turn off Hover on the Waypoints layer, the "No data available" will not appear when hovering over a point from the SPOT layer.

Any suggestions on how I can remove the "No data available" message?


texts: Coordinates freeze when I give a projection

texts:
I have recently begun to use ArcMap and needed some guidance on a problem that I am facing.

I have a district shapefile with 200 districts with geographic coordinate system: GCS_WGS_1984. I need to calculate the area of each district for which I am trying to project it with Asia Lambert Conformal Conic. But when I do this, the coordinate system displayed just stand froze i.e. it displays 95.001 30 decimal degree no matter where I put the cursor in the display.

Is there something wrong with the shapefile I have?


texts: Zonal histogram for each region (row) in a polygon (table) with python script

texts:
I have DEM.tif (symbology uses equale interval - 100m) and a polygon file (regions.shp) with many regions (names are defined in 'reg_name' column). I would like to use Zonal Histogram to produce tables (and later xlxs files) with two columns: first one should be class (200-300,300-400 etc.) and the second one the frequency.
I have DEM.tif and regions.shp in my map document and if I use code in python window for one region (without using a for loop) it works. But when I want to run the code to generate tables for each region it doesn't work!
The main problem is probably because after the lineextracted_raster=arcpy.sa.ExtractByMask(in_raster, lyr) the 'extracted_raster' does not actually show in table of contents or in data view. Thus arcpy.ApplySymbologyFromLayer_management("extracted_raster",in_raster) cannot be applied.
I noticed that if I run the code only up to arcpy.ApplySymbologyFromLayer_management("extracted_raster",in_raster) the 'extracted_raster' is shown in my TOC only after the last one iteration (the last region in the list). All 'extracted_raster' rasters before the last region are not shown... How can I make them appear? I'm using arcpy.RefreshActiveView() but it doesn't show the raster...
And I know that Zonal Histogram produces table with columns for all regions at once, but I'd like to get tables with actuall min and max in a table for each region (without classes with frequency zero before and after).
Here is my script:

#Set MapDocument to "CURRENT"...
mxd=arcpy.mapping.MapDocument("CURRENT")

#Defining parameters...
in_raster='DEM.tif' #Raster to be extracted and used for histogram data production...
zone_layer='regions.shp' #This is a layer with regions...
zone_field='reg_name' #The name of the field with names of the regions...

#Code for listing names of regions...
def getValueList (inputTable, field):
    values = set()
    rows = arcpy.SearchCursor(inputTable)
    for row in rows:
            values.add(row.getValue(field))
    return sorted(values)
arcpy.AddMessage("Koda getValuesList je definirana!")

#Listing names of regions...
list=getValueList(zone_layer, zone_field)
lyr=arcpy.mapping.Layer(zone_layer)

where_clause=(str(zone_field))+"="

#Loop for
for zone in list:
    lyr.definitionQuery=where_clause+"'"+str(zone)+"'"
    arcpy.RefreshActiveView()
    extracted_raster=arcpy.sa.ExtractByMask(in_raster, lyr)
    arcpy.RefreshActiveView()
    arcpy.ApplySymbologyFromLayer_management("extracted_raster",in_raster)
    output_table="histogram_data"+str(zone)+".dbf"
    arcpy.sa.ZonalHistogram(lyr,zone_field,"extracted_raster",output_table)
print "Done!

"


texts: Con expression in two rasters ArcGIS, raster calculator

texts:
I would like to write a condition in ArcGIS 10.1 in Raster calculator.

I have 2 overlapping rasters (Ras1.img, Ras2.img) with different attributes in attribute table:

Ras1.img: Class - results of supervised maximum likelihood classification
Ras2.img: Att   - landuse, raster converted from shapefile polygon format

I want to change pixel values Class in Ras1.img based on the values in Att, Ras2.img. something like:

if ((Ras1.img, Class == 20 & Ras2.img, Att == 5), then change Class value 20 in Ras1.img  to Class 25))

I can't figure that out. It is super simple in R, but I have no idea how to make it in ArcGIS 10.1.

I was trying to figure that out by : http://resources.arcgis.com/en/help/main/10.1/index.html#//009z000000zn000000

and applying something like:

Con(("landuse1" == 1) & ("landuse2" == 5), "landuse1" + "landuse2", 99)

but I don't understand where I specify the attributes of my rasters I would like to change?


texts: Database versus Files for Solo Developer

texts:
I'm just getting started in GIS and am interested in using it to create maps of mountain biking trails in my area and nearby. I've been learning QGIS and have made enough progress to generate a decent map. So far I've been storing my layers in Shapefiles, which seems to be working okay.

This past weekend I started looking at PostGIS and have it up and running and managed to import my trail shapefiles and connect QGIS to PostGIS.

I'm now wondering what the pros and cons of using a database versus files would be. I have some thoughts already, but am curious what enlightening opinions others may have.

Pros of Database

One table for all my trails, regardless of park, allows me to have a consistent set of columns for the information.
Easier to query/export my data to multiple formats
Multiple clients/tools allow for easier viewing and editing of data
Easier editing of table structure (compared to QGIS table editing)

Cons of Database

Standing up/maintaining a database adds complexity to my setup.
I have a desktop and a laptop I do my GIS work on and using DropBox allows me to easily sync files between the two machines. This likely won't be possible (or as transparent) using PostGIS as using DropBox with files.
Easier to backup/restore files than a database.


texts: How to insert a single point feature into a shapefile?

texts:
I have the following problem. I am doing an analysis which creates a shapefile "A" with one or more point features. In the next step I need to insert the individual point features of the shapefile "A" into another point shapefile "B" which stores the results from multiple analyses. How do I insert (append) features from "A" into "B"? Both shapefiles have the same attributes (shapefile "B" has one additional attribute) and I need to transfer them as well.

Here is the code I think should work but it does not:

analysis_res = fcn1(inp1, inp2) # analysis results - shapefile "A"
cursorA = arcpy.da.SearchCursor(analysis_res)
cursorB = arcpy.da.InsertCursor(outFile, ['SHAPE@'])
for row in cursorA:
     cursorB.inserRow(row) # here I would like to insert the poit feature of "A" into "B"
del cursorA, cursorB


texts: Update attribute table by values from other layer

texts:
let's imagine that I have 2 different layers with geometry. Some objects from the first layer are "within" the second layer's objects. I'd like to update firstLayer_value in attribute table by value of attribute secondLayer_value from second layer for objects which are within. Till now I was using refFunctions plugin by creating new column with query: geomwithin('layerName','columnName') but it works super slow for many objects. Is there a faster way to do this?
My layers are in .shp format, In main line layer i got vectors which represents dykes or cuttings (2100 records) which goes through the roads, and second layer is 2m buffors of roads(22000 records). Now I need to create new column on my dykes layer which have value of the road's width. Its stored in attribute width on my buffors layer.


texts: Data security for GIS Data (shapefiles,AutoCAD files)

texts:
This is regarding security for GIS Data (in Esri & open source) (Esri shapefiles, AutoCAD files etc.)

Does anyone know how to assign a password to or encrypt shapefiles OR CAD Data?

I have found an ESRI script on this site but I think its too old..So if anyone knows such information/tool then please let me know.


texts: Multiple labels for only selected features of one layer and multiple colors

texts:
I have a point layer in my QGIS project that displays a big number of addresses of a number of cities.
In the map for each feature, I would like to display two values: The housenumber and the number of units currently residing in that address. Therefore, I have two fields: Housenumber and Units.

I want the housenumber to be displayed above the point on the map in blue and the number of units below in red, similar to this:

Currently, I have solved this issue by the way that is explained in this thread.

That is a pretty good solution, but unfortunately using those text diagrams does not seem to allow using Symbology.
Now, I would like to use a rule-based symbology, to be able to create copies of the layer that each display only the addresses for one city.
Normally, I would just use the expression "City" = '...' and assign my symbology to it, and everything would be perfect.
But this does not seem to work in combination with the text diagrams, because no matter what rule I put in the symbology, it still displays me all features.

The only solution I could think of was creating two copies of the layer per city, and configure the symbology to once display the housenumbers and once the units, but I was wondering whether there is a more elegant way?

Thanks a lot!


texts: List for ui.Select from Modis land cover classes in Earth Engine

texts:
Previously, I have asked a question and received very valuable answers and useful solutions. But, still I could not get what I want. I really need support to solve the issue. I am still struggling on making a list of land cover classes in particular study area for ui.Select?

My question was: "I would like to identify the list of land cover classes in particular area and use them as an item for ui.Select. And I would like to use below given label for identified land cover values. I have drafted a code, but I am facing with an error. Is anyone can help me to improve it?"

Link for the question: List of land cover classes in particular area and use them as an item for ui.Select

var block =
    ee.Geometry.Polygon(
        [[[10.884499359130846, 50.55773169053401],
          [10.884499359130846, 49.37053314015084],
          [13.433327484130846, 49.37053314015084],
          [13.433327484130846, 50.55773169053401]]]);

var dataset = ee.ImageCollection('MODIS/006/MCD12Q1')
  .select('LC_Type1')
  .filterDate('2018-01-01','2018-05-01');

var landcover = dataset.first();
var block_landcover = landcover.clip(block);
var scale = block_landcover.projection().nominalScale();

print ("MODIS Scale", scale);

var Vis = { min: 1.0, max: 17.0,
  palette: [
    '05450a', '086a10', '54a708', '78d203', '009900', 'c6b044', 'dcd159',
    'dade48', 'fbff13', 'b6ff05', '27ff87', 'c24f44', 'a5a5a5', 'ff6d4c',
    '69fff8', 'f9ffa4', '1c0dff']};

Map.centerObject(block_landcover, 8);
Map.addLayer(block_landcover, Vis, 'Land cover');

var classes_list = ['Evergreen Needleleaf Forests', 'Evergreen Broadleaf Forests',
  'Deciduous Needleleaf Forests', 'Deciduous Broadleaf Forests',
  'Mixed Forests', 'Closed Shrublands', 'Open Shrublands', 'Woody Savannas',
  'Savannas', 'Grasslands', 'Permanent Wetlands', 'Croplands', 'Urban and Built-up Lands',
  'Cropland, Natural Vegetation, Mosaic', 'Permanent Snow and Ice', 'Barren', 'Water Bodies'];

var list = ee.List.sequence(1, 17, 1);

var areas = list.map(function(l){

    var iter = ee.Number(l).int();

    var areas_list = ee.List([]);

    var areas = block_landcover
     .eq(iter)
     .multiply(ee.Image.pixelArea().divide(10000))
     .reduceRegion({
        reducer: ee.Reducer.sum(),
        geometry: block,  // a geometry
        scale: 463.3127165275,   // scale for MODIS
        maxPixels: 1e9
    });

    areas_list = areas_list.add(areas.values()).get(0);

    return areas_list;

}).flatten();

var fus_list2 = ee.List(classes_list).zip(areas);

print(fus_list2);

Link for the solution: https://code.earthengine.google.com/3d42306c59c4c3df5794e1001a51f305


texts: What does "unit name Degrees not found" error mean in ERDAS?

texts:
I want to create a mosaic in Erdas. When I want to add an image file in the Mosaic Tool interface in Erdas, Erdas showed me the error saying "unit name Degrees not found". Then I tried to change the map projection from Geographic Lat/long to UTM. But this error still appears. What does this error mean?

Thank you


texts: Polygons split themselves when I zoom out

texts:
I've been having issues with polygons getting deleted and re-created in this layer, and also with polygons not merging the way I would expect. I thought I had everything fixed-- my topology check comes back clean (no gaps or overlaps) and when I fill the layer with a solid color there are no gaps-- but now my layers are behaving very strangely, and I'm worried I might not have fixed everything and I don't know what sanity check to use aside from my topology.

For instance, this polygon was originally split in two along a vertical line:

I merged it:

but sometimes when I zoom out, it looks like it is still in 2 pieces.

When I export it to a KML file and open it in Google Earth, it still seems to be in 2 pieces:

I feel like I'm pretty experienced with ArcMap but I have no idea where to start with trying to fix this!

I am using ArcGIS 10.6.1 on Windows 7.


texts: NO2 CSV downloaded files shows zero values for consider date range using Google Earth Engine

texts:
When I am downloading the CSV file, the results show zeros in the CSV file. How can I download exact NO2 values as shown in Chart after running the mentioned below code in Google Earth Engine?
// Create a FeatureCollection from the list and print it.
var BergamoGeometry = ee.FeatureCollection(BergamoGeometry);
print(BergamoGeometry);
Map.centerObject(BergamoGeometry);
Map.addLayer(BergamoGeometry);
//////////////////////////////////////

var dataset = ee.ImageCollection("COPERNICUS/S5P/NRTI/L3_NO2")
                 .filterBounds(BergamoGeometry)
                  .filter(ee.Filter.date('2019-03-30', '2019-04-29'));

var NO2 = dataset.select('NO2_column_number_density');
var band_viz = {
  min: 0,
  max: 0.0002,
  palette: ['black', 'blue', 'purple', 'cyan', 'green', 'yellow', 'red']
};
// Create and print the chart.
print(ui.Chart.image.series(NO2, BergamoGeometry, ee.Reducer.mean(), 30));


texts: Batch process for merging features in multiple layers

texts:
I have multiple (50) layers representing rivers and creeks which are polylines. However, they are broken into line segments so each river/creek layer might have 30 or so features within it.

I know how to select and merge features within 1 layer but is there a way in QGIS to do this as a batch process with all 50 layers?

Otherwise, I have to manually go through and merge each line layer separately.

Perhaps someones know of a script to perform this.
The process for 1 layer is:

open layer for editing
open attribute table
select all in table
merge selected features
close and save

Essentially I need to run these 5 processes on each of the 50 layers?


texts: Added new record to attribute table, but can't display it

texts:
I used the Editor tool to add a new record to the existing attribute table of my shapefile.

But when I save my edits, it doesn't appear on my map. In theory, it should be a point located to the upper left of Conway.

I read this answer that told me to calculate the geometry for my coordinate fields, but when I do that, I just get "null". What am I doing wrong...?


texts: How do I change from lat/long to a metric coordinate system in GRASS?

texts:
Yes, it's yet another transformation question.

I am working with ASTER DEM data, which is LatLong WGS84. I want to process the DEM map into an STL file for use in a CFD simulation. The script I am using to accomplish this requires metric coordinates, though.

Here's what my ASTER DEM mosaic looks like:

I have created a new location in GRASS and chosen UTM with the etrs89 datum and grs80 spheroid. (I don't know if this is the correct choice.) Following these instructions, I created a mosaic of some ASTER tiles using gdalwarp:

gdalwarp -t_srs "+proj=utm +zone=32 +datum=etrs89" -r bilinear -dstnodata -9999 -co "TILED=YES" AST_DEM.vrt AST_DEM_Rhein_Main.tif

and managed to successfully import them.

(EDIT: It would appear that proj4 does not support ETRS89 as a datum:

$ proj -ld
__datum_id__ __ellipse___ __definition/comments______________________________
       WGS84 WGS84        towgs84=0,0,0
      GGRS87 GRS80        towgs84=-199.87,74.79,246.62
                          Greek_Geodetic_Reference_System_1987
       NAD83 GRS80        towgs84=0,0,0
                          North_American_Datum_1983
       NAD27 clrk66       nadgrids=conus,ntv1_can.dat
                          North_American_Datum_1927
     potsdam bessel       towgs84=606.0,23.0,413.0
                          Potsdam Rauenberg 1950 DHDN
    carthage clark80      towgs84=-263.0,6.0,431.0
                          Carthage 1934 Tunisia
hermannskogel bessel       towgs84=653.0,-212.0,449.0
                          Hermannskogel
       ire65 modif_airy   towgs84=482.530,-130.596,564.557,-1.042,-0.214,-0.631,8.15
                          Ireland 1965

Apparently proj doesn't fail with an error if it doesn't recognize a parameter, it just silently proceeds without warning the user, even if the result is inaccurate. Caveat emptor!)

I have a couple of problems with the result. First, when I imported the resulting mosaic, GRASS complained that it didn't recognized the datum:

r.in.gdal input=/home/sfbosch/Documents/Quelle/OpenFOAM/GIS/ASTER/Data/AST_DEM_Rhein_Main.tif output=AST_DEM_Rhein_Main --overwrite
WARNING: Datum <unknown> not recognised by GRASS and no parameters found
Projection of input dataset and current location appear to match
r.in.gdal complete. Raster map <AST_DEM_Rhein_Main> created.

Next, the mosaic obviously looks quite, well, warped:

I expected some change with the reprojection, but nothing that stark; could this be due to the non-existent (in proj4) datum? I am having trouble recognizing the topographical landmarks I had picked out for myself. Finally, I had been using OpenStreetMap data, and now it won't import because the projections don't match:

ERROR: Projection of dataset does not appear to match current location.

GRASS LOCATION PROJ_INFO is:
name: Universal Transverse Mercator
proj: utm
zone: 32
no_defs: defined
datum: etrs89
ellps: grs80
towgs84: 0.000,0.000,0.000

Import dataset PROJ_INFO is:
name: Lat/Lon
proj: ll
datum: wgs84
ellps: wgs84
no_defs: defined

You can use the -o flag to v.in.ogr to override this projection check.
Consider generating a new location with 'location' parameter from input data set.

I tried it with the override, but ended up with a layer that wasn't visible. I've since learned that I also need to reproject the vector layers with v.proj.

Really, all I want is to work with metric coordinates; it almost doesn't matter which system, the topography is the important part. Was there something wrong with my gdalwarp approach? Am I making this too complicated? What's the best way to accomplish what I want?


texts:
I have a .ecw file that I want to convert to a GeoTIFF. The CRS is EPSG:2052, but it needs to be ESRI:102566.

Essentially both the X and Y value signs need to be multiplied by -1.
How do I go about fixing this?

Thanks
earthworkshop


texts: Calculating area of blocks in polygon intersect?

texts:
I need help with calculating the area of blocks in a circular polygon intercept. The circular polygon intercept is 2000m and arable field blocks are present within the intecept. In some of the cases these field blocks are being cut and I would therefore need to know how much of the area of each field block that is inside the polygon. I would like to know this for all the field blocks that are present within the circular polygon intercept.

I have a CSV file that tells me the area of all the arable field blocks present in my shapefile and I have made a circular polygon intercept but I'm unable to make ArcGIS calculate the area of each block inside the polygon, especially the ones being cut.

The intercept was perfomed by creating a buffer and then intercepting the buffer.

Here's a picture of the polygon and the arable field blocks in green and yellow. They are separate colour because of interest in a specific crop - yellow one.


texts:
I have a table of boundaries(polygons) which need to be served up through ArcGIS Server with their centroids(points). I'm using Oracle spatial and thought the most efficient way of doing this was to create a spatial view where the SDO_CENTROID function was used on the geometry.

The table definition:

create or replace view VW_LOCALITIES_CENTROIDS as
select
    l.IDENTIFIER,
    l.SUBURB,
    l.LOCALITY_NAME,
    l.LOCALITY_CODE,
    l.LOC_STATUS,
    SDO_GEOM.SDO_CENTROID(l.GEOMETRY, (SELECT DIMINFO FROM USER_SDO_GEOM_METADATA WHERE TABLE_NAME = 'VW_LOCALITIES_CENTROIDS')) GEOMETRY
FROM LOCALITIES l;

I have a metadata entry which looks like:
table_name: VW_DCDB_LOCALITIES_CENTROIDS
column_name: GEOMETRY
diminfo: ((X, -180, 180, 0.05),(Y,-90,90,0.05)(,,,)(,,,))
SRID: 8311

The coordinate system is GDA94.

I add this view to my ArcMap and it won't draw. I then try to browse the table and the error I get is ORA-13226: interface not supported without a spatial index. I can't put a spatial index on the view or the Oracle Spatial function as far as I know.

Does anyone know how I can publish a working service with the view I created?

I would like to avoid creating a materialized view..

There are 110 rows in the base table and in the view.

The purpose of this service is for someone to be able to query the attribute "SUBURB" and get back the centroid of its polygon. They are going to attach the centroid to multiple records. The service will NEVER be used to draw geometry.


texts: Auto Trace Installation Error in QGIS 3.8

texts: My QGIS 3.8 doesn't support Auto Trace 1.0.6 or 1.0.7. My pc has a 64bit Windows 10. How can I solve it?


texts: Auto incrementing new column using Field Calculator in QGIS modeler?

texts:
I am making a model which requires creating a grid and clipping various layers to it. Using the Create grid function does not provide an ID column, instead it creates 2 columns: "longitude" and "latitude".

Because of this, I would like the Modeler to insert a new ID column using either Field Calculator or Advanced Python Field Calculator but I have no idea what the expression is.

I am using QGIS 2.2.0 Valmiera.

I have noticed this topic: Auto-incrementing in Python Script with ArcPy cursor?

If using the Advanced Python Field Calculator, would the python coding be inserted to the Global expression or to the Formula?


texts: ECW Support removed QGIS - Can't reinstall QGIS on Ubuntu 18.04

texts:
I had QGIS 3.10 installed on Ubuntu 18.04.  Hoping for ECW support I ran the bash code located

here.  While this code executed without error QGIS 3.10 is not available on my computer anymore.

After running the installation again using

sudo apt-get install qgis python3-qgis qgis-plugin-grass

I get the following error:

The following packages have unmet dependencies:
 python3-qgis : Depends: python3-qgis-common (= 1:3.4.15+28bionic-ubuntugis) but it is not going to be installed
 qgis : Depends: qgis-providers (= 1:3.4.15+28bionic-ubuntugis) but it is not going to be installed
        Depends: qgis-common (= 1:3.4.15+28bionic-ubuntugis) but 1:3.10.2+28bionic-ubuntugis is to be installed
 qgis-plugin-grass : Depends: qgis-plugin-grass-common (= 1:3.4.15+28bionic-ubuntugis) but 1:3.10.2+28bionic-ubuntugis is to be installed
E: Unable to correct problems, you have held broken packages.

How do I install QGIS given this situation?


texts:
I am a beginner to QGIS. I have found the QEPANET plugin to be very interesting  with a lot of features. While locating a reservoir or tank there are no problems. When working with pipes, it works well but on saving, I receive the following message: "Could not commit changes to layer pipes" and "geometry type is not compatible with the current layer."
The pipe layer was created by QEPANET. What can I do?


texts:
I received a folder containing different .adf files ranging in size. My task is to turn this into a TIN file which then I can bring into Microstation.  I have searched the web for days and it all leads back to using 3D Analyst in ArcMap, something I do not have access to.  My ArcMap skills aren't the greatest as I was just thrown onto this program.

Is there any other way to create a TIN from these .adf raster files?

Does Bentley Map do such a thing?


texts:
Trying to create 9 separate leaflet legend variables for the 9 different overlays (graduated colour chloropleths) I have on a leaflet map. Each overlay will have 6 classifications of colour dependant on the attribute being visualised by the bounds of colour differ between them - for instance on 1 overlay a feature might have a value of 100 and will be complete black whereas on another overlay a value of 100 might be light grey. Once I know the structure for creating 1 legend I can then replicate it appropriate for the other overlays.

The legends I have seen so far all make use of a getcolors() function that is also used in the styling of the layer. I just want to form the table of the legend so i can plugin all of the necessary colours and bounds.

legend.onAdd = function (map) {

var div = L.DomUtil.create('div', 'info legend'),
    grades = [0, 10, 20, 50, 100, 200, 500, 1000],
    labels = [];

// loop through our density intervals and generate a label with a colored square for each interval
for (var i = 0; i < grades.length; i++) {
    div.innerHTML +=
        ' ' +
        grades[i] + (grades[i + 1] ? '&ndash;' + grades[i + 1] + '
' : '+');
}

return div;
};

legend.addTo(map);

Instead of the method above, how would I create a table with 6 colour boxes and 6 corresponding labels?

The layers are being style with their own style function that is called on when they're added:

        function doStyleGrandtotal(feature) {
    if (feature.properties.grandtotal >= 0.0 &&
            feature.properties.grandtotal <= 500.0) {

        return {
            color: '#000000',
            weight: '1.04',
            dashArray: '',
            fillColor: '#2b83ba',
            opacity: '0.6',
            fillOpacity: '0.6',
        }
    }
    if (feature.properties.grandtotal >= 500.0 &&
            feature.properties.grandtotal <= 1000.0) {

        return {
            color: '#000000',
            weight: '1.04',
            dashArray: '',
            fillColor: '#80bfab',
            opacity: '0.6',
            fillOpacity: '0.6',
        }
    }
    if (feature.properties.grandtotal >= 1000.0 &&
            feature.properties.grandtotal <= 1500.0) {

        return {
            color: '#000000',
            weight: '1.04',
            dashArray: '',
            fillColor: '#c7e8ad',
            opacity: '0.6',
            fillOpacity: '0.6',
        }
    }
    if (feature.properties.grandtotal >= 1500.0 &&
            feature.properties.grandtotal <= 2000.0) {

        return {
            color: '#000000',
            weight: '1.04',
            dashArray: '',
            fillColor: '#ffffbf',
            opacity: '0.6',
            fillOpacity: '0.6',
        }
    }
    if (feature.properties.grandtotal >= 2000.0 &&
            feature.properties.grandtotal <= 2500.0) {

        return {
            color: '#000000',
            weight: '1.04',
            dashArray: '',
            fillColor: '#fdc980',
            opacity: '0.6',
            fillOpacity: '0.6',
        }
    }
    if (feature.properties.grandtotal >= 2500.0 &&
            feature.properties.grandtotal <= 3000.0) {

        return {
            color: '#000000',
            weight: '1.04',
            dashArray: '',
            fillColor: '#f07c4a',
            opacity: '0.6',
            fillOpacity: '0.6',
        }
    }
    if (feature.properties.grandtotal >= 3000.0 &&
            feature.properties.grandtotal <= 3017.0) {

        return {
            color: '#000000',
            weight: '1.04',
            dashArray: '',
            fillColor: '#d7191c',
            opacity: '0.6',
            fillOpacity: '0.6',
        }
    }
    }


texts:
I'm looking for the simplest way to track changes ([delete, update, create] feature(s)) made to the DEFAULT version of a given feature class in ArcSDE(10) and Oracle(11g).

Qualifier: I am not a dba and I have only enough of an understanding of versions to get by as a self-proclaimed developer.

I'm trying to avoid:

Mucking around with state IDs, delta tables, and triggers if poss.
Running expensive feature-by-feature comparisons with target datasets.

I'm hoping to access this 'change log' table/query with python/cxOracle or .Net if necessary.

Update:
*Here are some additional items I'd like to avoid:*

Developing, deploying, and managing custom editing class extensions to all editing users' machines.
Upgrading to ArcGIS 10.1


texts: plot coordinates as points with spplot in R

texts:
I want to do something that I thought should be relatively simple, but it doesn't seem to work...

I have a SpatialPointsDataframe and a SpatialPolygonsDataframe.

I want to plot all the points I have in the SpatialPointsDataframe together with the lines I have in the SpatialPolygonsDataframe. The name of my SpatialPointsDataframe is UK.points and the name of my SpatialPolygonsDataframe is UK.map. They both use the same OSGB projection. I try the following:

# create outline from the UK.map object
UK.outline <- list("sp.lines", UK.map, col = "black")

# create a variable to plot (all the same value)
UK.points$var1 <- 1

# plot
spplot(UK.points, "var1", cex=0.2, sp.layout = UK.outline, par.settings =
     list(axis.line = list(col =  'transparent')))

This gives me the error:

Error in cut.default(as.matrix(z), cuts, dig.lab = 4, include.lowest = TRUE) :
  'breaks' are not unique

If var1 has 2 different values, than it does work, but it plots the two different values in 2 different colors. I just want points of all coordinates all black. How to do this?


texts: How to set paint effects programmaticaly

texts:
In my plugin I want to set up a 'halo'-effect like this:

So far I tried this:

setting up the paint effect:

    QgsPaintEffect * p_mPaintEffect(0);
    p_mPaintEffect= new QgsOuterGlowEffect();

    (static_cast<QgsGlowEffect*>(p_mPaintEffect))->setSpreadUnit(QgsSymbolV2::OutputUnit::MapUnit);
    (static_cast<QgsGlowEffect*>(p_mPaintEffect))->setSpread(2.0);
    (static_cast<QgsGlowEffect*>(p_mPaintEffect))->setBlurLevel(1);

assign it with the symbollayer:

p_Symbol->symbolLayer(i)->setPaintEffect(p_mPaintEffect);

After opening the layer properties dialog there is only one painteffect not like in the image above a list of 5. But I don't know how to set more effects. Where can I get the list? In API is only documented to set one effect.


texts: MapInfo/Mapbasic - test if a mapper is open

texts:
I have a custom menu that I created which opens a mappable table and I'm trying to work out how I open it when

there is no other .tab files or mapper windows open
there are other .tab files open which are in a mapper.

My question relates to the syntax of how I test if there is a mapper open. I then need to use the results of the test to work out whether I use 'map from' or 'add map'.


texts: How to remove attribute fields with _1 or _2 suffixes from shapefiles which have be created due to merge process?

texts:
I have a large number of archive shapefiles which need their attribute fields tidying. During their history they have gone through various merge processes, which have resulted in some of the shapefiles having duplicate attribute fields where the system has auto assigned the suffix '_1', '_2' etc to the attribute field name. The attribute fields vary between shapefiles, and so can the duplicated attributes. Therefore I am thinking a wildcard solution rather than a definitive list removal solution is the way to go.
I am looking for an arcpy script which will run through multiple shapefiles from input source folder, remove all duplicate attribute fields suffixed with _1, _2, _3 etc and save the updated shapefiles to an output folder.
I have managed to perform the required task on a single shapefile using the acrpy code below. I use the "1" as a wildcard which seems to work in removing all attribute fields with suffix of _1, and the same  for *_2 and *_3 attributes:
import arcpy

fc = 'd:\\Input\\Buildings.shp

One = [f.name for f in arcpy.ListFields(fc,"*1*")]
arcpy.DeleteField_management(fc, One)

Two = [f.name for f in arcpy.ListFields(fc,"*2*")]
arcpy.DeleteField_management(fc, Two)

Three = [f.name for f in arcpy.ListFields(fc,"*3*")]
arcpy.DeleteField_management(fc, Three)

However, I understand this is very clumsy and it falls over if the input shapefile does not contain any  attribute filed that contains an underscore digit suffix. So I am looking for a more elegant solution for removing these duplicated attribute fields using ArcPy.


texts: Trouble with if/raise SyntaxError and ArcPy?

texts:
I am having trouble because when I enter a get parameter it takes me to the correct city name but it displays my SyntaxError also, but I want the SyntaxError to just be displayed when a city name is misspelled or not within the attribute table. Any help?

import arcpy
from arcpy import mapping
#set workspace
mxd = arcpy.mapping.MapDocument("CURRENT")
df = arcpy.mapping.ListDataFrames(mxd,"Layers")[0]
#set the dataframe
fc = ("N:\Lab13\cities.shp")
cities = arcpy.mapping.ListLayers(mxd, "Cities")[0]
cities.showLabels = True
def cityname():
    cityName = arcpy.GetParameterAsText(0)

    cities_layer = arcpy.MakeFeatureLayer_management(fc, "cities_lyr")
    arcpy.SelectLayerByAttribute_management(cities, "NEW_SELECTION", "CITY_NAME = '{}' ".format(cityName))

    arcpy.mapping.ListDataFrames(mxd)[0].zoomToSelectedFeatures()
    arcpy.RefreshActiveView()
    arcpy.Delete_management("cities_lyr")
    if cityName != "CITY_NAME":
        raise SyntaxError("Check your spelling or whether it is a European city.")
    arcpy.AddMessage("You zoomed in to " + str(cityName))
try:
    cityname()
except SyntaxError as se:
    arcpy.AddError("There was an error. " + str(se))


texts: Batch importing all vector types of hundreds of DXF files

texts:
When importing dxfs, you get a window asking you to select vector layers to add. With all the dxfs I've imported so far, that has had two layers called entities, one of which is points and one of which is linestrings.
I select all, click OK, and then the file imports. This has to be done for each dxf.
Today I had to import approximately 200 dxfs (which I eventually combined using the data management tools.)
Having to click on select all and then on OK, to select the vector layers for each and every one of those dxfs was very tedious!
Is it possible to set QGIS to select all automatically as it imports dxfs, without my having to sit there for half an hour repeating those clicks?
The dxfs were created from NTF tiles, using the LL2DXF converter. Unfortunately, many of them will not open in Autocad due to an invalid block name so I can't combine them in CAD if there are lots. I do have a script that can strip out the blocks, but it is a little slow and has to be used one file at a time. That's OK for a small number, but today I had a very large area of map to convert and there are other towns to do over the coming weeks.


texts: Only Black & White text in Python editor in ArcGIS for Desktop?

texts:
When opening IDLE or Pyscripter in ArcGIS by right clicking on a toolbox and Edit... the specified editor opens with the code - BUT all text is only black and white not the context colors of blue, red, green...
I've tried it with both editors to no avail.
When opening either editor outside of Arc the context color is as usual.

Has anybody experience this? How can we change this?

Windows 7, 32bit, ArcGIS 10.3, Python 2.7.8, Pyscripter 2.6.0


texts: Google Earth Engine interpolate produces strange results

texts:
I am trying to interpolate to increase the resolution of an image, but the results are not what I expected. I can't figure out what I am doing wrong.

var ols1 = ee.ImageCollection("NOAA/DMSP-OLS/CALIBRATED_LIGHTS_V4")
 .filterDate('2009-01-01', '2017-01-01');

var ols2 = ee.ImageCollection("NOAA/DMSP-OLS/CALIBRATED_LIGHTS_V4")
 .filterDate('1997-01-01', '2002-01-01');

var nl2012 = ols1;
var nl1997 = ols2;

var median1 = nl2012.reduce(ee.Reducer.mean());
var median2 = nl1997.reduce(ee.Reducer.mean());

var median = median1.subtract(median2);

var clipped = median;
var clipped1 = median1;

var zones = clipped.gt(1).add(clipped.gt(2)).add(clipped.gt(3)).add(clipped.gt(4)).add(clipped.gt(5)).add(clipped.gt(10));

// Define the visualization parameters.
var vizParams = {
  min: 1,
  max: 3};

var zonesMasked = zones.updateMask(zones.gte(4));

var image_10m = zones.resample('bilinear').reproject({
  crs: zones.projection().crs(),
  scale: 50
});

Map.addLayer(image_10m, vizParams, 'image_10m');

// Display the thresholds.
Map.setCenter(-74.6225, 4.712, 6);
Map.addLayer(zonesMasked);
Map.addLayer(zones);


texts:
I'm using OpenLayers - 2.11,GeoServer 2.13, ExtJS 3.4, and Java for developing GIS based web applications.

Recently i'm getting below error message in my application

"A script on this page may be busy, or it may have stopped responding. You can stop the script now, or you can continue to see if the script will complete." Mostly it showing OpenLayers.js file.

I've done"about:config" and set "dom.max_script_run_time" to 20 and more and nothing seems to work.

Please guide me to solve this issue..


texts:
I have a dataset, in shapefile format, for agricultural land grades in England. It is divided into 7 categories in its attribute table. As the file is a large file and I only want a small section I want to select a rectangle out of it. This was done for all other layers with the select rectangle tool. However it seems to select the entire data when i try to select the small area.

Is there a way to cut the data down to this small area which is undefined in the attribute table???


texts: Create XS cut lines having the coordinates of the cross sections

texts:
I am working with HEC-GeoRAS in ArcMap 10.3.1. and I am trying to create the cross section cut lines. I have the coordinates (X,Y) of the cross sections of the river.

Is it possible to create automatically the cross sections in HEC-GeoRAS having the coordinates or is it necessary to create manually?


texts: How can I avoid full scan big PostGIS table in QGIS

texts:
I try to open big PostGIS table in QGIS, with 500+M rows. I zoom in to small area, and limit zooms, so the query response is expected to be reasonably small. Tested this also directly in DB. The table has indexes, so I can see that raw query is fast. The table has already unique index for two fields ("partition","partition_offset")

Problem: When I add the table as layer, then QGIS hangs. I see in DB a long running query: SELECT count(distinct ("partition","partition_offset"))=count(("partition","partition_offset")) FROM "public"."mytable_geo" and from query planner I see that this does full scan to table, which takes possibly hours in my case. It looks like that QGIS wants to be really sure that my "unique" constraint makes the values really unique. Is there a way to avoid it?

QGIS version: 3.6.1-Noosa

Update: GDAL has parameter checkPrimaryKeyUnicity='1'. If this could be turned to "0" value in QGIS then it should help, but have not figured out how/if this is possible. Cannot find UI option for this.


texts: Can Spatialite SQL statements be issued with a python script in QGIS?

texts:
In my previous scripting environments (ms-access and Manifold GIS), it was possible to create and execute SQL queries from withing VBA or VB Script procedures.  This is really useful for things such as reloading a combo box on a form based on an entry in a previous widget/control. It is also good for creating update queries on the fly.

Given that spatialite can be run from a command line, I wonder if it is possible to use Spatialite SQL statements from within Python scripts in QGIS.  With ms-access and Manifold, they are both SQL aware and procedures will wait until a recordset is returned.  I wonder if this is the case with Python scripts in QGIS?

Andrew M


texts: Polygon Nodes do not have separate coordinates

texts:
I'm attempting to measure the distance of a series of points from a smaller number of polygons. As a first step, I've changed the polygons into points along their border, via the

Vector> Geometry Tools> Extract Nodes

command. This gives me the following result, after I've manually deleted irrelevant borders.

It's now my intention to export these points' coordinates into R to do some rudimentary neighbor and distance calculations. But when I open the attribute table, all the points representing the borders of the separate regions share a single geographic coordinate (which appears to approximate the polygon's centroid, but that's a guess). It's like the following:

Is there a different way to convert a polygon's borders to points, and recover the points' coordinates themselves?

Thanks


texts: Postgis reorder mixed up linestring chunks / max_segment_length(linestring)?

texts:
I imported ~1000 paths from somewhere into my postgis database to a linestring field.

I had the problem that the paths were divided into chunks, and those chunks were mixed up in some cases.

Suppose a linestring that was divided at point number 50 and 70:

Chunk A: points 1-50
Chunk B: points 51-70
Chunk C: points 71-100

When I migrated that into my database, they got mixed, so the resulted linestring could end up like this:

Chunk A: points 1-50
Chunk C: points 71-100
Chunk B: points 51-70

So that produces a jump from 50 to 71 and another one from 100 to 51

I want to be able to reorder those chunks so I would like to construct a SQL query to detect which paths are mixed, then I could manually rearrange them.

It would be desirable to have a SQL update query to solve this problem, but I think that detection is easier (I presume there are ~5% or less of paths with errors)

EDIT: I think the script for detection could check if a path contains a pair of consecutive points too far away. Maybe a SQL that orders paths from the path that contains the longest segment, would be good.

How could I make a function to get the length of the max segment in a linestring?


texts: Creating an app using the webmap app qgis plugin

texts:
I recently attempted to make an app using opensource software.

I tried re-creating an app following the http://suite.opengeo.org/opengeo-docs/qgis/webappbuilder/tutorials/schools/index.html tutorial.

However after saving it,instead of getting a "do you want to open it in a web browser" option i get "application files have been correctly generated. Use boundlessWebSDK for building the final web app from them."

i have opengeosuites installed on the computer but the boundless isn't showing anywhere in qgis.

How do i get the app to run?


texts: Select By Atribute using python script on .Shp file getting error 000840 value is not raster layer

texts:
I have a script that calculating vertices angles on Polyline layer and the selecting features in that layer at a specified angle. at the end of the script i'm using select by attribute
and getting error 000840 the value is not raster layer, the value is not mosaic layer

SelectLayerByAttribute_management(my_shp_name, 'NEW_SELECTION', myquery)


texts: How can I create a new layer from a one-to-many Relation in QGIS?

texts:
I am using QGIS 2.14 and am working with two layers: a layer of US Census Blocks and a layer of information on internet providers within each block. The internet provider layer is a CSV flat file, and the census block layer is a shapefile from the US Census TIGER site. I'd like to be able to create a unified layer of internet providers with information on each so I can, say, click "Verizon" and see everywhere Verizon offers service, or select just those blocks where the minimum speeds are in excess of X.

I understand I can't do a simple Join because it only joins one row of the CSV to each census block. I have set up a Relation in the Project Properties, which reflects the data structure I want, but I can't figure out how to make it just display in the map as a new layer.

I also tried creating a new layer by joining the CSV to the census block layer, but QGIS complained it wasn't a valid shapefile.


texts: Postgresql Connection issue "unexpected EOF on client connection"

texts:
I'm using Tomcat, Apache (for Load Balencing),OpenLayers - 2.11,GeoServer 2.13, ExtJS 3.4, PostgreSQL 8.4,PostGIS 1.5 and Java for developing GIS based web applications.

I have 2 Server.

Server1 contains

Apache(for Load balancing. This will divert all GIS links to ServerIP1:8081 and ServerIP2:8080)
Tomcat(ServerIP1:8080)
GIS Application
Geoserver Files
GIS Database(PostgreSQL 8.4,PostGIS 1.5)

Server2 contains

Tomcat(ServerIP2:8081)
GIS Application
Geoserver Files

All Map link will directly go to Apache from there after the load balancing process that will divert into Server1 or Server2.

LoadBalancer code in Apache

<Proxy balancer://clusterABCD>
    BalancerMember http://ServerIP1:8080
    BalancerMember http://ServerIP2:8081

    Order allow,deny
    Allow from all
</Proxy>
ProxyPass / balancer://clusterABCD/

ServerIP1 and ServerIP2 contains same GIS applications and the same Geoserver files.
Both Server accessing datas from ServerIP1:5433 Postgres Database Files.

While running select count(client_addr),client_addr from pg_stat_activity group by client_addr; this query
i'm getting from ServerIP1 Postgres counts are increasing continuously(count upto 50 or 100 or 200 like) and after some time intervals(not a fixed time)
that will come down to 1 then 2 then starting to increase the count.

If i stop ServerIP1:8080 tomcat then there is no problem.
Then i have checked which query is taking this much of connections in PostgreSQL But, i'm getting results only IDLE.
And getting only

2013-10-24 18:23:43 IST ERROR:  relation "public.gt_pk_metadata" does not exist at character 15

2013-10-25 06:15:31 IST LOG:  unexpected EOF on client connection

these errors in postgresql logs.

Geoserver Store:

If anyone faced a problem like this please guide me to solve this problem.


texts: ArcGIS: how to clear selection in Attributes Table with a keyboard shortcut?  I'm getting tired of going to Options, Clear Selection;  I work with Attribute Tables frequently.


texts: Problem with adding new field to attribute table after performing join by attributes

texts:
I am studying ArcMap 10.2 and concentrates on different join operations. I am working on a Personal Geodatabase.

I have digitized some polygons from a raster representing houses and have created a feature class under a certain feature dataset for this purpose. I have further created a table (not within any feature dataset) contains number of people occupy each house. I've used the ID in the polygon layer as a primary key and field named house_number in the table i've formed as a foriegn key.

After joining those two tables - I've received the attached new table:

After that I've tried to add a field (double precision) to use the field calculator, but it caused data to shift from one field to another, messing the attribute table. While in editor mode data is organized as it should be, but after exiting the editor mode it is being shifted back. Photos are attached:


texts: Exporting Selected Feature to New Shapefile using ArcPy?

texts:
I currently have a search cursor iterating through a shapefile in ArcGIS 10.1 that selects a feature and runs a viewshed analysis on that feature (and only that feature).
What's the easiest way to also export that same feature to a shapefile with the same name in Python?

fieldFID = 'FID'
arcpy.CheckOutExtension("Spatial")

arcpy.MakeFeatureLayer_management (inPoints, "pts")

with arcpy.da.SearchCursor('pts',[fieldFID]) as cursor:
    for row in cursor:
        fid = str(row[0])
        print fid
        arcpy.SelectLayerByAttribute_management ("pts", "NEW_SELECTION", '"FID" = {}'.format(fid))
        outViewshed = Viewshed(inDEM,"pts",1,"CURVED_EARTH",0.15)
        outViewshed.save("C:/temp/output/viewsheds/"+fid)


texts: Spatial join attributes based on closest point in PostGIS

texts:
I have two point tables from which I want to extract a value from one table to the other, based on the nearest distance of those two points.

t1               t2
id | value       id
1  | 100         1
2  | 105         2
3  | 120         3
4  | 102         4
etc

So if t2 with id 1 is closest to t1 with id 2 it should get the value 105. The t2 table has more entries than t1 and all points should get a value.

Using the following query I am able to get the distance to the nearest point:

SELECT round(ST_Distance(
 (SELECT ST_AsText(t2.wkb_geometry) FROM t2 WHERE id = 1),
 (SELECT ST_AsText(t1.wkb_geometry) FROM t1
ORDER BY ST_AsText(t1.wkb_geometry) <->
    (select ST_AsText(t2.wkb_geometry) FROM t2 where id = 1) LIMIT 1 )
 ))
as Distance;

How would I add the corresponding t1 values bases on the closest point to t2 for all the t2 rows?


texts: Misalignment of landsat images in GRASS?

texts:
I am new to Landsat images processing, and I'm trying to create land cover change maps over an area covering around 15 landsat scenes.

I downloaded all scenes available, unfortunately some layers don't align perfectly over time, even having same path and row, and same coordinate system (UTM22N). At first, I just deleted the misaligned scenes, because I had enough images available to get a correct cloud free map for each year I wanted.

But in a few cases (p227r057, p227r058 for example, but from what I've seen until now, I'm facing the problem for probably half of the landsat locations), I cannot find more that 2 or 3 correctly aligned images over more than 20 layers. It is probably just a vertical and horizontal shift, but the problem is that this shift is not constant in distance (from a few tens of meters to a few kilometers) and in direction.

Is there an easy way to cope with that problem, or is it necessary to proceed to the orthorectification of all (a few hundreds scenes multiplied by the number of bands) images?

I'm working with GRASS GIS.


texts: How to put an SVG symbol in my rules-based symbology?

texts:
I have added a points layer in which I have to put different SVG icons, according to the column level.

So far so good, the separation is according to the level using a default symbol, but I want to change this default symbol to be an SVG symbol.

my_rules = ( ('INICIAL','regexp_match("NIVEL",\'INICIAL\')','#B18904'), ('PREESCOLAR','regexp_match("NIVEL",\'PREESCOLAR\')','#00FF40') )
symbol = QgsSymbolV2.defaultSymbol(vlayer.geometryType())
renderer = QgsRuleBasedRendererV2(symbol)
root_rule = renderer.rootRule()
for label,expression,color_name in my_rules:
    rule = root_rule.children()[0].clone()
    rule.setLabel(label)
    rule.setFilterExpression(expression)
    rule.symbol().setColor(QColor(color_name))
    root_rule.appendChild(rule)
root_rule.removeChildAt(0)
vlayer.setRendererV2(renderer)

I am trying to create a symbol QgsSvgMarkerSymbolLayerV2,
but I do not know how to use it within the rules-based symbology?


texts: I was wondering if anyone is aware of any tools that are available to measure building lean in ortho images from within ArcMap. I am aware of calculations and processes that can be used to measure and calculate an angle, but not a specific tool to assist within ArcMap.


texts: QGIS set default style for layers based on name

texts:
I know that it is possible to apply default fill/line/point and colour ramp styles for layers within a project, However...

I have a model which spits out raster layers with defined names (e.g. NAME_XXXX_wd).  I don't want to apply a colour ramp but instead an actual style file.  Is it possible to set up a project that would automatically apply a style file to an added raster layer based on the name and how might you do it?

I'd also be happy to select the layer in the TOC and hit a toolbar button but would prefer something completely automatic.


texts: Unable to select FLOAT precision when adding column

texts:
I'm trying to do the lessons in the training-manual at qgis.org.

I,m having some problems with  3.3 Ratio Classification http://www.qgis.org/en/docs/training_manual/vector_classification/classification.html

according to the training-manual its supposed to look like this

but mine look like this

I am unable to select precision and width of FLOAT.

I'm using QGIS 2.2 Valmiera on an iMac.

Please exude my bad english.


texts: Using .tif files to get elevation and slope in ArcMap?

texts:
I'm using ArcMap 10.3 on Windows 10 and I have the 3D analyst, ArcScan, Geo-statistical analyst, network analyst, publisher, schematics, spatial analyst and tracking analyst extensions.

I have a number of .tif files from ASTER and I need to get a Digital Elevation Model  (DEM) of Ireland from them. The end goal is to get the altitude and slope for a number of points across Ireland. Here is a screenshot of these images in ArcMap:

These files are in a different coordinate system to everything else that I have - I have read that it's best to project them to the correct coordinate system after I change them to raster. Is this the correct thing to do? Or should I project them first?

I've used 'Mosaic to new raster' to try and bring these images together into something I can use. I had a value of 1 in the number of bands and selected 16 bit unsigned. This did not work. Here is an image of the result:

Am I on the right track with mosaic to new raster and just made a mistake, or is this not the way to go about this?

I have tried to project these tifs to the Irish Grid coordinate system (Project Raster), like the rest of my data. However, this doesn't seem to work - what's going on here?

in terms of trying to get the slope, I have been using the spatial analyst tool Slope - every time I use it I get this error, why?

Summary: What are the steps I need to take to use these .tif files to obtain altitude and slope?


texts: Duplicate Layer Names shown on Arctoolbox

texts:
I have a python script that generates a polygon feature class file and it automatically displays the layer in the table of contents of ArcMap. So whenever I select the layer file in any Arctool as a input, it is showing me duplicate names of the layer. And also, whenever I remove the layer from ArcMap, the lock file of the shapefile doesn't get released unless I close ArcMap down. The code are as follows:

arcpy.CreateFeatureclass_management(outShapefile, shapefileName, "POLYGON", "",
                                "", "", spatial_reference)
arcpy.MakeFeatureLayer_management(outLayer, shapefileName1)
arcpy.ApplySymbologyFromLayer_management(shapefileName1, outLine_Layer)
arcpy.SaveToLayerFile_management(shapefileName1, outLayer2, "ABSOLUTE")

addLayer = arcpy.mapping.Layer(shapefileName1)
arcpy.mapping.AddLayer(df, addLayer, "TOP")

Do you have any idea why its doing this? The most frustrating part is that the locked file persist once the layer is removed from ArcMap.(20171109_2323.shp.SABBILGEDUMP.5028.13292.sr.lock)


texts: What units are used for coordinates in geographic and projected coordinate systems?

texts:
I'm trying to grasp the basics of GIS but my literature seems to leave out some of the basic concepts (or it's just very vague).

What units are usually used for coordinates in a geographic coordinate system (such as WGS84) and what units are used for coordinates in a projected coordinate system?

I believe the answer to the first part of the question is degrees, minutes, seconds or sometimes decimal degrees (depending on the level of precision).

What about projected coordinate systems?


texts: Connect ArcGIS 10.6 to Oracle 12c Enterprise Edition

texts:
I want to connect ArcGIS 10.6 to an Oracle 12c Enterprise Edition. However, I found different informations on the Esri websites on how to do so.

I followed this website (http://desktop.arcgis.com/de/arcmap/latest/manage-data/gdbs-in-oracle/connect-oracle.htm) and executed it step by step. All worked fine, as far as I can tell. However, I got an Error message telling me that:

  Failed to connect to specified server. Do you want to continue? Cannot
  connect to database because the database client software failed to
  load. Be sure the database client software is installed and configured
  correctly.

I found this ESRI website https://desktop.arcgis.com/de/arcmap/latest/get-started/introduction/whats-new-in-arcgis.htm where it states that "Oracle Instant client will not work.", directly contradicting the instructions on website 1.

So what is right, respectivley, how can I connect ArcGIS 10.6 to Oracle 12c?


texts:
I am running into an error when running pdal translate on Mac OS High Sierra:

PDAL: Argument references invalid/unused stage: 'filters.voxelgrid'.

This is the command I am using based on a PDAL tutorial on thinning point cloud data:

$ pdal translate -i /Users/aaron/temp/1958-09-54.laz -o /Users/aaron/temp/1958-09-54-thinned.laz --filters.voxelgrid.leaf_x=4.5 --filters.voxelgrid.leaf_y=4.5 --filters.voxelgrid.leaf_z=4.5

Essentially, this command thins the point cloud data with a voxel grid algorithm.

I see that the filters.voxelgrid driver is enabled when I call pdal --drivers. I'm using pdal 1.6.0 which was installed using brew install pdal. I am also running pcl 1.8.1_2, which was installed by default during the brew install of pdal.

What might be causing this error and how can I resolve it?


texts: Mapping/visualising data with Long Tail?

texts:
I work with a dataset that gathers stats from all across a city. Each datapoint is of interest and worth visualising to our users, but there are a few parts of the city that significantly deviate from the average into hostspot clusters.

For example many wards throughout the city would have a value of [0-5] for any given month but one particular ward might be [300-400+]. Additionally there's not much mid range.

I'm trying to think about clever ways of visualising this that demonstrate both the data and it's month-by-month changes for most wards throughout the city, but also clearly demonstrate how much of an outlier particular wards can be.

I've tried:

Heatmaps - Radius needs to be set to a narrow zone to prevent spill over small boundaries, but due to the massive deviation of our outliers I end up with very pale, widely distributed wards across the city and then one tight dark blotch in the middle. It looks like a pitri dish and isn't very effective. Increasing the radius improves readability but means that datapoints now spill over ward/borough lines.
Choropleth - Much more visually appealing but again in a single swatch of colours (yellow for our [0-5]'s and dark red for our [300+]) I think it's still difficult to demonstrate the magnitude of the skew. My mind draws a linear line between each colour step.

Really I need a good idea for geographically representing data on a Logarithmic scale when I have a massive deviation in my data. Additionally data for the mode is quite low significance (0-5) and so is subject to a high variability...

Really struggling to think about how best to plot this. Perhaps the answer not to try and instead break it into separate graphs demonstrating each aspect (skew, significance of mode and peaks?).

Has anyone tried mapping a dataset like this before?


texts: Calculating field after arcpy.Append_management()?

texts:
I have some AutoCAD layers (f.e. Cad1, Cad2) and I want to export them into feature class in a geodatabase (f.e. feature_class). Also I want to calculate field in "feature_class" for objects from CAD.

I have Feature class "GISLayer" with "Name" field in my geodatabase. I want to add to this feature objects from CAD. And populate field "Name" in feature class with identical value "CAD_line" for all new features,that have been added from CAD

I use:

arcpy.Append_management([Cad1,Cad2], feature_class, NO_TEST)

It's ok, but I have no idea, how to calculate the field for CAD objects using python.
Any idea?


texts: CARTO Builder: show categorical data by no. of casualties

texts:
I have a categorical Carto map that shows the different types of industrial accidents, but I also want to differentiate the points based on the number of deaths i.e. the higher number of deaths, the bigger the marker width.

This is the CSS style, but I am not sure why I keep encountering an error when I apply the style.

 #layer {
   marker-width: 6.5;
   marker-fill: ramp([type_of_accidents], (#311271, #1D6996, #129C63, #73AF48, #EDAD08, #E17C05, #C94034), category(7));
   marker-fill-opacity: 0.5;
   marker-allow-overlap: true;
   marker-line-width: 1;
   marker-line-color: #FFF;
   marker-line-opacity: 1;
 }
 #no_of_deaths [ value <= 5] {
 marker-width: 5.0;
 }
 #no_of_deaths [ value <= 20] {
 marker-width: 15;
 }
 #no_of_deaths [ value <= 40] {
 marker-width: 25;
 }


texts: Show QGIS layer definition (qlr) filename in my project not the source file it was created from

texts:
I have a QGIS project with multiple qlr (layer definition files) in it which all source the same original data but all have different styling, labeling and filters applied. But when I bring these qlr files into my project the filename I get in my 'layers window' is the original source file which they are all sourced from!
So I have 4 files all named the same. Can I bring in qlr files and preserve the names they have and not display the source file?


texts:
I'm updating on Ubuntu 14.04 , but I have these tracebacks for QGis :

Preparing to unpack .../qgis-provider-grass_2.8.4+20trusty_amd64.deb ...
dpkg-divert: error: `diversion of /usr/bin/qgis to /usr/bin/qgis.bin by qgis-provider-grass' clashes with `diversion of /usr/bin/qgis to /usr/bin/qgis.bin by qgis-plugin-grass'
dpkg: error processing archive /var/cache/apt/archives/qgis-provider-grass_2.8.4+20trusty_amd64.deb (--unpack):
subprocess new pre-installation script returned error exit status 2

and it finishes like this :

Errors were encountered while processing:
 /var/cache/apt/archives/qgis-provider-grass_2.8.4+20trusty_amd64.deb
E: Sub-process /usr/bin/dpkg returned an error code (1)


texts: Distance measurements across UTM zones: use geographic or planar approaches?

texts:
I have a survey grid spread across 3 UTM zones (36N, 36S, 37S). I want to find the nearest (or shortest) distances of the centroids of these grids to roads and various points between.

It seems like there are too many compromises when using any kind of planar projection (read: with regards to preserving distance between any number of points on the map). Should one just forget about using projections in this case and go for goedesic, or ellipsoidal (read: Geographic) techniques?

Is there to anyone's knowledge a planar technique that will preserve distance between any number of points on the map? It does not seem like I can use an equidistant projection with exception of the gnomonic projection. Is this correct?


texts:
I was trying to convert a roadnetwork table.Its geometry contains both MULTILINESTRING and LINESTRING.I converted all the MULTILINE string to Linestring using the query

update all_india_merge set the_geom=ST_GeometryN(the_geom,1) where astext(the_geom) <> 'LINESTRING';

After conversion Some data has been lost.How can I recover those geometries.
Please Help me.


texts: SQL request with a shapefile data with PYTHON

texts:
I have a shapefile with some attributes (Hylac_id, Lake_name, ...) and i want to catch the polygone with a specefic Lake_name. The database than I use is available on the website https://www.hydrosheds.org/page/hydrolakes
I try with this method but if you have a completly different method I open to change my coding strategy:

I open the shapefile with the module ogr of python
I give a sql request to the method ExecuteSQL to shortlist my shapefile

this is the code that scheaches for example the lake Baikal
from osgeo import ogr

file2 = "/home/surface4/bernus/Base_Lac/HydroLake/HydroLAKES_polys_v10_shp/HydroLAKES_polys_v10.shp"
driver = ogr.GetDriverByName ("ESRI Shapefile")
ogr_ds = driver.Open(file2)
sql = "SELECT * FROM HydroLAKES_polys_v10 WHERE Lake_name='Baikal'"
layer = ogr_ds.ExecuteSQL(sql)

And when I print print layer[0]["Lake_name"], I get the first line of the database which is the Caspian Sea.
I try other request to see as:

"SELECT * FROM HydroLAKES_polys_v10 WHERE Hylac_id>5"
This do not work, here i do not only have the first line but all the database

"SELECT * FROM HydroLAKES_polys_v10 WHERE Hylac_id<5"
This work well

I have: "SELECT * FROM HydroLAKES_polys_v10 ORDER BY Lake_name"
This work well

So I do not understand what I do in bad way, I try my request in online SQL database to check if my syntax is correct.


texts: What is the meter square value of 1 pixel at a scale 10 in Google Earth Engine?

texts:
I am looking at the NDVI of various mangrove forests in the UAE.
Using the .gt function I have binarily separated forests into not canopy versus canopy, based on an NDVI-value threshold.
Then using ee.Reducer.sum() I managed to get the number of pixels in each mangrove site that are classified as "canopy".
If I wanted to know the size of each of those pixels in order to find out the total canopy size of each site, what should I do? I've tried the ee.Image.pixelArea() function but it doesn't seem to work or maybe I don't fully understand it.
My main train of thought was that just knowing the scale which is set at 10, shouldn't that be enough to find out the square meter coverage of 1 pixel?


texts: I am trying to get to my Google Earth extent in ArcMap without have to manually pan/zoom around. I do not need Google Earth imagery, just need to view the extent in arcgis. What would be the best way to do this?


texts: How can I see which layers my style is assigned to in Geoserver?

texts:
How can I easily tell which layers in Geoserver (2.3.5 at the moment) are assigned a particular style?  For instance, I have a style called "buildings"... How can I tell which layers have a "Default Style" of 'buildings'?

The only way I've been able to find is to go in and click on "Remove Style", which will then pop up a box that says "if you delete this, these layers will be affected..."   But if I accidentally click OK instead of Cancel, bad things would happen... is there an easier/different way?


texts:
I'm using EPANET in QGIS (extension QWater).
I have different layers :

pipe network (line vector)
household (point vector)

Firstly, I'd like to connect household automatically and perpendicularly to pipe network (create another line vector).

Is it possible?


texts: comparing rasters, cells that have changed value

texts:
I am attempting to compare two rasters representing a species distribution model for now and the future.

The rasters have the exact same extent and cell sizes.

here is when is gets tricky, they are multi-band rasters, so to create an attribute table I must first convert them to a binary raster (maybe there is another way around this?). Doing this I lose a lot of the information in the rasters.

I want a quick way to compare how many cells have changed value (using python probably) as I have many sets of these to compare.


texts: How to obtain the width and height of an image using arcpy

texts:
This Stack Overflow question shows how to obtain the pixel dimensions of an image using Python.

Is it possible to find an image's height and width using arcpy and ArcGIS?

arcpy.Describe would seem like a logical place to start, but I can't see whether it returns the height and width properties. I looked at File Properties and Raster Dataset Properties - am I missing something obvious?


texts: I'm trying to split polygons with overlaping line features, using Split Polygons tool from Topology toolbar. The problem I encounter is that only one from new polygons have old attributes and others only null values. I'm using arcgis 10.1.


texts:
I have two vector shapefiles: One file is points representing manholes and the other is a polyline file representing pipes leading to and from the manholes.  I need to snap the lines within a specified distance to the closest manhole.

Is there an easy way to accomplish this without doing it manually?

Update

As per Alexandre, I have exported my shapefiles to spatialite format.  Not sure what to next.

Update #2
I created a new db file and imported shapefiles into the new db using spatialite-gui,  I adjusted my SQL query but as you can see I get an error "no such column: f.geom"


texts: Can't get GRASS 7 V.kernel.vector to work on QGIS 2.18.14

texts:
I have reinstalled GRASS multiple times, with all of the necessary listed programmes to install beforehand (e.g. GDAL), and it won't work.

The error comes up - GRASS was not found in '/Users/Shared/unix/grass-qgis-mtlion/grass-7.0.4' (GISBASE), provider and plugin will not work.

I have tried to redirect it to multiple file locations - e.g. the application in the applications folder, and it still will not work.


texts: TypeError with UpdateCursor on a Field with a Coded Text Domain

texts:
Trying to run an arcpy.da.UpdateCursor on a field that has a Coded Text Domain.
I keep getting a typeError.
Not sure if there's a certain syntax that I need to be using when trying to update the row or what, but Esri documentation doesn't have anything.
Here's what I've got:
edit = arcpy.da.Editor(workspace)
edit.startEditing(False, False)
edit.startOperation()

# Create update cursor
with arcpy.da.UpdateCursor(fc, fields) as cursor:
    # For each row, evaluate the DateGroomed value, and update Status
    for row in cursor:
        if (row[0] != None):
            if (row[0] < datetime.datetime.now() and (row[1] == "Groomed" or row[1] == "Groomed Classic")):
                row[1] = 'Groomed Previously'
                cursor.updateRow(row)

edit.stopOperation()
edit.stopEditing(True)

The error flags on line 11 here (row[1] = 'Groomed Previously') with:
"TypeError: cannot update the table"
The domain value and description are both "Groomed Previously" for the value I'm trying to fill in.


texts: Launching high-detailed Ordnance Survey maps in QGIS 3x

texts:
I would like to have the Ordnance Survey maps tiles in QGIS.

Most of the sources from the web, from where we can import the maps are not in high detail.

The scale, which I am interested in is 1:500 (zoom level 15).
One reasonable source, from where I can enter these maps is:

mapping.cityoflondon.gov.uk/geocortex/mapping/?viewer=compass&runworkflowbyid=Switch_layer_themes&LayerTheme=Show%20the%20Explore%20The%20City%20layers

I have got the option to export them as the GeoTIFF, but I am restricted to the viewing area only.

Is there another source, where I can get these maps possibly for free as XYZ tiles of alternatively the WMS server?


texts:
I have been using gdal from the command line to convert an asc file to a GeoJSON output. I can do this successfully:

gdal_polygonize.py input.asc -f "GeoJSON" output.json

Now I wish to use Python and follow this process for a range of files.

import gdal
import glob
for file in glob.glob("dir/*.asc"):
    new_name = file[:-4] + ".json"
    gdal.Polygonize(file, "-f", "GeoJSON", new_name)

Hpwever, for exactly the same file I get the following error TypeError: in method 'Polygonize', argument 1 of type 'GDALRasterBandShadow *'

Why does the command line version work and the python version not?


texts: Plotting a csv file containing latitude and longitude into DEM (geotiff file)

texts:
I produced an elevation map. I want to plot my csv file named 'Coordinates.csv' that contains latitude and longitude points to my DEM (geotiff file) and the elevation map I produce will turn into latitude (y) and longitude (x) plot. When I import my DEM it is flip. I used gdal and matplotlib for plotting. How can I resolve this?
Here's my sample python code

from osgeo import gdal, gdalconst
import numpy as np
import matplotlib
import matplotlib.pyplot as plt
import elevation
import affine
import cartopy.crs as ccrs

filename = "Sample-30m-DEM.tif"
gdal_data = gdal.Open(filename)
gdal_band = gdal_data.GetRasterBand(1)
nodatava1 = gdal_band.GetNoDataValue()

# convert to a numpy array
data_array = gdal_data.ReadAsArray().astype(np.float)
data_array

# replace missing values if necessary
if np.any(data_array == nodatava1):
    data_array[data_array == nodatava1] = np.nan

# Plot out data with Matplotlob's 'contour'
fig = plt.figure(figsize = (12, 8))
ax = fig.add_subplot(111)
plt.contour(data_array, cmap=plt.cm.terrain,
            levels = list(range(0, 2500, 100)))

ax.set(xlabel='Longitude (Degrees)', ylabel='Latitude (Degrees)', title='Sample Plot')
ax.xaxis.grid(color='gray', linestyle='dashed')
ax.yaxis.grid(color='gray', linestyle='dashed')

cbar = plt.colorbar()

plt.gca().set_aspect('equal', adjustable='box')

plt.show()


texts: MakeLine with aggregate based on line id

texts:
I want to create a line between points 104 and 74. This works :

select 1 as axe, make_line(p1.geometry, p2.geometry) as geometry
from MR200_amers_plan p1, MR200_amers_plan p2
where p1.matricule = 104 and p2.matricule = 74

Now I can't get the same result with the aggregate function. I added a column ligne holding the line number, and tried:

select p.ligne, make_line(p.geometry) as geometry from
(select ligne, geometry from MR200_amers_plan where matricule
in ('104', '74') order by gid DESC) as p group by p.ligne

I can't figure out the reason even according to http://postgis.net/docs/manual-2.1/ST_MakeLine.html or http://www.gaia-gis.it/gaia-sins/spatialite-sql-4.3.0.html#p5


texts:
Isn't there a save method for QgsVectorLayer created for example like this:

i_out = processing.run("native:buffer",
    {'INPUT':'C:/GIS/data/testdata/ak_riks.shp',
    'DISTANCE':100000,'SEGMENTS':5,'END_CAP_STYLE':0,
    'JOIN_STYLE':0,'MITER_LIMIT':2,'DISSOLVE':True,
    'OUTPUT':'TEMPORARY_OUTPUT'})['OUTPUT']

print(type(i_out))
<class 'qgis._core.QgsVectorLayer'>

The way I save it now is:

processing.run("gdal:convertformat",
    {'INPUT':i_out,
    'OPTIONS':'','OUTPUT':'C:/GIS/data/testdata/fdsfdsds2.shp'})

Does not seem to be a save method:

for m in dir(i_out):
   if m.startswith("save"):
      print(m)

saveDefaultMetadata
saveDefaultStyle
saveNamedMetadata
saveNamedStyle
saveSldStyle
saveStyleToDatabase


texts: Choosing map projections specifically for San Francisco Bay Area?

texts: I am asking because I am not familiar with specific map projections for the San Francisco Bay Area.


texts:
In ArcGIS I have a set of routes, some of whose direction I change using the flip edit so that they are in the direction I require. I am able to make route event layers. However, why is it that after creating a route event layer, these routes whose direction was changed are found to be with their original direction.

workflow process is as follows; add autocad drawing to arcMap,edit properties of polyline features by setting layer=tracks; assign road codes to the tracks and create routes using roadcode attribute; edit the direction of the routes created;export as a layer; make route event layer with the exported layer using linear referencing tool and an excel table with road code attribute.

I need to have route events with the direction after flipping. How do I do this?


texts: XMLHttpRequest cannot load http://localhost: /geoserver/wfs. Origin http://localhost is not allowed by Access-Control-Allow-Origin

texts:
Possible Duplicate:
  Setup own Openlayers Proxyhost
  Error WFS Request, PROXY problems

I have this error:

XMLHttpRequest cannot load http:// localhost: /geoserver/wfs
Origin http:// localhost is not allowed by Access-Control-Allow-Origin

in my google chrome console when i make a WFS request. I use apache tomcat 6 and geoserver which run on 8081 and 8080 port. I read that this error is related to proxy.

How must to resolve my problem? I run my application on Windows 7 machine.


texts: How to calculate area based on a certain value range in Google Earth Engine?

texts:
So I'm trying to calculate the area with certain color, which those colors are representing certain range of value.

Here's the map

var indon = ee.FeatureCollection('users/sdwijati/IDN_adm2')
var Sesar = ee.FeatureCollection('users/sdwijati/Sesar_Jawa')
var nama_prov = ['Jawa Barat', 'Jakarta Raya', 'Jawa Tengah',
                  'Jawa Timur', 'Yogyakarta', 'Banten']
var prov = indon.filter(ee.Filter.inList('NAME_1', nama_prov));
var s = ee.Image("OpenLandMap/SOL/SOL_SAND-WFRACTION_USDA-3A1A1A_M/v02")
var savg = s.reduce(ee.Reducer.mean());
var gtsavg = savg.updateMask(savg.gte(45))

var n = ee.Image("OpenLandMap/SOL/SOL_WATERCONTENT-33KPA_USDA-4B1C_M/v01")
var navg = n.reduce(ee.Reducer.mean())

var xx = gtsavg.add(navg)

var sesar = Sesar
var bufferBy = function(size) {
  return function(feature) {
    return feature.buffer(size);
  };
};

var intervals =
  '<RasterSymbolizer>' +
    '<ColorMap type="intervals" extended="false" >' +
      '<ColorMapEntry color="#0CE84A" quantity="45" label="40" />' +
      '<ColorMapEntry color="#FFF300" quantity="62.658" label="62.658" />' +
      '<ColorMapEntry color="#FF1500" quantity="80.316" label="80.316" />' +
      '<ColorMapEntry color="#000000" quantity="97.974" label="97.974" />' +
    '</ColorMap>' +
  '</RasterSymbolizer>';

var sarbuff = sesar.map(bufferBy(25000*4))
Map.addLayer(xx.sldStyle(intervals).clipToCollection(sarbuff).clip(prov));
Map.centerObject(prov, 7)

I'm trying to calculate the area of the color black (value 80.316-97.974), the code below doesn't work how it's supposed to, my guess is because instead of polygon for the geometry, it's a feature collection (which is kinda doubtful because the shapefile itself is basically a collection of polygons).

var sgttinggi = xx.updateMask(xx.lte(97.974))
var sgttinggi2 = sgttinggi.updateMask(sgttinggi.gt(80.316))
var stinggi = sgttinggi.clipToCollection(sarbuff).clip(prov)

var area_pxa = stinggi.multiply(ee.Image.pixelArea())
                    .reduceRegion(ee.Reducer.sum(),prov,250,null,null,false,1e13)
                    .get('mean')

 area_pxa = ee.Number(area_pxa).divide(100000)
  print ('Area Dengan Potensi Kerentanan Sangat Tinggi (km²)', area_pxa)

Yes, the code runs and it got a result, but the actual area of the island itself is about 138.000 km sq, and the result is way bigger than that. So, any idea on solving the problem? (Anyway shout out to Solved Scholar on Youtube for the code)


texts:
var map = L.map('map').setView([55.42, 53.50], 6);

var wmts = L.tileLayer('C:/Program Files(x86)/GeoServer 2.15.1/data_dir/gwc/web/{z}/{x}/{y}.png', {
    wmts: true,
});
wmts.addTo(map);

The data I cache with geowebcache is stored C:/Program Files(x86)/GeoServer 2.15.1/data_dir/gwc/web/{z}/{x}/{y}.png.
How can I view this cached data with Leaflet my file formats image or png.


texts:
I have written the following code filter out some power data in overpass turbo.
[out:json][timeout:1000];
//fetch area "Tokyo" to search in
{{gecodeArea:Tokyo}}->.searchArea;
//gather results
(
  // query part for: "power=tower"
  node["power"="tower"](area.searchArea);
  //query part for: "power=line"
  way["power"="line"](area.searchArea);
  // query part for"power=substation"
  node["power"="substation"](area.searchArea);
);

//print results
out body;
>;
out skel qt;

However I get the following 2 errors,
Error: line 8: parse error: Unknown type "->"

Error: line 8: static error: For the attribute "type" of the element "query" the only allowed values are "node", "way", "relation", "nwr", "nw", "wr", "nr", or "area".

Any ideas what I am doing wrong?
Note to create -> symbol I am using the minus sign and greater than sign from the keyboard.