texts: Limit allowed request bbox of GeoServer layer

texts:
I am cascading a nationwide WMS with GeoServer and creating a local GeoWebCache layer from it.
I only want to make the local region available when the cascaded wms is requested and thus only chache whats within my own defined bbox. Is that possible?
I tried to limit the bbox by self defining Lat/Lon Bounding Box for the layer, but that did not work. Will try to do the same but set the Native BBOX, but I am unsure if that would do any difference.


texts: Calculate area of ellipsoidal polygon with QGIS

texts:
There is online service of GeographicLib, by Karney, where one can calculate the area of ellipsoidal polygon without need of reprojection.

Can one calculate the area of geometry in QGIS so as ellipsoidal polygons (in order to get the most accurate area), and not with $area that gives planar area after a equal-area projection of a geometry?


texts: Are .shp files coordinates in a certain direction?

texts: I'm working on a project where I'm importing a .shp file into SQL Server.  I'm using the Catfood.dll and I'm wondering if each point in the polygon are sequential (going in the same direction)?


texts:
In QGIS 3.12.2. I created a virtual layer to use for the sort-order in a print layout as described here:
Sort attribute table by hidden column in Print Composer

That works, but I have a problem with the geometry. The virtual layer does not seem to have a geometry and I was not able to create one in the edit field of the virtual layer (I tried set automatically as well as defining a point-geometry and the CRS of the original layer it is created from). In the layers panel, the virtual layer shows with a point-rendering symbol, however does not show on the map. Also zooming to the layer has no effect.

I need the geometry to insert the page-number of the atlas sheet each point apears on, see: QGIS expression to access page-number of atlas for certain features

Here is what I tried, based on the first link above:


texts: RuntimeError: Underlying DBMS error ORA-00936: missing expression schema.tableName STATE_ID is not column in table?

texts:
I created a search cursor to retrieve some data from a feature class.

The broken example shown below is used successfully in other scripts I created. The only difference that comes to mind is I'm retrieving records from a database table, not a feature class, in the successful cases.

The broken code throws an error on the line, for row in inCursor:. The message is

RuntimeError: Underlying DBMS error [ORA-00936: missing expression
] [schema.tableName][STATE_ID = 718361]. `STATE_ID` is not a column in my table.

I prefer to explicitly send a list of fields to the cursor as there are 50+ columns in this table.

The working example shown below is a result of a post I found. I tried it to ensure my database connection was good. The code prints out row data as expected.

What is wrong with the broken code sample? Any ideas why it works with database tables but not feature classes? Which is the better way to write a cursor? I'm not opposed to going back to rework some code using the now-broken cursor.

This code works:

fields = [f.name for f in arcpy.ListFields(inFC)]   # get field list
features = list(cursor_to_dicts(arcpy.SearchCursor(inFC), fields))

for feature in features:
    print(feature)

This code DOES NOT work:

inputFields = [ "Field1", "Field2", "Field3", "Field4", "Field5", "Field6" ]

with arcpy.da.SearchCursor(inFC, inputFields) as inCursor:
    for row in inCursor:
     # do stuff

Update #1

The feature class is in an enterprise geodatabase (Oracle). We may be using versioning but I can't confirm at the moment.


texts: What should be the spatial resolution of the classification output if I am using all the bands of Sentinel-2 sensor?

texts: Sentinel-2 sensor have 10 m, 20 m and 60 m pixel resolution depending upon the wavelength band. If I use all the available bands in my classification, what should be the resolution of the output classified raster?


texts: I would like to implement an arc node topology in arcmap. Do you know how can i build node topology from line layer in ArcMap 9.3.1?


texts: Geotools access attribute table on Geoserver

texts:
I want to access the attribute table of a SimpleFeatureCollection I have created within a WPS process by clipping it out of a Vectorlayer stored on Geoserver.

How can this be done? Is using an AttributeDescriptor the right approach? It's quite difficult to google and I can't find answers on this one.

Edit:
I tried it like this (my SimpleFeatureCollection is called test):

SimpleFeatureType type = test.getSchema();
AttributeDescriptor attributeDescriptor = type.getDescriptor(0); // just fixed index 0 for testing purpose
attributeDescriptor.getName();

Obvously the last line would need to get printed out in order to look at the result ..

But then I realized that this is probably a very wrong approach to get the attribute table.


texts:
I generated a rectangular grid using the vector tools available and added a attribute in this grade (sum of lentghs of a line).
Now I want to convert this rectangular grid to a matrix (if possible in a text, because I would put it in my thesis, so I wouldn't need to type cell by cell), in which the values of the cells would be that attribute.
I couldn't find a tool in QGIS to do that.
I'm using QGIS 3.8.3 Zanzibar.
This is the grid that I have generated, I did set the style according to a attribute value.

This is the attribute table.
I want to extract the values of the last column in a matrix form, in which the value of each cell would be the value of the attribute in that cell.

As suggested by Gabriel De Luca, I exported the vector in a CSV format and opened it in LibreOffice Calc. Now I'm closer to what I need. I have the column with 2500 rowns, but need to convert it to a 50x50 matrix.


texts: How to make a line shapefile converted from a gpx track file display over basemaps such as OpenStreetMap?

texts:
I have created a line shapefile from a Garmin gpx track file using QGIS 1.12.3 Lyon.  Resultant vector file displays properly with other vector shapefiles.  Problem is that created vector file (line shape) will not display over basemaps such as 'OpenStreetMap'.

After additional investigation, I have determined that created line shapefile has the correct projection but it is always overlain by any basemap.  For instance, if I set transparency of USATopo map to 50%, I can see the track lines displaying correctly through the map.  If transparency  is set to 0%, the track lines are blocked, they no longer display.  I have moved track line shapefile to the top of the Layer List in the Map Legend and still the tracks are overlain by any basemap.  What's really interesting is that labels for the track lines such as name and length display when the track lines themselves do not.


texts: How do I make equidistant circles covering the whole globe in QGIS?

texts:
What's the simplest way to make equidistant circles in PostGIS/QGIS? (my software versions are Postgres 9.3, PostGIS is 2.1, QGIS 2.10)

I tried this in PostGIS, but the circles that come near or across North pole get broken:

WITH series AS (SELECT generate_series(1, 10, 1) AS i)
SELECT
    st_buffer('point(82.9216 55.0292)'::geography, i * 1000000) buff,
    i AS id
FROM series;

Another attempt was to use MultiRingBuffer plugin, but it makes circles only in the selected crs, and this makes ellipses: 1,000 km wide east-west and 2,000 km wide north-south.

The horizontal segment is about 2,700 km, the vertical is 4,000 km, but they appear near the same ring. This ellipse has north-south "radius" (in real ground distance) twice bigger than east-west, in other words, it's not just a projection issue.

[edit] SOLUTION (see details in the answer below)

create a custom projection, here's my one

+proj=aeqd  +R=6371000 +lat_0=<city latitude> +lon_0=<city longitude>

create a ESRI Shapefile layer in this custom projection, with your city in the center (you might need another layer in usual projections with your city there, to know where to put the point).
create buffers with Multi Ring Buffer or something else (check for units, I got metres)
change project's projection to Pseudo Mercator (3857)
transform polygons into lines (vector tools) and cut out the part near 180 longitude, so that they don't cross the entire map


texts: Making a Feature Form for a standalone PyQGIS application as in QGIS

texts:
In QGIS, there is a method named openFeatureForm defined in QgisInterface class. It opens a form that shows value of attributes for a given feature.
layer = iface.activeLayer()
feature = layer.getFeature(0)

iface.openFeatureForm(layer, feature)

This method doesn't work in a standalone PyQGIS application. Because iface is predefined variable that refers QGIS interface.
How can I make a form like that showing the attribute values of a given feature? I also would prefer to filter the attributes.


texts: Why does QGIS spontaneously corrupt shapefile data?

texts:
I'm using QGIS 2.14 (Mac), but I have had the same problem with 2.8 (Ubuntu). When I am editing a shapefile, QGIS seems to spontaneously corrupt the data in unpredictable and catastrophic ways that include:

Duplication of features
Deletion of features
Scrambling of attribute data

Most recently, this occurred while using the Table Manager plugin to rename a column of the attribute table. Upon saving the change, one of the columns of the attribute table (not the one being renamed) was suddenly scrambled such that features were unpredictably reclassified, completely corrupting the file and forcing me to reload a backup version.

The problem seems to be associated with saving, though I usually don't discover the corruption until some time after it occurred, so I can't say for sure what the circumstances were in most cases, the exception being the case described above.

I realize this might be the same problem described here

QGIS corrupting shapefiles?

If so, it is evidently still a problem in 2.14, and there was never a satisfactory answer provided. If not, does anyone have a clue as to what is going on? This makes QGIS completely unusable.

-------------- EDIT ---------------

I was able to replicate the result by reloading a backup version and repeating the renaming and saving that precipitated the corruption the first time. The problem seems to be specifically with the Table Manager plugin, similar to what was reported here: https://hub.qgis.org/issues/13316

Apparently it remains unresolved, so I guess one should stay away from Table Manager for the time being.


texts: Create or move points at specific distances perpendicular to a line?

texts: I'd like to create new points or move existing points at specific individual distances (according to a field in the attribute table), perpendicular to a line. The points are evenly spaced out at 100 m.


texts: Calculating perpendicular line in ArcPy only works for shapefiles?

texts:
I'm trying to use this field calculator expression on field Shape to generate a perpendicular line, of 0.01m, at the midpoint of each existing line feature in ArcGIS:

CALCULATE FIELD (data management):
Input table: Clip_point_buffer
Field name: Shape
Expression: RotateExtend( !Shape!, 0.01)
Expression type: Python
Code Block:

Pre logic script code:
def RotateExtend(plyP,sLength):
l=plyP.length
ptX=plyP.positionAlongLine (l/2).firstPoint
ptX0=plyP.firstPoint
ptX1=plyP.lastPoint
dX=float(ptX1.X)-float(ptX0.X)
dY=float(ptX1.Y)-float(ptX0.Y)
lenV=math.sqrt(dX*dX+dY*dY)
sX=-dY*sLength/lenV;sY=dX*sLength/lenV
leftP=arcpy.Point(ptX.X+sX,ptX.Y+sY)
rightP=arcpy.Point(ptX.X-sX, ptX.Y-sY)
array = arcpy.Array([leftP,rightP])
section=arcpy.Polyline(array)
return section

Shape:
RotateExtend( !Shape!, 0.01)

(the code block i'm using comes from: https://gis.stackexchange.com/a/201871/63418)

It works if I do it with a shapefile, but if I have the same file as feature class in a database (.gbd) it doesn't. The database has a tolerance 0.001 and the default resolution (0.0001).
Why the process does not work when the feature is in a database and how can i do this process within a database.
This is the error I get:

I've read that it might have something to do with the properties of "Shape" type of field, because it's locked when in a database. If that's the error, how could I obtain the same results working with a feature class in a database?


texts: How to clip a raster using rasterio without lose resolution

texts:
I'm using rasterio to clip a large Tiff image using a shapefile, generating one raster for each feature.
However, the output resolution is low(96 dpi) and I would like to obtain a better resolution.

How should I configure it? Is it possible?Bellow is my function.

def mascada(pathIn, areas, pathOut):
    '''
    função que aplica a mascada na interpolação
    :param path:
    :param areas: gdf
    :return:
    :param pathin : caminho do raster
    :outmask: src.meta = traz o metadado da img
    :np.where: diz que quando nao tem dado considera 0
    :ff : file format
    '''

with rio.open(pathIn) as src:
    extent_geojson = mapping(areas['geometry'][0])
    outMask, crop_affine = mask(src,
                                shapes=[extent_geojson],
                                crop=True)
    soap_lidar_meta = src.meta.copy()
    return outMask

with rio.open(raster_entrada) as src:
        imagem = src.read
        out_meta = src.profile

Clipping process

raster_entrada = r"C:/Users/Ana/OneDrive - idgeo.com.br/Icrop/RESULTADOS/Mosaico_Geral_Entrega_01/Soub_Grupo_02.tif"
raster_saida = r"C:/Users/Ana/OneDrive - idgeo.com.br/Icrop/RESULTADOS/Clipadas/"

contar = 0
for index,poligono in gdf2.iterrows():
    pivo = gpd.GeoDataFrame(gpd.GeoSeries(poligono["geometry"]),columns = ["geometry"])
    pivo["id_pivo"] = poligono["id_pivo"]
    #pivo["fid"] = poligono["fid"]

#     print(poligono.Cultura)
    print(pivo.head())
    result = mascada(raster_entrada, pivo, raster_saida)
    meta = ({"driver": "GTiff",
         "height": result.shape[1],
         "width": result.shape[2],
         "count": result.shape[0],
         "dtype": 'uint8',
         "crs":crs,
         "transform" : transform
            }
            )
    with rio.open(raster_saida + str(pivo["id_pivo"][0]) + ".tiff", 'w', **meta ) as ff:
        ff.write(result)

Thank you very much for your time.


texts: WMS 1.3.0 invalid-format test failed even with proper/required code returned

texts:
I got the following test failed, though it looks that WMS 1.3.0 server responded with proper error code. See below:

Test Name:invalid-formatTest
Description:When a GetMap request contains a Format that is not advertised, then the server throws a Service Exception (code=InvalidFormat).Test Result:FAILEDTest Requests 1:
Method: GET
URL:    http://<host>:7003/mapviewer/wms?FoRmAt=UnknownFormat&VeRsIoN=1.3.0&LaYeRs=WMS_BASICPOLYGONS&CrS=CRS%3A84&BbOx=-0.0025,-0.0025,0.0025,0.0025&HeIgHt=100&ReQuEsT=GetMap&WiDtH=100&StYlEs=

Outputs:    ?

    <cnt:rest><response xmlns:xi="http://www.w3.org/2001/XInclude">
         <parser local-name="HTTPParser" namespace-uri="http://www.occamlab.com/te/parsers" prefix="parsers"/>
         <content>
            <response>
               <status code="200" protocol="HTTP/1.1">OK</status>
               <headers>
                  Thu, 11 Mar 2021 15:29:30 GMT
                  chunked
                  text/xml;charset=UTF-8
                  SAMEORIGIN
                  7040519b275824a90951f6bad613ae83b836a7c333932294
                  1
                  X-CSRF-MVS=7040519b275824a90951f6bad613ae83b836a7c333932294; Path=/mapviewer; HttpOnly; SameSite=Strict
               </headers>
               <content>
                  <ServiceExceptionReport xmlns:ns2="http://www.w3.org/1999/xlink" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" xsi:schemaLocation="http://www.opengis.net/wms http://schemas.opengis.net/wms/1.3.0/exceptions_1_3_0.xsd">
                     <ServiceException code="InvalidFormat">[wms] invalid FORMAT parameter.</ServiceException>
                  </ServiceExceptionReport>
               </content>
               </response>
         </content>
      </response></cnt:rest>

Reason of Failure:

Do you see any hints that failed the test?


texts: Create polygon label on both sides of polyline overlaping polygon

texts:
The image below represents areas which are as a final product exported to PDF as individual PDF's cut on the black line. The labels are created in a usual way.

I need to generate labels and after converting labels to annotation but in a way that I have labels on both side of black lines.

Every polygon that has black line over it should have two (or more) labels on each side of black line. The image below has been photoshopped.

Is this possible to do the same thing in ArcMap 10.3 using Maplex engine?


texts:
From the Wikipedia article:

  In geometry, an affine transformation or affine map or an affinity (from the Latin, affinis, "connected with") between two vector spaces (strictly speaking, two affine spaces) consists of a linear transformation followed by a translation.

Can anyone give a real-world example of when and how this is used in GIS?


texts: ArcPy ArcGIS Basic Standard: Get symmetric difference with partially contained, uneven numbers of polygons

texts:
I already looked up topics for rewriting the erase and symmetric difference functions, but I don't get it running.

I got the following situation:

FC1 and FC2 are my polygon input feature classes. FC1 consists of two polygons, FC2 only of one. I want following output:

So I want two feature class outputs: OUT1 and OUT2, which looks pretty much like a symmetric difference. So I did arcpy.Union_analysis((FC1,FC2), union). After this I worked with arcpy.SelectLayerByLocation_management(union, "CONTAINS_CLEMENTINI", FC1, "", "NEW_SELECTION") (also tried "CONTAINS"), but this gives me an empty output feature class after executing arcpy.CopyFeatures_management(union, FC1_only). Normally I would proceed with arcpy.SelectLayerByLocation_management(union, "", "", "", "SWITCH_SELECTION") and then already get OUT2. Why does this not work? And is this a right way to rewrite the Erase function only available for Advanced ArcGIS?

EDIT:
Ok, instead of going with the arcpy.SelectLayerByLocation_management stuff, I used the arcpy.Select_analysis on the union output with the query = '"FID_fc1" > -1 AND "FID_fc2" > -1'. Now I got the purple part shown above.
What is the best way to get two separate feature classes, one with the leftovers of fc1 and one with the leftovers of fc2?

EDIT2:
This is the code I have so far... I thought of maybe replacing the search cursor with an update cursor and handling the feature class output while still in the Loop by deleting the rows I don't want.... or do I have to create two queries? One with FID_fc1 > -1 AND FID_fc2 = -1 and one with FID_fc1 = -1 AND FID_fc2 > -1 ?

 elif desc.shapeType == "Polygon":
        # Initiate variables
        poly_union = outputFolder + "poly_union.shp"
        select_poly = outputFolder + "select_poly.shp"
        fid1 = ""
        fid2 = ""
        # Perform union on both polygon inputs
        arcpy.Union_analysis([new_data,old_data], poly_union)
        # Get the names of the FID fields for building the query later on
        with arcpy.da.SearchCursor(poly_union, ["NAME", "NAME_1"]) as cursor:
            for row in cursor:
                # If the NAME field is empty, just do nothing
                if row[0] == "" or row[0] == " ":
                    pass
                else:
                    # Use the name for building the FID_ field name
                    fid1 = "FID_{0}".format(row[0][0:6])
                # If the NAME_1 field is empty, just do nothing
                if row[1] == "" or row[1] == " ":
                    pass
                else:
                    # Check for the same names for the polygons to get correct FID_ field names
                    if row[0][0:6] == row[1][0:6]:
                        fid2 = "FID_{0}_1".format(row[1][0:4])
                    else:
                        fid2 = "FID_{0}".format(row[1][0:6])
        # Build the query with the extracted FID_ field names
        query = '"{0}" > -1 AND "{1}" > -1'.format(fid1,fid2)
        # Perform selection of the intersection of the polygons
        arcpy.Select_analysis(poly_union, select_poly, query)
        # ToDo: Delete intersection to get symmetric difference and output 2 feature classes


texts: Script for populating field in ArcGIS 10 with consecutive numbers (where start numbers is decided by me)

texts:
Possible Duplicate:
  How to add automate values to attribute field?

I need a script (for Arc GIS 10) by which I can populate a field against some selected polygons features.
 For instance, i selected several polygons serially caring the plot number serial(for the polygon features) then i separated the selected polygons in the attribute tables.After all i want to apply  a field calculator to populate the field by consecutive numbers where the start number is as i wish, say 1283128.


texts: Is it possible to render several DEMs simulatneously with QGIS2threejs?

texts: I'd like to know if anyone knows a way of including several DEMs in one representation using QGIS2threejs. If I have two separate DEMs, each representing one modeled surface, both from the same location (so, one is above the other), how can I get a visualisation where both surfaces are included, one "floating" above the other?


texts: Looking for high quality Koppen climate maps

texts: I wonder if anyone can point me in the direction of high resolution Koppen climate maps for use in academic research. I want to accurately place c. 100 study sites from around the world in their respective climatic zones for a meta analysis I am conducting.


texts: Problem with making a piechart in ArcMap , because of  values in attribute table?

texts:
I have a problem when trying to make a pie chart in ArcMap,

In the middle-right polygon there is no piechart of age groups. I think the reason is that in attribute table, there are some NULL values for this polygon (area)

The 80 and 90 age category in the piechart is probably the reason that there is no piechart in the  middle right polygon, and because there are NULL values instead of 0 in attribute table. Is there any way to change the NULL to 0, and to make the piecharts complete over the entire map?


texts: How to add text to LineString in OpenLayers?

texts:
We use OpenLayers 4.6.5.

The code looks as following:

const line = new Feature({
  geometry: new LineString(points)
})
line.setStyle(
  new Style({
    fill: new Fill({ color: 'black' }),
    stroke: new Stroke({ color: 'black', width: 2 }),
    text: new Text({
      text: '1',
      font: '18px "Roboto", Helvetica Neue, Helvetica, Arial, sans-serif',
      fill: new Fill({ color: 'black' }),
      stroke: new Stroke({ color: 'black', width: 2 })
    })
  })
)

this.vectorLayer.getSource().addFeature(line)

But, I'm getting the error:

Uncaught TypeError: textStyle.getFill is not a function
at _ol_render_canvas_TextReplay_.setTextStyle (textreplay.js?3128:470)
at _ol_renderer_vector_.renderLineStringGeometry_ (vector.js?98c2:183)
at Object._ol_renderer_vector_.renderFeature_ (vector.js?98c2:121)
at Object._ol_renderer_vector_.renderFeature (vector.js?98c2:94)
at _ol_renderer_canvas_VectorLayer_.renderFeature (vectorlayer.js?c03d:429)
at _ol_renderer_canvas_VectorLayer_.eval (vectorlayer.js?c03d:378)
at _ol_renderer_canvas_VectorLayer_.prepareFrame (vectorlayer.js?c03d:395)
at _ol_renderer_canvas_Map_.renderFrame (map.js?1103:183)
at _ol_Map_._ol_PluggableMap_.renderFrame_ (pluggablemap.js?0fa0:1180)
at _ol_Map_.eval (pluggablemap.js?0fa0:87)

Without text settings, there is no error, just a line.


texts: Calculating Net Primary Production (NPP)?

texts:
When I calculate Net primary production (NPP) for a local scale, what kind of software should I use?

I always use ENVI, ERDAS IMAGINE, ArcGIS for Desktop and IDRISI.

Is it possible to use these programs for NPP calculation?


texts:
The geoJSON-String:

{
  "type": "FeatureCollection",
  "features": [
    {
      "type": "Feature",
      "properties": {"title": "Standort 1 "},
      "geometry": {
        "type": "Point",
        "coordinates": [
          13.444175720214844,
          52.580522509085554
        ]
      }
    }
    ,...
  ]
}

This is my code to read and create the features to markers in openlayers:

var geo = new OpenLayers.Format.GeoJSON();

var obj = JSON.parse(featureJson);
var features = geo.read(obj, "FeatureCollection");

alert(features[0].attributes.length);

This works fine and I'm able to create the markers.

But what am I missing on this line?

for(var i=0;i < features.length;i++) {
    alert(features[i].attributes.title);
}

The output = "undefined".


texts: Where can i find oceanic borders defined as a set of lat/long points?

texts:
I'm trying to write some code that, given a latitude and longitude, will report which ocean (if any) this point is in.

Is there a public resource available that defines the borders of the oceans as a set of lat/long points that I can download?


texts: calculate turning acceleration from a compass

texts:
I have degree and timestamp and want to calculate the turning acceleration. What to trigger if the compass is turning fast.
Is there a mathematic way to get an acceleration parameter?


texts: Calculating polygon area within polygons in QGIS?

texts:
I am new to QGIS - I want to estimate the area of the reef under a polygon grid.
My shapefiles are already in meters, however I am having a hard time calculating the area under the grid layer. Attached are the photos.
I want to calculate the total area of the bold black polygons (reefs) under each grid square (image 2). How much reef is inside each square, for example, grid1: 25 sq.km, grid 1 sq.km, etc. Each reef polygon has its corresponding area in the layer attribute table.

ADDITIONAL PHOTO: RED IS THE REEF

THE ATTRIBUTES OF THE REEF LAYER


texts: Conditional statements in Field Calculator QGIS

texts:
Is it possible to modify the following syntax somehow that it would be applicable in Field Calculator in QGIS? The input layer is point layer and I would like to do following calculations which I can't do in Excel due to large amount of data.

if(  "VALUE1" < 0.01, 2,
  if( AND( "VALUE1" >= 0.01, "VALUE1" < 0.1 ), 3,
    if( OR( "VALUE1" >= 0.1, "VALUE2"  >= 2 ),  4, 0 )
  )
)

Any ideas?


texts: Color gradient in symbology in ArcGIS 10.1 for interpolation raster

texts:
I have measured noise and want to create with this measured points a noise map with interpolation....is there any way in symbology to create a flowing transition between the different colors for the created interpolation raster?

I used color ramp, so the colors have a gradiation but the borders between these colors are not fluent

or in simple words: i want to have a result like a rainbow, where the colors are fluent but without recognizable borders

thanks in advance!


texts: Encrypt geographic data in WFS using GeoServer

texts:
I have a layer containing polygon data which is supposed to be used in Google Earth. We want the users to see the polygons without getting their coordinates. Right now I solved this using a WMS service, so the user just gets pictures of the data. But all the scaling in Google Earth makes those WMS-layers, containing a lot of small and edgy polygons look pretty bad. If I load the data as polygons inside a KML the visualization is much, much better. Also the WMS-layer is not clickable, I would like the user to be able to click on the polygons and receive additional metadata.

If I use the KML-output of the GeoServer both is possible, but the KML contains the coordinates in plain-text. I know that it is not possible to encrypt kml-files, I searched around a lot, but maybe GeoServer offers a different solution?

I know that of course anything I transfer to the client can be copied somehow (he could also just copy the wms-images and recreate the polygons), but I need it at least safe enough that only "bad ass programmers" can get the data :)

Is there any way to hand out the data as polygons to the user, without him being able to just copy the coordinates of the polygons?

Or maybe a way to output the WMS in Google Earth with less "destroying" of the image?


texts: Copying feature classes from Personal Geodatabase to File Geodatabase using ArcPy?

texts:
I have huge amount of personal geodatabases that contain dataset and feature classes (one dataset with four feature classes and one feature classes outside of the dataset in each mdb).

I want to create a file geodatabase for each of mdb and copy the dataset and feature classes in the file geodatabases.

For creating file geodatabse i use the below code but i don't know how to copy  the datasets and feature classes to the file geodatabase programmatically.

import arcpy,os
for file in os.listdir(r"X:\test\2"):
    if file.endswith(".mdb"):
        arcpy.CreateFileGDB_management(r"X:\test\2",file[:-4],"CURRENT")


texts: How can i create points per 10 and 100 meters for several lines with ArcPy?

texts:
I have 10 lines and i want to generate points per 10 and 100 meters for each one of them. In ArcGis there is the tool construct points(COGO) but is not scriptable. Is there any suggestion?
Also i don't want to be saved in memory but in Hard disc, also i don't have the points i want to create them, with a loop, for several lines (10).


texts:
I'm trying to load tiles from the below server to QGIS by GDAL TMS support (like here: https://www.3liz.com/blog/rldhont/index.php?post/2012/07/17/OpenStreetMap-Tiles-in-QGIS), but it doesn't work - any hints on what's wrong with my XML?

<GDAL_WMS>
    <Service name="Kompass.de">
        <ServerUrl>http://ec2.cdn.ecmaps.de/WmsGateway.ashx.jpg?Experience=kompass&MapStyle=KOMPASS%20Touristik&TileX=${x}&TileY=${y}&ZoomLevel=${z}</ServerUrl>
    </Service>
    <DataWindow>
    <UpperLeftX>-20037508.34</UpperLeftX>
    <UpperLeftY>20037508.34</UpperLeftY>
    <LowerRightX>20037508.34</LowerRightX>
    <LowerRightY>-20037508.34</LowerRightY>
    <TileLevel>18</TileLevel>
    <TileCountX>1</TileCountX>
    <TileCountY>1</TileCountY>
    <YOrigin>top</YOrigin>
    </DataWindow>
    <Projection>EPSG:3857</Projection>
    <ImageFormat>image/jpeg</ImageFormat>
    <BlockSizeX>256</BlockSizeX>
    <BlockSizeY>256</BlockSizeY>
    <BandsCount>3</BandsCount>
     <Cache />
</GDAL_WMS>


texts: What are the Pros and Cons of hstore data type in postgis?

texts:
PostgreSQL has a module named hstore that implement the use of key/value pairs (it's something like the OpenStreetMap model, but in a single table).

What are the Pros and Cons of hstore regarding GIS compatibilities, data adding and retrieving?


texts: Jump to specific folder in catalog pane using Arcpy

texts:
I am developing a Python add-in for ArcGIS Desktop.

As a part of this add-in there is a button that I would like to jump me to specific folder in the Catalog Window.

Is there any way to do this?


texts:
I am trying to use Mapbox to display my maps.
a) I need to bind a popup on each polygon within a Tileset. I have done it in Leaflet using the onEachFeature function. The data should come from the Tileset's vector fields. Please focus on the Tileset as the source of information.
b) I also need to label the polygons using a field within the Tileset.
I simple script will do me just fine. This is where I am right now

<script>
        L.mapbox.accessToken = 'Token';
var map = L.mapbox.map('map')
.setView([-1.3328338, 36.8831222], 17)
.addLayer(L.mapbox.styleLayer('mapbox://styles/map/Tileset'));

var featureLayer = L.mapbox.featureLayer({
type: 'vector',
features: [{
type: 'Feature',

    onEachFeature: function (feature, layer) {

        var content = '<a "href="PlotMainDetails.apx?"' + feature.properties.plotno_1 + '" </a>"';
        layer.bindPopup(content);

    },
}]
})
.addTo(map);


texts: How to Set the Vertical Coordinate system when creating 3D Shapefile?

texts:
I have data of wells in CSV format, which I'm converting to a 3D shapefile, using ogr2ogr.
The latitude and longitude are in WGS84, while the depth is from the WGS84 geoid. I have the OSGeo4W installation of ogr2ogr.

The process completes successfully, and Points with XYZ are created;

However I can't figure out how to set the Vertical Coordinate system when creating this shapefile.

My CSV file is like this:

id, name, Lat, Long, depth
1, A, 12, 35, -10
2, B, 13, 40, -15
3, C, 15, 45, -20

My vrt file is like this:

<OGRVRTDataSource>
    <OGRVRTLayer name="data">
        <SrcDataSource>data.csv</SrcDataSource>
        <GeometryType>wkbPoint25D</GeometryType>
        <LayerSRS>WGS84</LayerSRS>
        <GeometryField encoding="PointFromColumns" x="Long" y="Lat" z="depth"/>
    </OGRVRTLayer>
</OGRVRTDataSource>

I'm using the following command to create the shapefile:

ogr2ogr -f "ESRI Shapefile" input.shp data.vrt

How to Set the Vertical Coordinate system when creating 3D Shapefile with ogr2ogr?


texts: Where can I find accurate DEM data of Brighton UK with an attribute table to be used in ARC GIS, being used for a university GIS project and have struggled to find this data.


texts: Converting TotalStation native field data to coordinates

texts:
I'm trying to convert field data from an older Sokkia TotalStation machine, to coordinates.  I've looked at Copan but as far as I can tell it doesn't support the file format.

There are some outdated software that referred to it as "reducing field book" with reference to "Booker Intl".

The field data is in tab delimited in the following format:

Can anyone advise me what format this is, and how I can convert this into coordinates?

If there are no off-the-shelf tools out there that can do it I can write something if someone can point me to the calculations.

Ultimate aim is to load it into QGIS and generate contours.

[Update 23/12 - Below some further details to clarify, in response to answer from @jbgramm]

I only need to resolve these into relative XYZ coordinates - they don't need to be georeferenced.  I can identify all the setup points on the geographic map, but that's not really relevant.
The screenshot I pasted is only an extract of the file, there are setup records for all the stations, and lots more observations from each station.
This specific file is just an example - I need to find a tool or approach to resolve other files like this one
I believe the records are indicated as follows: { = setup, % = backsight, e = observation
I suspect the "setup" record fields are: [station_id] { [date] [time] [reference bearing] [instrument height]
Fairly confident the "backsight" fields are: [backsight_target_id] % [horizontal angle] [vertical/zenith angle] [not sure] [target height]
Fairly confident the "observation" fields are: [obs id] [note] e [horizontal angle] [vertical/zenith angle] [slope distance] [target height]
All distances and heights are in meters
I believe you are correct that all "angular data appears to be Degrees.MinutesSeconds (DDD.MMSS)"
You are spot on regarding the backsight and checkshots.  It might be worth adding that the survey is in the southern hemisphere and the reference bearing is due South.

Ideally I would prefer an off the shelf tool that can read this format as is and convert it to X,Y,Z points for me.  But lacking that I wonder if you can direct me at material that reference the calculations.

Does it come down to calculating the polar for all observations?


texts: Zonal Statistics As Table with more than one zone?

texts: I am trying to sum up the cells of a raster that fall within a polygon of multiple habitat types and ten management authority zones. Is there a way for Zonal Statistics As Table in ArcMap 10.5 to use both the habitat type (gridcode) and the management authority (CMA_NO) as the Zone Fields (see image below)? I could split up the polygon into ten separate polygons, one for each management authority, then rejoin all the tables at the end, but I have about 50 rasters I need to do this for, so I would like to find some way of minimizing the number of steps.


texts: Clipping tiled raster with polygon using PostGIS?

texts:
I have a table in postgis with a tiled raster and a table with polygons overlaying the raster. I want to get the number of pixels for each category (pixel value) that fall inside a specific polygon.

If I use the following query the result I am getting is wrong because, from what I understand, it hasn't unioned the tiles of the raster:

SELECT (pvc).value, SUM((pvc).count) As total
FROM (
  SELECT ST_ValueCount(ST_Clip(myraster.rast, myvector.geom)) As pvc
  FROM rastertable As myraster, vectortable As myvector
  WHERE ST_Intersects(myraster.rast, myvector.geom) AND myvector.gid = 6
) As f
GROUP BY (pvc).value
ORDER BY (pvc).value;

With the following query the result I am getting is the right one, but it takes longer than it should, since it unions all the tiles before clipping

WITH
  unionedraster AS (SELECT ST_Union(rast) As rast FROM rastertable)
SELECT (pvc).value, SUM((pvc).count) As total
FROM (
  SELECT ST_ValueCount(ST_Clip(myraster.rast, myvector.geom)) As pvc
  FROM unionedraster As myraster, vectortable As myvector
  WHERE ST_Intersects(myraster.rast, myvector.geom) AND myvector.gid = 6
) As f
GROUP BY (pvc).value
ORDER BY (pvc).value;

I thought the following would work, but the results are the same with the first query (like no union took part).

SELECT (pvc).value, SUM((pvc).count) As total
FROM (
  SELECT ST_ValueCount(ST_Union(ST_Clip(myraster.rast, myvector.geom))) As pvc
  FROM rastertable As myraster, vectortable As myvector
  WHERE ST_Intersects(myraster.rast, myvector.geom) AND myvector.gid = 6
) As f
GROUP BY (pvc).value
ORDER BY (pvc).value;

What is the correct way to do this?

Upon further investigation, it seems that the problem lies in the clipping of the raster tiles. I managed to export the result of the clipping as an image using this query:

COPY (
WITH clipraster As (
    SELECT ST_Clip(myraster.rast, myvector.geom) As pvc
    FROM rastertable As myraster, vectortable As myvector
    WHERE ST_Intersects(myraster.rast, myvector.geom) AND myvector.gid = 6
)
SELECT encode(ST_AsPNG(
    clipraster.pvc
), 'hex') AS png FROM clipraster
) TO 'clip.hex';

and the resulting image was the following

The blue line is the polygon to which the raster tiles should be clipped.

However, it is obvious that some of the tiles where not clipped.

Any idea to what is the problem?


texts: GDAL reclassify and polygonize Esri ASCII grid files

texts:
I realized that reclassification is the tool to go to reduce the polygon count Creating simplified (binned) polygons from raster file? and implemented it in Java using geotools. Unfortunately this is really frustrating due to JAI https://stackoverflow.com/questions/44341018/jai-cant-execute-in-native-spark-only-in-sbt-and-as-a-separate-scala-function

So I want to try it instead using GDAL. Here is a dummy file. How can I interlink gdal_polygonize.py and http://www.gdal.org/gdal_calc steps is it required to write output to file?

ncols 4
nrows 3
xllcorner 0
yllcorner 0
cellsize 1
nodata_value -999
-1 -2 -3 -4
-4 -21 -3 -3
-3 -2 -31 -1

There should be 2 classes [-50 -10[, [-10 0[


texts: Convert raster outputs to data.frame format

texts:
I've moved over to version 2 of the lidR package. In previous versions, the output of grid_metrics was a data.frame. In version 2.0, the outputs are stored in raster format. I have additional code, such as running a PCA that works with the metrics in a data.frame.

Is there a suggested method for converting the raster products back into a data.frame?


texts: Leaflet: Container coordinates vs Layer coordinates?

texts:
What is the difference between map container coordinates and layer coordinates?

From what I gathered, the map container coordinates are always relative to the map viewport, and while layer coordinates always start out with the same values as the container coords, they change once the map is panned.

APIs that refer to coordinates/point:

LatLng to Layer Point
LatLng to Container Point


texts:
I have two MXD files, both of them have one feature layer with unique value renderer and multiple values.

Symbols are completely different but all of them are approximately 500 x 500 and 72 DPI. when I open them using ArcMap their memory usage are completely different, one of them use 450 MB and the other one user 1 GB of RAM.

What exactly is happening  behind ArcMap visualization, and what parameter is affecting the memory usage?


texts: Replicating the Layers Panel in my app, with sub-layers

texts:
I've followed this post, using the QGIS Desktop Python Console, to create color ranges for my vector layer. QGIS displays a nice selectable tree in the Layers Panel.

Now I want to develop my own app, instead of playing in the console. How can I replicate this tree? If I only had top-level layers, I think it would be straightforward to get the list of layers from QgsMapLayerRegistry and populate my own widget. But the ranges add complexity. As far as I can tell, the sub-layers exist only in the renderer, so do I really need to dig through the renderer's items somehow?

I'm stuck with QGIS v2.18 and PyQt4 in case it matters.


texts: Incorrect measurement of intersection between polygons

texts:
I have a map divided in parts using polygons, and I want to measure the length of the boundary between each pair. I am using this script from a tutorial to find the actual neighbors, and then measure the shared boundary with this code:
ring_f = QgsGeometry.fromPolyline(f.geometry().asPolygon()[0])
ring_intersecting = QgsGeometry.fromPolyline(intersecting_f.geometry().asPolygon()[0])
ring_intersection = ring_f.intersection(ring_intersecting)

Where f is the active feature and intersecting_f is iterated in a list of
all features intersecting f.
However, the script is giving incomplete measurements in some features. For example, below are, whose shared border is near 367m as measured manually using Qgis, but the script returns 224.91. Wierdly enough, what the script returns coincides with measuring between the red crosses encircled in green by me (it would be nice to now what this crosses are, by the way).

Does anyone have any idea of what might be happening?


texts: Storing a Mapbox-hosted vector tileset in a variable (Mapbox-GL-JS)

texts:
I'm trying to make a map using Mapbox GL that is basically identical to this example: https://www.mapbox.com/mapbox-gl-js/example/filter-markers/.

The main difference that source I'm using is a vector tileset (that I've uploaded to mapbox.com) with a specific source-layer that needs to be split up.

The problem is the first line of this part:

markers.features.forEach(function(feature) {
    var symbol = feature.properties['marker-symbol'];
    var layerID = 'poi-' + symbol;

    if (!map.getLayer(layerID)) {
        map.addLayer({

Since my data isn’t stored in a variable what do I put in place of markers.features?

Do I have to create a variable that references my source, source-layer, and map ID and use that?  How would that be structured?

Apologies if this is a total noob question, I tend to over think these things.


texts:
I am following the PostGIS Cookbook's code to create an NDVI. I created the function but when using the example code but when I callback the function there seems to be an error... and I cannot figure out what is wrong.
This is the code I used:
SELECT ST_MapAlgebra(rast, ARRAY[1, 3, 4]::int[], -- only use the red, blue and near infrared bands
                     'raster.landsat8(double precision[], int[], text[])'::regprocedure,
                     -- signature for callback function '32BF',
                     -- output pixel type 'FIRST',
                     NULL, 0, 0, '1.', -- L '6.', -- C1 '7.5', -- C2 '2.5' -- G
                    ) AS rast
FROM raster.landsat8 m;

And this is what I got:
ERROR:  syntax error at or near ")"
LINE 6:      ) AS rast
             ^
SQL state: 42601
Character: 330


texts: How to move address search box in qgis2web plugin?

texts: I am using 2.22.0 version of qgis2web plugin and when the plugin provide the address search box that is placed on the right side of the screen. I am using Leaflet map. Now I want to edit the script and move my searchbox to the left side of screen. My question is which script and which part of code I need to change in order to solve my problem? I searched for it in js folder but with no success.


texts: Various parameters in submitting SNAP InSAR jobs

texts: I am a newbie to the InSAR technique, and I am working with this currently. When submitting jobs, I have to set the values of Azimuth Coherence Window Size, Range Coherence Window Size, Azimuth Multilook Factor and Range Multilook factor. Indeed, I am not really clear about the meaning and the effect of those parameters on the results. Additionally, how can I know what the good values for those are?


texts:
I struggle with creating a deep copy of a memory layer in QGIS. Whenever I modify the layer copy, the original gets modified, too. Can anybody point me to the mistake I am making in the following code:

def createMemoryLayer(self, layerName, sourceLayer):
    # create empty memory layer
    memoryLayer = QgsVectorLayer(
        QgsWkbTypes.geometryDisplayString(sourceLayer.geometryType()) +
        "?crs=" + sourceLayer.crs().authid() +
        "&index=yes",
        layerName,
        "memory"
    )
    memoryLayerDataProvider = memoryLayer.dataProvider()

    # copy the table structure
    memoryLayer.startEditing()
    memoryLayerDataProvider.addAttributes(
        sourceLayer.fields().toList()
    )
    memoryLayer.commitChanges()

    # copy the features
    #memoryLayerDataProvider.addFeatures(
    #    [QgsFeature(feature) for feature in sourceLayer.getFeatures()]
    #)

    # DEEP copy the features
    features= []
    for feature in sourceLayer.getFeatures():
        f = QgsFeature(feature)
        f.setGeometry(QgsGeometry(feature.geometry()))
        features.append(f)
    memoryLayerDataProvider.addFeatures(features)

    return memoryLayer

see the commented-out lines for what I had thought to fix with the lines thereafter.
According to the API docs, QgsFeature(QgsFeature) and QgsGeometry(QgsGeometry) both create deep copy of the original object.

Any clues?


texts: While True: loop hangs after a few seconds

texts:
I am brand new to PyQGIS and am at the very beginning of a project to create a Python script that continuously logs data to a text file using the location from the GPS toolbar to look up the value of an attribute from a vector layer that corresponds to the coordinates of the GPS.  Ultimately the script will log data from the underlying vector layer based on the GPS coordinates as a tractor drives across an agricultural field.  Because I am at the beginning of the project, most of the functionality has not been implemented.  So far I have just tried to create the infinite loop that continuously looks up a value from the attribute table on the vector layer and prints it to the console.  Right now the XY coordinates are static.  What I am working on now is trying to get the loop to keep running until I hit the ESC key.  The loop works for about 5 seconds until QGIS "goes donuts" (i.e., becomes unresponsive while the blue donut circles on the screen).  Also, the ESC functionality to end the loop does not work.  I am cobbling together Python code from different Stack Exchange conversations, so maybe some of what I am doing with the While True: loop is not appropriate for QGIS.

import time
import msvcrt
layer = QgsProject.instance().mapLayersByName('RS18_19-CCS Plot Map 14 plots with Treatments and Rates Joined')[0]
while True:
    MyPnt = QgsGeometry.fromPointXY(QgsPointXY(-77.926819,40.724967))
    feats = layer.getFeatures()
    for feat in feats:
        if MyPnt.intersects(feat.geometry()):
            rate = feat['StarterRat']
            print(rate)
    time.sleep(0.5)
    if msvcrt.kbhit() and ord(msvcrt.getch()) == 27:
        break


texts: Make geographic squares that follow parallels in SQL Server

texts:
I have a list of points (with latitude and longitude) in an SQL Server database, and for each of these points I have to generate a 50 x 50 degrees square (the point being the South West angle of the square).
For instance, if I have the point 0,0 I create the polygon with the following SQL query:
SELECT geography::STPolyFromText('POLYGON((0 0, 50 0, 50 50, 0 50, 0 0))', 4326)
I expected that any point with a latitude greater than 50 would not be part of the polygon generated above. But if I check if the point 25,52 is within the polygon, the query returns true
SELECT geography::STPolyFromText('POLYGON((0 0, 50 0, 50 50, 0 50, 0 0))', 4326).STIntersects(geography::STPointFromText('POINT(25 52)', 4326)) -> true
The reason is because when I display the polygon, the shape is rounded at the top:

In this case, how can I generate a geographic polygon where the upper edge is straight and follow the 50th parallel? Should I use another projection?
Also Geometry data type is not an option for this project, I need to stick to Geography.


texts: How are cost and reverse cost computed in pgrouting?

texts: I downloaded the nyc osm metro extract from Mapzen. I then used osm2pgrouting to import that data into postgis, pgrouting and now I have all these columns cost and reverse_cost, but I'm not sure how they were computed. I want to do a special form of routing by modifying the cost columns, but first I want to understand how these quantities were computed. Is it based on the length of the streets?


texts: GetRasterBand() method for gdal in Python

texts:
I was going through a tutorial book called Python Geospatial Development. On the chapter on using working with geospatial data in python there is an example of a script meant to handle and analysis raster data for the height values. I ran the following code:

import sys, struct
from osgeo import gdal
from osgeo import gdalconst

minLat = -48
maxLat = -33
minLong = 165
maxLong = 179

dataset = gdal.Open("l10g")
band = dataset.GetRasterBand(1)

t = dataset.GetGeoTransform()
success,tInverse = gdal.InvGeoTransform(t)
if not success:
    print("Failed!")
    sys.exit(1)

x1, y1= gdal.ApplyGeoTransform(tInverse, minLong, minLat)
x2, y2= gdal.ApplyGeoTransform(tInverse, maxLong, maxLat)

minX = int(min(x1, x2))
maxX = int(max(x1, x2))
minY = int(min(y1, y2))
maxY = int(max(y1, y2))

width = (maxX - minX) + 1
fmt = "<" + ("h"* width)

for y in range(minY, maxY+1):
    scanline = band.ReadRaster(minX, y, width, 1, width, 1,
                               gdalconst.GDT_Int16)
    values = struct.unpack(fmt, scanline)

    for values in values:
        try:
            histogram[value] += 1
        except KeyError:
            histogram[value] = 1
for height in sorted(histogram.keys()):
    print (height, hsitogram[height])

but I got the following error:

Traceback (most recent call last):
  File "D:\Python\Progies\Geospatial\l10g\histogram.py", line 12, in <module>
    band = dataset.GetRasterBand(1)
AttributeError: 'NoneType' object has no attribute 'GetRasterBand'


texts: pgRouting - How get all pubs that are less than x meters from a point

texts:
I'm building a Mobile App to discover all pubs that are less than 2 km away from my position!

First I obtained all the roads of a town using Open Street Maps "osm2po"
Next I obtained all coordinates of the pubs using Qgis (Desktop app)

Now, I need to find all pubs that are less than 2km away from my position.

How can I do this with pgrouting?


texts: Getting what seems to be incorrect data from GeoNames

texts:
I have a hiking website.  I am trying to query for parks big and small around the world where people could potentially go walking or hiking.

I was trying to use the GeoNames data and was using query arguments from this table:

http://www.geonames.org/export/codes.html

Here is an example query I made trying to get parks:

http://api.geonames.org/search?maxRows=10&username=demo&q=PRK

and it gives me some things which are I have no idea what.  Here is also another thing where I try to query for arctic spaces:

http://api.geonames.org/search?maxRows=10&username=demo&q=LAND

And that gave places in the Caribbean :)

So I must be doing something totally wrong.  My goal was to look up all parks in the world, and also cities (with state/country info) so that I can make that connection on my site for my users.

Is my goal possible to achieve?  What was I doing wrong in my queries?


texts: Keeping live update coming from WMS using PyQGIS?

texts:
I'm using a WMS that contain weather radars that is update every 10 min.
I would like to know how to keep the map updating every 10min.

I tried:

from PyQt4.QtCore import *
timer = QTimer()
timer.setInterval(600000)
QObject.connect(timer, SIGNAL("timeout()"), qgis.utils.iface.mapCanvas().refresh)
timer.start()

but my WMS doesnt update.


texts: Adding drawing effects to layer using query in QGIS

texts:
I am working with a layer of administrative areas. The administrative area in focus needs to have no label and a drop-shadow around it; surrounding administrative areas do need to have the label. E.G. Andorra would have a drop-shadow and only 'France' and 'Spain' labels would be visible, 'Andorra' being filtered out.
Currently, I am achieving this by duplicating the layer and applying to one layer just drop-shadow whilst hiding surrounding polygons and in the other layer I am hiding the polygon in focus and show the surrounding labels.
How, and where can I run a query like
if (name = 'a',
    show labels for (name != 'a') AND apply drawing effects,
   )

This will avoid duplicating the layer for every administrative area I have to apply my styling.
Do you have any suggestions?


texts:
I get elegant transparent hillshades via a combination gdal and convert. Compare to grey-based hillshade commonly used, such transparent hillshades are very cool because they can be placed between the map background and other upper layers (roads, buiding) to provide a 3D feels whatever the background's type and color.

How it works

The trick: Starting with a grey scale hillshade produced by gdal hillshade, the trick is to take the grey channel's values, invert each value, and flow this result into a new opacity channel. Black pixel [0,0,0] become [0,0,0,255] (opacity=255), grey pixel [120,120,120] get lower becoming  [120,120,120,135] (opacity=135, aka 255-120), white pixel [255,255,255] become transparent [255,255,255,0] (opacity=0, aka 255-255) and so on. Hills' shades are opaque and black, plains progressively become (white) transparent. Conceptually, the pixel's band equation is something such :

alpha_value([0])    = 255-0 = 255;   // black pixels get full opacity
alpha_value([Grey]) = 255-Grey;
// grey pixels get 255-GreyValue opacity.
alpha_value([255])  = 255-255 = 0  // white pixels get null opacity !

For a laid back video on this approach, explained by a Photoshop designer, see Adding Shaded Relief in Photoshop (16mins).

Question

Taking an ETOPO or SRTM derived grey-based hillshade (gdaldem hillshade crop_xl.tmp.tif shadedrelief.tmp.tif -s 111120 -z 5 -az 315 -alt 60 -compute_edges, file available to download here) as input...

...How to do the trick cited upper via gdal or an other non GIS-destructive way on such .tif files ?

Note that I wish to keep GIS properties (geolocalisation).


texts: Creating point in ArcPy to use in Select by Location query?

texts:
I have a section grid for a US state. I also have a quarter/quarter section grid which splits the section grid into 16 equal squares.

I have written a Python script that finds the XY coordinates of the centroid of each quarter/quarter section, so 16 of them.

I need to create a temporary point feature at that centroid and then run a Select by Location on that point to see which quarter/quarter section it falls within and then I will be adding an attribute to that quarter/quarter section the point falls in.

How can I create this temporary point feature to be used in a Select by Location? Is there a way to create a single point in memory?


texts:
For a suitability analysis I need to define the suitability of slope from DEM values (in %) to build new trails.

Knowledges says that flat and rather steep areas are not suitable due to water runoff. Unfortunately, it is not possible to define threshold like

<5% nogo, 5%<x>60% go, >60% nogo.

So I was thinking about a fuzzy set. I use QGIS and do not want to install GRASS GIS (due to security settings on my macOS). SAGA GIS fuzzify allows to define gradients between 0 (fully flat) and 1, see

http://www.saga-gis.org/saga_tool_doc/2.2.3/grid_calculus_11.html

and

Basically, this works for me. But: For further calculations my output should be in the range from 1 to 5, whereas

0 - 1 = NoGo
1 - 2 = not good
2 - 3 = take care
3 - 4 = good
4 - 5 = must go

So I wanted to recalculate the 0 to 1 gradient from SAGA GIS into my gradient.

0   =! 1
...
0.5 =! 3
...
1   =! 4

How can I do this? Is there a better fuzzy logic algorithm available than this one from SAGA GIS (apart from GRASS GIS)?

If I didn't find any solution, I would probably create a discrete assignment between the slope and my categories.


texts:
During the process of doing kernel density in R, using UK crime data as data set.

Below is the work I have done so far but encounter some issues:

library(ks)
library(raster)
data1<-subset(data,select=c(Longitude, Latitude))
data1<-data[(complete.cases(data1)),]
Hpi <- Hpi.diag(x=data1,binned=TRUE,bgridsize=c(10000,10000))

Error in chol.default(S) :
 the leading minor of order 1 is not positive definite
In addition: Warning message:
In var(x) : NAs introduced by coercion

fhat.diag <- kde(data1, H=Hpi, binned=TRUE)

Error in as.integer(x) :
cannot coerce type 'closure' to vector of type 'integer'

str(data)

'data.frame':   6773 obs. of  13 variables:
 $ Crime.ID             : Factor w/ 6656 levels "00019a191b438730dfbccc03fc6494758f84bc095f1e0959450b40fef2586d48",..: 1453 866 100 969 3904 6480 4193 4657 1533 4786 ...
 $ Month                : Factor w/ 1 level "2014-02": 1 1 1 1 1 1 1 1 1 1 ...
 $ Reported.by          : Factor w/ 2 levels "City of London Police",..: 2 2 2 2 2 2 2 2 2 2 ...
 $ Falls.within         : Factor w/ 2 levels "City of London Police",..: 2 2 2 2 2 2 2 2 2 2 ...
 $ Longitude            : num  0.141 0.148 0.139 0.139 0.129 ...
 $ Latitude             : num  51.6 51.6 51.6 51.6 51.6 ...
 $ Location             : Factor w/ 4800 levels "On or near A1",..: 1769 3287 4795 4795 3327 3429 33 3808 236 157 ...
 $ LSOA.code            : Factor w/ 3218 levels "E01000002","E01000006",..: 15 15 16 16 16 16 16 17 17 17 ...
 $ LSOA.name            : Factor w/ 3218 levels "Barking and Dagenham 001A",..: 1 1 2 2 2 2 2 3 3 3 ...
 $ Borough              : Factor w/ 37 levels "Barking and Dagenham",..: 1 1 1 1 1 1 1 1 1 1 ...
 $ Crime.type           : Factor w/ 1 level "Burglary": 1 1 1 1 1 1 1 1 1 1 ...
 $ Last.outcome.category: Factor w/ 16 levels "Awaiting court outcome",..: 16 6 3 6 16 16 16 16 16 16 ...
 $ X                    : logi  NA NA NA NA NA NA ...

 class(data)

 [1] "data.frame"


texts: CARTO: centered icons instead coloured countries

texts:
In CARTO I'm making a tourists per country map

Instead of colouring countries, I'd like to put a centered Bagpack icon in the middle of every country, of size depending on the amount of tourists.

When I use Centroid analysis, it deletes my amount_of_tourists column. So how to do that?


texts:
I'm using gdal lib for Python in a Jupyter Notebook environment.
I'd like to retrieve the resolution of a raster to use it in a gdal.warp call.
At the moment I do the following :
src = gdal.Open(mask)
ulx, xres, xskew, uly, yskew, yres  = src.GetGeoTransform()

but I get too much information. Is there a way to only retreive xres and yres?


texts: CRS of layers and google areal base map in QGIS 2.18- Error

texts:
I try to load layer of population in U.S.A on base map of google earth satellite. this this part of the attribute table:

My problem is that the point located in worng place:

I tried to change the CRS in order to locate the point in the right place- but with no success. Now the CRS is EPSG:3857

I red Can someone clear up a Google Maps/ArcMap lat/long issue? but didn't understand how can i solved my problem.

But when i load the layer (csv file) and change the layer CRS into EPSG 4326 the point located in the right place (i think...)

When i open the project properties there is WGS84 grid:

I really confused from this issue- how can i change the CRS so it will be the same as the project + layer+base map of google


texts:
I have about 10,000 points that I want to apply a partial circular buffer to based on an angle,  from -30 to 135°.

A similar question for polygons seems to have been answered (How to create an oriented buffer using arcpy?), but I cannot make sense of how I would go about implementing that for my scenario.

What steps would I take in ArcMap?

Converted points to raster, ran focal statistics with wedge from 210 to 315.


texts: Leaflet image overlays and iOS 12: partial rendering

texts:
My PWA (using Leaflet 1.5.1) is currently running (since 2 years) on Android and iOS without problems, until iOS 12 was released. From that point on I have very weird issues on iPhone and iPad.

Main problem is that image overlays are not fully rendering on iOS devices.

I am setting up an array of overlays:

ovl = loadOverlay(url);
overlay.push(ovl);
where loadOverlay() is:

function loadOverlay(imgurl) {
    bounds = new L.LatLngBounds(
        new L.LatLng(30, -10),
        new L.LatLng(50, 36)
    );

    var ovl = new L.ImageOverlay(imgurl,
        bounds, {
            pane: 'imagebg',
            attribution: "my attribution"
        });
    return ovl;
}

Then, each overlay is displayed this way:

startMap() {
 mymap = new L.Map('CopMap', {
 maxZoom: 6.50,
 minZoom: 4.75,
 zoomSnap: 0.25,
 zoomDelta: 0.25,
 zoomControl: false,
 maxBounds: [
        //south west
        [32.5, -4.],
        //north east
        [46.8, 31]
     ],
        crs: L.CRS.Simple
 });
    mymap.setView([41.5, 13], 5.25);
    mymap.attributionControl.setPrefix('');

    mymap.createPane('imagebg');
    mymap.getPane('imagebg').style.zIndex = 1;

    mymap.addLayer(overlay[iDX]);
}

This is working fine on Android, while in iOS the image is not fully rendered (not always, but very often). I have checked that the image file is correctly downloaded in the cache, but not fully displayed


texts: Mapping Superfund Sites using ArcGIS for Desktop?

texts: I am trying to create a map of all Superfund sites in the U.S., and of course it isn't as easy as I thought it would be. I think I have downloaded all the necessary information, but can't figure put how to create this map. I have a basemap, I downloaded the shapefiles from the EPA, but now I am stuck. When I unzip the shapefiles, I have a bunch of excel tables, some of which do contain latitude, longitude information. I just need to do some preliminary mapping for a future project.


texts: PYQGIS 3 Using the ouptut from GRASS v.distance algorithm

texts:
I have a script which is trying to run v.distance on two existing shape files,
buildings.shp which containing some polygons and singleasn.shp which contains a single point. This will be run from the console using QGIS 3.4

Using the algorithm via the Processing toolbox works fine, but the code has a problem.

I'm using the following python code:

e = buildings_layer.extent()
extents_str = str(e.xMinimum()) +"," + str(e.yMinimum()) +"," + str(e.xMaximum()) +"," + str( e.yMaximum())

parameters = {'from':'C:\QGIS\Projects\QGIS3 Dev\area 1\buildings.shp',
             'from_type':[0,1,3],
             'to':'C:\QGIS\Projects\QGIS3 Dev\area 1\singleasn.shp',
             'to_type':[0,1,3],
             'dmax':-1,
             'dmin':-1,
             'upload':[0],
             'column':['calculated'],
             'to_column':None,
             'from_output':'memory:nearest',
             'output':'memory:distance',
             'GRASS_REGION_PARAMETER':extents_str,
             'GRASS_SNAP_TOLERANCE_PARAMETER':-1,
             'GRASS_MIN_AREA_PARAMETER':0.0001,
             'GRASS_OUTPUT_TYPE_PARAMETER':0,
             'GRASS_VECTOR_DSCO':'',
             'GRASS_VECTOR_LCO':'',
             'GRASS_VECTOR_EXPORT_NOCAT':False}

algout = processing.run('grass7:v.distance', parameters)
for ft in algout['output'].getFeatures():
    #do something here#

The algorithm appears to run ok, but when it reaches the for loop i get the following error:

  AttributeError: 'str' object has no attribute 'getFeatures'

Similar code to this worked for a native algorithm but not for a grass one - is there a difference in the output and if so how do I address it?


texts: I'm looking at plant biomass values in 2m intervals (center to 30m) around several waypoints. I've extracted the zonal stats as a table, however each of the buffers also includes the smaller increment buffer values (i.e the 8-10m interval also includes the center-8m values). Is there any way to create multiple annulus around each point so that I'm isolating the values per interval and not the total? Using ArcMap/GIS 10.3.1.


texts: How do I aggregate multi-part features based on a field value?

texts:
I have a city park layer that has several multi-part features. I would like to dissolve these into single part features using the Aggregate Polygon tool. I cannot use the Dissolve tool because most of these multi-parts are not directly adjacent to each other. There are highways and waterways separating the parts by up to 300'. But the individual parts have the same ParkID value, so I would like to use this field to identify the aggregated shapefile records.

In the second image, you see green polygons that are not park of the red park. I would like to model or script a tool that would aggregate these too, based on their respective park names/IDs


texts: display lon and lat of current selected CRS in Map Tip

texts:
I would like to set Map Tips to display coordinates of some feature in custom selected CRS. When I do:

<b>X: </b>[% $x %]<br>
<b>Y: </b>[% $y %]

Tip displays coordinates in WGS84, but not in my current canvas CRS.


texts: Can I open remote datasets with Rasterio?

texts:
I know I can open remote datasets with Gdal, I was wondering if I can do the same with Rasterio?
Gdal bindings code of what I mean by opening remote dataset:
import gdal
import numpy
import matplotlib.pyplot as plt
from gdalconst import GA_ReadOnly
mmap_name = '/vsicurl/https://esgf.nccs.nasa.gov/thredds/fileServer/BioClim/WRE/bio11_equiv_2053.tiff'
dataset = gdal.Open(mmap_name)


texts: Creating frequency graph in layout view of ArcGIS for Desktop using Quantile symbology?

texts:
I'm interested in creating a graph showing the frequency of my data based on the categories of my Quantiles symbology. The Graphs Manager does not seem to allow me to do that, so I was wondering if there is another way to do so.

For example this is an image of the different classes I have using the Quantile symbology,

and I would like a histogram/vertical bar graph showing a count by class.

The current "Create Graph" wizard only allows me to create graphs based on my original dataset in the attribute table. My original attribute table does not show the breakdown by class.

Does anyone know a solution to this?


texts: Return column number of min value in dataframe

texts:
So as a result of a bunch of preceding calculations I end up with a dataframe of varying numbers of rows and columns.

I want to find the minimum value of the df and return the row and column index. How can I do this ? For example the below data is a typical output:

row.names   2838        2849        2865
2768        9118.4184   14548.0008  6118.441
2769        10749.0109  18873.3781  14490.520
2775        11336.9757  19230.1397  13515.399
2780        3535.5576   11619.9425  5497.814
2842        2946.4059   10614.0673  5994.932
2843        127.0874    5250.5626   6012.107
2844        11684.4736  3851.1980   13670.723
2846        449.1957    149.3105    6171.846
2862        8574.5068   8143.8069   14397.451

This line will correctly give me the row index of the lowest value (127.0874):

    > which(df==min(df))
    [1] 6

But what command will return me the column index for the minimum value?


texts: Join GPS data and video from UAV to get frames from specific locations

texts: Is there any piece of software capable to present gps track simultaneously with video ie. google maps data and allowing to clip a frame on selected locations on track? Both data and video are from Phantom 3 Advance


texts: ArcPy: Distortion towards edges when exporting to TIFF

texts:
I am writing a function to zoom to selection (based on a criteria) and Exporting it to TIFF.

The problem I am getting is that the when I load the exported TIFF image in Map, I am seeing distortion towards the edges (perfect near the center). Can someone help me with what am i doing wrong in this case. Attaching my code below:

def SelectDataZoomSelected(mpath, fpath, cond, opath):

# Import system modules
import arcpy

mxd = arcpy.mapping.MapDocument(mpath)
df = arcpy.mapping.ListDataFrames(mxd, "Layers")[0]
lyr = arcpy.mapping.ListLayers(mxd, fpath, df)[0]

# Use the SelectLayerByAttribute tool to select and zoom to the selection
arcpy.SelectLayerByAttribute_management(lyr, "NEW_SELECTION", cond)
df.zoomToSelectedFeatures()

# Export the map to a TIFF
arcpy.mapping.ExportToTIFF(mxd,opath,df,640,480,96,'TRUE','24-BIT_TRUE_COLOR','NONE','TRUE')

# Clear the selection and refresh the active view
arcpy.SelectLayerByAttribute_management(lyr, "CLEAR_SELECTION")
arcpy.RefreshActiveView()

del mxd


texts: How to write a file in an Server Object Extension (SOE)

texts:
I would like to create a file with an ArcGIS Server Object Extension and return the url.  Is this possible?  Where would I write the file?

I'm working with .Net, IIS, and AGS 10.1.


texts: Is it possible to open MapGIS files ( .WL, .WP & .WT ) in QGIS or other GIS?

texts:
We received some files from clients located in China created in a program called MapGIS. I was wondering if anybody knows how to make QGIS or any other GIS program, excluding MapGIS, open these files?

Our office currently does not own a license of MapGIS.

Each map comes with the following four file types: .MPJ, .WL, .WP & .WT.  Internet research has told the .WL is linework files, .WP is points files.

A link to the MapGIS program is here http://www.mapgis.com/index.php/index-show-tid-113.html


texts: QGIS 3x Viewshed plugin no result after re-projection

texts:
Continuing the query here:
QGIS 3x viewshed plugin returns wrong results for visibility map

I was trying to set proper CRS for my DEM on the area, where I am currently working on.

Unfortunately, the reprojection wasn't successful.

After changing from EPSG:4326 (WGS84) and EPSG:27700 (1936 British National Grid), the QGIS says:
Cannot use preferred transform between EPSG:27700 and EPSG:4326

In the details I can find, that this transform is not available.

Basically, the layer is created, although I am getting nothing.

Is there any way to fix it?

EDIT:

I found, that I am not only one with this problem, which shows the question below:

QGIS 3.10 projection error using BNG27700

Sticking to the edvices there I downloaded the relevant file from the OS website

https://www.ordnancesurvey.co.uk/business-government/tools-support/os-net/for-developers

and next, I have done the following steps:

Reprojection of the layer
Saving the SHP points again under the British grid EPSG 27700, making them visible for the plugin (otherwise I was getting the error stating, that "there are no points here", whereas there were projected for EPSG:4376
Generated the view with changes in the data attribute table (image below).

Unfortunately, the result is not acceptable, since I get the blank black canvas with 2 white smallish rectangles showing the location of my points.

What can I do to get my result right?


texts: Reclassifying field using Python wildcard string search?

texts:
I have a point layer with 900K Tree Point records.  In that layer I have field (d_Species) that contains the tree species of each record. I've created a field (SpecieValu) to assign a value (0-3) for each record based upon the species of the tree point.  There are hundreds of unique species and so I need to assign values based on a substring in the species field like "oak" or "cedar".

How do I use a wildcard search to find a substring and return a value?

My new field (SpecieValu) is a Double data type and the field I'm trying to reclass (d_Species) is a String data type. This is what I have so far but doesn't return anything but else 0:


texts: Programmatically determine if a local geodatabase is File or Personal

texts:
Possible Duplicate:
  How to distinguish between File and Personal Geodb workspace?

How can you determine whether a workspace is a file or personal geodatabase?  Using IWorkspace::Type you can determine wether the geodatabase is a local or remote, and using IDatabaseConnectionInfo2::ConnectionDBMS you can determined the DBMS for an SDE geodatabase.  However, I can find no class members that tell me directly whether a local goedatabase is File or Personal.  The only way I can see to do it is to retrieve the property set from IWorkspace.ConnectionProperties and evaluate the DATABASE property to see if the path string has a .gdb or .mdb extension.  My guess is that there is better way to do this, and I just haven't found the right interface.


texts: How can I split a polygon set with a line set using PostgreSQL recursive queries?

texts:
I have a simple polygon layer GEOMETRY(POLYGON,4326) and a simple line layer GEOMETRY(LINESTRING, 4326). I want to find all split areas. I am using PostGIS 2.0.

I am using following PL/pgSQL procedure:

CREATE TEMP TABLE cursor_tmp (geom GEOMETRY(POLYGON, 4326));

DO $$
DECLARE
  _geom GEOMETRY(LINESTRING, 4326);
  _curs CURSOR FOR SELECT geom FROM blades;

BEGIN
  OPEN _curs;
  LOOP
    FETCH
    _curs INTO _geom;
    EXIT WHEN NOT FOUND;

    TRUNCATE TABLE cursor_tmp;

    INSERT INTO cursor_tmp
      SELECT (ST_Dump(
          ST_CollectionExtract(
              ST_MakeValid(ST_Split(poly.geom, _geom))
              , 3))).geom :: GEOMETRY(POLYGON, 4326) geom
      FROM poly
      WHERE st_intersects(_geom, poly.geom);

    DELETE FROM poly WHERE st_intersects(_geom, geom);

    INSERT INTO poly SELECT * FROM cursor_tmp;

  END LOOP;
  CLOSE _curs;
END $$;

This gives me the desired results:

A colleague suggested that there might be a more elegant solution to this problem using PostgreSQL recursive queries. The problem is that I have zero knowledge of recursive queries.

Is it worth investigating? Will a recursive query be able to do this process?


texts:
I've added vectorlayer from database using the following code
# Create New Qgis Data source uri
uri = QgsDataSourceUri()
uri.setConnection(database_connection['host'],
                  database_connection['port'],
                  database_connection['database'],
                  self.db_user, self.db_password)

uri.setDataSource(schemaname, tablename, geometrycolumn, query)
vlayer = QgsVectorLayer(uri.uri(False), tablename, "postgres")
vlayer = self.iface.addVectorLayer(uri.uri(False), tablename, "postgres")

I'm unable to edit the attribute or geometry of the layer that is getting added. The below is the capability string for the layer
print(caps_string)
# Add Features, Add Attributes, Delete Attributes, Rename Attributes,
# Fast Access to Features at ID, Presimplify Geometries,
# Presimplify Geometries with Validity Check, Transactions, Curved Geometries

Could you please guide me on how to change the capability of the layer to include changing attributes and geometry?


texts: Hopefully this isn't too dumb a question, but is it possible to define projection at the record level, instead of a whole layer or table? I want to create a "project" database for all the work my company has done, and currently have a table of lat longs. It'd be nice to somehow create a dataset of actual site boundaries (usually in UTM NAD83 projection, various zones), without having to convert everything over to lat long. Ideas?


texts: Is there an online catalog for declassified Soviet Union Spy Satellite images (Zenit-4, Resurs-F1, Okean)?

texts:
I'm a glaciologist user of declassified US Spy satellite images (Corona Key Hole satelites, Hexagon, etc.), that provide an extraordinary record of glacier extent at resolutions equivalent to about 4 meters per pixel (see this amazing example). There are also ex Soviet Union declassified images from the '60s onward from satellites like Zenit-4, Resurs-F1, Okean-O1, etc. with similar technical specification than the US images or maybe better (they had color and NIR capability). However all my attempts to access such images end up in pages written in Russian that I can't understand or navigate trough. Some interesting information about sources of imagery is here (search for "Subject: Russian Imagery"), but with my null knowledge of Russian I haven't been able to find where to order such imagery, not an online catalog, nor to find if such catalog even exist.

So: Does anyone know how to access/order declassified ex Soviet Union imagery?

[EDIT January 25 2018: I've added and answer below with information I found regarding an archive containing part of the images but not all.]

As an example, here is a 1973 Zenith image taken over antarctica:

And here is a crop of the above image compared with Landsat 7 ETM+ at 30 meters per pixel. However, according to specifications, properly scan film should provide a much sharper and detailed image:


texts:
I have some .tif world image with 1km x 1km spatial resolution. I also have .shp file with 2km x 2km grid data. I need to convert .shp file into .tif in the same spatial resolution and projection as .tif world image.
Is it possible to process in QGIS?


texts: Creating a random GDAL raster with random predefined values

texts:
So far I have developed the below code but for some reason it only works for say one class of land-use instead of two. Would you please tell me where I have gone wrong?

import os, numpy, osgeo.gdal, math, random, sys
from osgeo.gdal import *
import  osgeo.gdalconst, osgeo.ogr
osgeo.gdal.AllRegister()
driver = osgeo.gdal.GetDriverByName('RST')
driver.Register()
ds = osgeo.gdal.Open('landii.rst')
cols = ds.RasterXSize
rows = ds.RasterYSize
ds.GetProjection()
print 'no of cols is:', cols
print 'no of rows is:', rows
print 'Raster has %s bands' % ds.RasterCount
geotransform = ds.GetGeoTransform()

out = driver.Create ('chandal.rst', cols, rows, 1, osgeo.gdal.GDT_Byte)
driver = ds.GetGeoTransform
raster = numpy.zeros((rows, cols), dtype=numpy.int32)

for i in range(rows):
    for j in range(cols):
        coll = 0
        raster[i][j] = coll
out.GetRasterBand(1).WriteArray(raster)
perc = [0.65, 0.45] # perc is the percentage of the pixels to be allocated to land use 1 and two respectively

nclasses = 2
for ii in range(nclasses):
    l = perc[ii]
number = int(cols * rows * l)
while (number > 0):
    rand1 = random.randrange(rows)
    rand2 = random.randrange(cols)
    k = raster[rand1][rand2]
    if k == 0:
        raster[rand1][rand2] = ii
        number = number - 1

out.GetRasterBand(1).WriteArray(raster)

out = None


texts: Cannot See a Web Map on ArcGIS Collector

texts:
I have created a web map and it appears that all settings and permissions should allow me to see it.  While I have other maps which I can see in Collector, when I open collector and sign in to this organization account I cannot see it.  The settings seem to be the same as other maps which I can see.

Here are the settings for the feature layer that I am trying to collect data for:

Any other suggestions to consider?


texts: Map extent matching permanently main canvas extent

texts:
Refreshing the map extent on QGIS's composer requires to click on Set map extent to match main canvas extent

I was wondering if there was a possibility to have a map extent (in the composer) that constantly refreshes itself depending on the main canevas location. That means, that no manual action has to be done.
In comparison to ArcGIS and if I'm not mistaken, the solution is symbolised by a locker located right bottom corner (at least for ArcMap 10).


texts: Overwrite existing Map Service through ArcPy keeping the existing cache tiling scheme

texts:
This is a follow-up question to Overwrite existing map service in ArcGIS Server 10.1 (through ArcPy)

The issue arises when I try to overwrite a map service that has a cache already defined.

If I overwrite the service through ArcGIS Desktop the existing cache tile scheme is recognized and preserved through the overwrite.
On the contrary, if I redeploy the service through Arcpy, that is create the .sddraft, transform it in .sd and upload it:

arcpy.mapping.CreateMapSDDraft(mapDoc, sddraft, service, 'ARCGIS_SERVER',
                                      con, True, None, summary, tags)
arcpy.StageService_server(sddraft, sd)
arcpy.UploadServiceDefinition_server(sd, con)

the existing tiling cache schema is not preserved and it's overwritten with the Bing standard cache tile scheme (I can see it in the generated .sddraft file).

Is there a way to preserve the existing schema?


texts: Working with Street Centerlines and a Address Point Layer

texts:
My Task was to update the Street Centerline Layer for our county and it wasn't snapped to intersection to intersection and I merged all that together. We at the county here keep a updated address point layer and I want to have each corner address show on the Street Centerline data.

I'm using ArcGIS 10.4.1 with the Advanced License without any extensions added.


texts: GeoServer WMS layer doesn't see file content

texts:
I've created a new layer, and styled with an .sld file (categorized by value). It works, if I check with GeoServer OpenLayers Preview, but when I check for this in a .html it doesn't appear. Is the problem with the bounds? I copied it from the URL:

...wms?service=WMS&version=1.1.0&request=GetMap&layers=foldrengesekGIS:mercatorvetulet&styles=&bbox=139.0697,36.60050000000007,143.85570000000007,41.53800000000007&width=496&height=512&srs=EPSG:3857&format=application/openlayers

Actually there is only one data, but the place is also wrong.


texts: Why doesn't this expression filter in my Value Relation?

texts:
Why won't this filter expression work?
I have followed all the examples I can find but it allows no items to pass through the filter.


texts: How to add images containing attributes in QGIS?

texts:
I am new to QGIS.  I know how to add vector and raster files into QGIS.  I have some satellite images taken at different times of the same location.  In addition, there are some annotation data associated with each file.

I also find the eVis plugin can integrate images on locations of any vector file.  The eVis display window, however, is small and cannot be zoomed in/out.  Does anyone know any way that I can add my images into QGIS on the "Map View" area and be able to edit my attribute information for each image?

Thanks!


texts: I'm working with ArcGis 10 and I need to show the name of the streets along the street line. I want to increase the distance of the label from the line and I saw in the Label Tab I can control this by changing the offset in Placement Properties, but I fear this option does not work when the feature is in a GCS. Is that right?


texts:
Since QGIS 3.6 there has been introduced a new script style. I am fairly new to this scripting and I would like to ask for the code to do an fix of two layers and after this to make an intersect of the fixed layers.

   from qgis.processing import alg
   from qgis.core import QgsFeature, QgsFeatureSink
   @alg(name="ex_new", label=alg.tr("Example script (new style)"), group="examplescripts", group_label=alg.tr("Example Scripts"))
   @alg.input(type=alg.SOURCE, name="INPUT", label="1. itsLayer")
   @alg.input(type=alg.SOURCE, name="INPUT2", label="2. itsLayer")
   @alg.input(type=alg.SINK, name="OUTPUT", label="Output layer")
   def testalg(instance, parameters, context, feedback, inputs):
       """
       Intersect two layers und make some SQL-Statements with the result
       """
       source = instance.parameterAsSource(parameters, "INPUT", context)
       (sink, dest_id) = instance.parameterAsSink(
           parameters, "OUTPUT", context,
           source.fields(), source.wkbType(), source.sourceCrs())
       total = 100.0 / source.featureCount() if source.featureCount() else 0
       features = source.getFeatures()
       for current, feature in enumerate(features):
           if feedback.isCanceled():
               break
           out_feature = QgsFeature(feature)
           sink.addFeature(out_feature, QgsFeatureSink.FastInsert)
           feedback.setProgress(int(current * total))
       return {"OUTPUT": dest_id}


texts: Issue to represent the palette in Google Earth Engine

texts:
I am trying to represent a raster giving to the palette the "max" and "min" of the raster according to the actual raster values. Instead of giving a minimum value of -1 and maximum value of 1 for the NDVI, I want that the minimum real value of the image represents the minimum value of the palette and the same for the maximum value. For instance, if the minimum value of a NDVI image is -0.7 and the maximum is 0.8 I want those values as "min" and "max" as parameters for the palette.
I thought that looking for the minimum/maximum value of the image with a reducer and transforming them to a Number I would be able to do it. However, when I run the code, I get the following error:
NDVI: Layer error: Can't encode object:
abs() Computes the absolute value of the input.
Args:this:input (Number): The input value.
The code is the following:
// Load image.
var img = ee.Image('MODIS/006/MOD09GA/2012_03_09');

// Use the normalizedDifference(A, B) to compute (A - B) / (A + B)
var ndvi = img.normalizedDifference(['sur_refl_b02', 'sur_refl_b01']);

var Min = ee.Number(ndvi.reduceRegion({
     reducer: ee.Reducer.min(),
     //scale: 250,
     maxPixels: 1e9,
     geometry: geometry,
     }));
  print(Min, 'NDVImin');

var Max = ee.Number(ndvi.reduceRegion({
      reducer: ee.Reducer.max(),
     //scale: 250,
     maxPixels: 1e9,
     geometry: geometry,
     }));

print(Max, 'NDVImax');

var BandCompViz = {min:Min, max:Max, palette:
  ['ffffff','b7f0ae','21f600','0000FF','FDFF92','FF2700','d600ff']};

// Center the map
Map.setCenter(-94.84497, 39.01918, 8);

Map.addLayer(ndvi, BandCompViz, 'NDVI');

https://code.earthengine.google.com/384e02f60104985848eb90fac875dfa6


texts: Finding shortest routes between multiple origins and destinations?

texts: For example, I have 100 origin points and 100 destination points. All points are connected to each other. I want to know the route between each origin and each destination. How do I do it all at once?


texts:
I am trying to develop an application, firstly web, later mobile, that builds a route based on points of interest, basically it should show the route and the points of interest and when the user clicks it will display some kind of information but I'm really struggling starting to create this...I already have my POI spatial database and I haven't really tried anything, because I don't know where to start. Is there any example you can give me?

Thanks in advance.


texts: What changes will be made to arcpy in ArcGIS 10.1?

texts: I heard that two big changes to arcpy will be add-ins and new, faster cursors. Also, new features include "the ability to administer your enterprsie geodatabase (ArcSDE) using Python scripting" and "Dynamic legend creation." I was wondering if anyone knew of other specific changes, maybe something overheard at the 2011 user conference. A changelog of all the changes would be ideal, but I haven't had luck finding anything.


texts: Create points with attributes from Cursor

texts:
I have a table imported from Excel that has three attribute columns (City, State, ZIP) as well as X/Y columns that have been added.

I have a cursor that goes through the records in the table and creates points

point = arcpy.Point()
pointGeometryList = []

with arcpy.da.SearchCursor("in_memory/xxExcelOutput", ["POINT_X", "POINT_Y", "CITY", "STATE", "ZIP_CODE"]) as cursor:
    for row in cursor:
        point.X = row[0]
        point.Y = row[1]

        pointGeometry = arcpy.PointGeometry(point)
        pointGeometryList.append(pointGeometry)

arcpy.CopyFeatures_management(pointGeometryList, "in_memory/pointFC")

How can I add the attributes from my table when I create the points?  I know the following doesn't work, but is there a way to do it within the same cursor?

    for row in cursor:
        point.X = row[0]
        point.Y = row[1]
        point.ZIP = row[4]

The only way I've come up with so far is to AddField() to the new feature class, and then UpdateCursor() the extra attributes with another search cursor.  While this way works OK for my test data (I have some data in my table, for proof of concept), I am concerned that I may end up with multiple points in the same location that have different attributes - I need to ensure that each point is created with the matching attributes.


texts: Copying and renaming files on SDE connection

texts:
I have a script that renames and copies files (feature classes) on an SDE connection. This script works just fine if the feature class name does not start with an 'S'. If the feature class name does start with an 'S', it gets removed. As you can see in the script below I use Python's lstrip to remove the 'SDE.' from the file name to rename it. The 'SDE.' gets added again automatically after copying. So a feature class named HAM would be renamed SDE.R_HAM_QC. If the feature class is named STREETS it get renamed to SDE.R_TREETS_QC. See how the first S in STREETS is removed? Does anyone have a possible solution?

    #Import modules
    import arcpy, os, sys
    from arcpy import env
    arcpy.env.overwriteOutput = True

    #Set variables
    outWorkspace = "Database Connections\Connection to 10.2.5.468.sde"
    arcpy.env.workspace = outWorkspace
    fileName = "SDE.STREETS"

    #Describe the file name
    desc = arcpy.Describe(fileName)

    if desc.datasetType == "FeatureClass":

          #Rename and copy data if it is a feature class
          name1 = fileName.lstrip('SDE.')
          print name1
          newName = "R_" + name1 + "_QC"
          outFile = os.path.join(outWorkspace, newName)
          arcpy.Copy_management(fileName, outFile)


texts: Converting shape layer with POLYGON and MULTIPOLYGON to POINT geometry

texts:
I'm trying to convert shape file to POINT geometry using QGIS 3.10.5.

Here is info of my shape file:

Here are my conversion settings:

Error that I get:

  Export to vector file failed. Error: Feature write errors: Feature
  creation error (OGR error: Attempt to write non-multipoint
  (MULTIPOLYGON) geometry to multipoint shapefile.)

I can successfully load data in MSSQL server which is recognizing all shapes nicely so shape file is OK.
I need to convert to POINT geometry to work with visualizations in Tableau.
Any help or pointers is highly appreciated.

EDIT:
My actual problem was that I had MULTIPOLYGONS which Tableau doesn't support. Tableau do support POLYGONS.


texts: How to fill fields with layer name in PyQGIS

texts:
I have hundreds of shapefiles that I have added to a QGIS project. I am using QGIS 3.0.0. I have created a new text field for each of them named 'Enumertr' in PyQGIS and I need to populate this field in each of the shapefiles with the name of the respective layer. I have the below script, but it is bringing an error. I also need to adapt it to iterate through all of the layers in my QGIS project. I am relatively new to Python but am hoping using PyQGIS will save me a substantial amount of time if I can finalise this script.

from PyQt4.QtCore import QVariant
from qgis.core import QgsField, QgsExpression, QgsFeature
vl = iface.activeLayer()
vl.startEditing()
idx = vl.lookupField(‘Enumertr’)
e = QgsExpression (‘vl.name’)
e.prepare(vl.fields())
for f in vl.getFeatures():
    f[idx] = e.evaluate( f )
    vl.updateFeature( f )
vl.commitChanges()


texts: Clipping features with the attributes of another feature

texts:
I have a vector polygon for a catchment which contains numerous subbasins. In addition, I also, possess another vector polygon showing the impervious areas of the catchment.

What I want to do is determine just how much of each subbasin is impervious. To do so, I planned to clip the impervious area vector with the boundaries of the subbasins contained by the cat catchment vector. How can this be achieved?

I've tried using the clip function with both datasets as they were, but that hasn't worked. I also attempted to use the clip function with the vector files converted to polyline, and still not the result I was looking for.

I'm using ArcMap 10.3.


texts: Reprojecting point data from EPSG:3078 to EPSG:4326

texts:
There are services to project from one coordinate system to another, GIS software packages that perform projections, but it is all black box, as far as I've been able to find.

I want to see the calculations/equations necessary to calculate the reprojection from NAD 83 Michigan GeoRef (EPSG:3078) to WGS84 (EPSG:4326). Does anyone know where I can find the equations?

I am aware, that one is a Projected Coordinate System and the other is a Geographic Coordinate System, so my use of "project" may not be the correct word choice.


texts: Analyzing Shapefile and Checking if particular coordinate is located in a region in shapefile

texts:
I am writing a Simulink/Matlab Application that takes in a shapefile. When given a coordinate(lat,long), I want to see if the coordinate is inside a polygon in shapefile in real-time.

My current option:
1) Use OGR tool to read the shapefile and use its feature deterimine if pointi in polygon. I do not know if it is easy to integrate OGR library into Simulink environment.

2) Use Python library to first read shapefile. Pass the read information into Simulink application to perform points in polygon tests.


texts: T-SQL with SDE and Web Map ArcGIS Server 10 MSSQL 2008

texts:
I need to build a system whereby the user of a web map (the standard web map that ArcGIS Server Manager builds) can edit a feature class's attributes.  The attribute will be a 'Comments' field.  The feature class will be replaced quarterly but the 'Comments' field must persist across quarterly updates.

I am thinking to build a stand alone table in sde (MSSQL 2008) and let users of the web map update the feature class's 'Comments' field directly then have a db trigger fire which will write to the stand alone table.. the unique id of the record in the feature class and also the comments edit.  This would be to retain the 'Comments' field apart from the feature class thus persisting the comments across quarterly updates of the feature classes.

I have already built the standard web map from AGS Manager and created the feature service of the feature class and tested the editing of a 'Comments' field.  Now I am starting to build the persisted 'Comments' architecture.

Also, there may be multiple 'Comments' attributes inserted over time, I am uncertain as to whether I should have multiple 'Comments' fields in a stand alone table OR have multiple records in the stand alone table with each comment and the unique id and do a one-to-many join to return them in the web map.  I also think deleting comments will be important.

Does anyone have any suggestions or tips or alternative approaches?  If the T-SQL approach, are there any tips on T-SQL and SDE?

Thanks.


texts: How to compare string fields in the same table with ArcGIS 10

texts:
I don't know exactly if there is a tool for this or a maybe a python script could be more suitable.
I want to compare string elements from a field (POI_NAME) with another string element in another field (ROTULO). Find attached an screenshot.

My goal, AT LEAST, is to have a YES/NO matching, but It could be fantastic having a % of matching. Do anyone know anything about?

I tried something like this:

But it doesn't work.

After editing the code, to:

We can see that it does not working well (BECAUSE IT HAS SPECIAL CHARS AND SPACES):


texts: Preview-problem when editing small polygon-features in QField

texts:
I have the following problem.
When surveying small objects on QField, such as structures with a diameter of 20 to 100 cm (with external GNSS-receiver, +- 1 to 3 cm accuracy),  a problem with the preview became apparent.
The preview is not displayed correctly during the measurement and the shape of the preview does not correspond to the object to be measured.
But after saving the object and finishing the editing process, it will be displayed correctly on the map.
The same problem reappears in the post-processing of the object, for example when I would like to edit the vertices again (see image below).
I suspect that Qfield has maybe problems with displaying the preview in larger scales, because the smaller the scale and the larger the object being edited, the better the preview image.
Unfortunately, this problem makes the measurement of smaller areas more difficult,
since during the measurement process, you have no control over whether the measurement is correct or not.
First I thought this might be a problem with the GNSS-signal, but this works fine. And after finishing the edit-process the feature looks correct. The problem also occurs, when the GPS-tool is turned off and a small polygon is drawn by hand (like in the image below). I also tested other polygon-files and formats for example the "simple_form"-layer in the tutorial.
Is anyone aware of this problem or knows how to fix it?
The project:

we try test qfield on archaeological site for the documentation-process
among other things, even small structures such as pits or post-holes need to be measured (20-100 cm and bigger)

Here some technical information:

tested with QGIS-Version 3.10 and 3.16 and QField-Version 1.7.6 (on Android-Version 8 Smartphone and Tablet)
for the measurement in the field we use a Leica Zeno gg04, which works together with an android-tablet,
the GPS-signal comes via the "Zeno-connect"- app which is activated in the mock-location


texts: Umlaute not displayed correctly when loading .mdb MS Access into QGIS (1.9.0-Master) with eVis

texts:
I am loading an .mdb (MS Access 97 or 2003) to QGIS (1.9) via the eVis database conection. Connecting to the mdb,  selecting a table and displaying it as a point-layer through the appropriate x/y-fields works fine.

However in the attribute table the German Umlaute ä, ö, ü and also ß are displayed as inverted "?". After days of searching I don't really know how to tackle this one. I am grateful for any recommendation on how to fix this.

The umlaute in the mdb / attribute table are displayed correctly in ArcMap. I would love to get this running in QGIS.