texts: QGIS 2.14.2 Difference Tool Returning blank layer

texts:
I'm using the Vector Geoprocessing/Difference Tool to cookie-cut out of one polygon using several polygons from another layer, but it is returning a blank layer.  I note this query has been asked before and the problem can be that the layers have different CRS's.  That is not the case here.  I have uploaded a sample of the data set here:

https://www.dropbox.com/sh/rovtjvngrzow82z/AABG3yM9o3-JtRwuc2Wfjld2a?dl=0

I want the 'building footprints' cut out of the 'negative area' polygon.


texts: GeoJSON Layer - ArcGIS JS API - External URL with oAuth2.0

texts:
I am using ArcGIS JS API and connecting to an external GeoJSON service for my GeoJSON Layer. However, for security reasons oAuth2.0 has been enabled on the service.

So now I am lost as to how to connect to this service.

esri/identity/OAuthInfo has URL option only for ArcGIS Portal URL.

I believe it should be a combination of esri/identity/IdentityManager and esri/identity/Credential but it would be great if there is some step-by-step guidance.


texts:
I have this problem with my images of the Ross Ice Shelf:

This is an image that has been Radar - Radiometric- Calibrate- Amplitude_HH - Gamma0
followed by Speckle Filter.
This is what happens when I add a Geometric - Geolocation - Ellipsoid Correction - Geolocation Grid

This is what happens when I add Geometric - Geolocation - Ellipsoid Correction - Average Height Range Doppler

This is what the image looks like on Copernicus preview:

Any suggestions? I'm at a loss


texts: Extracting only some properties of each feature from a feature collection and creating a new a feature collection from them using Google Earth Engine

texts:
I have a feature collection but I need to make a new feature collection from it by taking only some columns , I am new to google earth engine and am unable to do so ,the code is

NDVI=NDVI.map(function(
    f=ee.Feature(f)
    return ee.Feature(f.toDictionary(dates).values(dates))
 })
print(NDVI)

here NDVI is the main feature collection and
dates is a list of properties that I want to extract per feature, the error I am getting is

FeatureCollection (Error)
Collection.map: A mapped algorithm must return a Feature or Image.


texts: Store the filenames of .kml files into the attribute table of a shapefile

texts:
I have to be joining .kml files in a single shapefile. The problem is that the kml files comes with a name that when I convert it into a shapefile it is not kept in the attribute table. thus, I lose a lot of time by filling the attribute table manually.
I would like that when I receive the files, 1234.kml
and 2345.kml, the shapefile's attribute table appears with "1234" and "2345".


texts: Fragstats input data for whole study area

texts: I have a land cover map covering an area of 1560 square km. I am interested in landscape fragmentation analysis using Fragstats. But I am confused about input data for the whole study area. I do not have problem on data preparation but problem about whole study area. I am not sure about how can I use whole study area as input data. Could anyone please let me clear about how can I use my whole study area for the analysis?


texts: I am displaying three different layers on a map, essentially three separate habitat types. Each have their own colour.  In cases, 2 or 3 of the layers overlap in part (essentially showing where habitat types overlap).  How can I display this overlap in the legend?


texts: Formating GeoJSON file with Leaflet Ajax-plugin

texts:
First of all I apologize for my basic english. I have to say, too, that I'm really a newby in coding.
I'm trying to create an interactive map using Leaflet (you can see it here). I was able to show my geojson file using Ajax-plugin. Now I want to change the style of my GeoJSON file, but I really don't know how I can do that. Can someone explain me?

This is my code:

<!DOCTYPE html>

<meta charset=utf-8 />
<title>Geojson data</title>
<meta name='viewport' content='initial-scale=1,maximum-scale=1,user-scalable=no' />

            <link rel="stylesheet" href="http://cdn.leafletjs.com/leaflet-0.7.2/leaflet.css" />

Thanks


texts: Coordinate drift when using osm2pgsql and specifying EPSG code

texts:
We had a PostGIS database which was initialized quite some time ago using an osm.pbf file from Geofabrik, without a specified coordinate system. It defaulted to using EPSG 900913 (Google Maps Global Mercator)

Recently, we initialized a new database the same way, but this time we used the flag -E 31370, which initalized the data to EPSG 31370 (Belge 1972 / Belgian Lambert 72).

When we compare the two databases, we can see a very clear coordinate drift (see pic)

Red = EPSG 31370 (Belge 1972 / Belgian Lambert 72)

Blue = EPSG 900913 (Google Maps Global Mercator)

Geodesy is really not my specialty. A colleague and I spent some time going over it trying to figure out which one is right.

Obviously the blue polygons are in the right position according to the basemap. But the basemap is mostly likely also derived from a PostGIS database that has been initialized to a Web Mercator projection for the sake of visualization, whereas the database for the red polygons was initialized to a local projection immediately.

It is my understanding that the OSM dumps from Geofabrik are in WGS84 natively so I am confused as to why they are ending up in two different places.

This is quite a big issue for us because a previous analysis workflow did all the coordinate transformations manually, while a new workflow is trying to initialize everything in the same coordinate system and just use it for the entire duration.

How can we stop this coordinate drift?


texts: Getting the elevation of 400 million points in the US, how to get the elevation of these points from a dataset?

texts:
I'm a software engineer working on a game concept that would create the United States in the form of a grid, each grid being 200m x 200m in size with an elevation attached to each corner to create a polygon-looking isometric view. Regardless, I have been looking into the USGS datasets for elevation, but they're currently very fragmented, as there are probably about 2,000+ datasets for the nation.

My question is: How would I be able to get the elevation for these territories? It would be around 400 million points, so I need to have the data available on my local instance. Is there a national data file?


texts: Measuring distance from observations to coastline

texts:
I am looking for a way to use a computer software to measure distances from observation points to the coast of Iceland. I have over 51,000 observations as a points layer and I have the coastline of Iceland as a polygon. Is it possible to do this calculation in QGIS and add a column (Distance) to the attributes table of the points layer?
I've found no instructional YouTube video to solve this particular exercise, though I've found many that explain how to measure distances between points. I've also found one thread where the same question was asked here last year, but the only answer was that using a different program would be the only way. The person asking that question has the same problem I do: I don't know how to code, so PostGIS is not a viable option.


texts: Latitude and longitude from .shp file with QGIS

texts:
I opened with QGIS the shapefile in the following archive  http://data.openstreetmapdata.com/coastlines-split-4326.zip and I want to convert this shapefile to a list of latitude and longitude values for every coastline point. To do so, I opened the .shp file with QGIS, then clicked on 'lines' on the layer panel, selected 'Save as...', selected the following options , and clicked on 'OK'.

I obtain a .csv file that looks like this

"WKT,FID,
"LINESTRING (-4.94237 55.725449,-4.941922 55.725585, ... "

May you please confirm that the pairs of numbers (e.g. -4.94237 55.725449) correspond to the GEODETIC latitude of a point and to its longitude (both in degrees), respectively?


texts: Make a raster with pixel scale 100 by 100 meter with gdal_rasterize

texts:
I'm asking here, because, I've got great results last time. Here it is the next of my adventure with gdal_rasterize, I want to make a raster picture where each pixel correspond to a square of 100 by 100 meters from a table of points of a Postgis DB (EPSG:3857).

From the manual : http://www.gdal.org/gdal_rasterize.html I found the option -tr normally is made to do it. I'm executing the following command :

gdal_rasterize -a value -tr 100 100 PG:"host=localhost dbname=myDB user=postgres" -sql "SELECT the_geom, value FROM table" out.tif

The resulting tiff is an irregular raster : sometime pixel size correspond to 50m sometime 150m (I calculate that thanks to the rule of Google earth)

What I did wrong ?

[Edit] As asked in comment, here it is the content of gdalinfo of the resulted tiff.

gdalinfo content for the file generated with the option "-tr 50 50" :

Driver: GTiff/GeoTIFF
Files: out_50_50.tif
Size is 22018, 15699
Coordinate System is:
PROJCS["WGS 84 / Pseudo-Mercator",
    GEOGCS["WGS 84",
        DATUM["WGS_1984",
            SPHEROID["WGS 84",6378137,298.257223563,
                AUTHORITY["EPSG","7030"]],
            AUTHORITY["EPSG","6326"]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433],
        AUTHORITY["EPSG","4326"]],
    PROJECTION["Mercator_1SP"],
    PARAMETER["central_meridian",0],
    PARAMETER["scale_factor",1],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    EXTENSION["PROJ4","+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +
x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext  +no_defs"],
    AUTHORITY["EPSG","3857"]]
Origin = (-1045006.900697760100000,5457301.723357450200000)
Pixel Size = (50.000000000000000,-50.000000000000000)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (-1045006.901, 5457301.723) (  9d23'14.84"W, 43d56'50.54"N)
Lower Left  (-1045006.901, 4672351.723) (  9d23'14.84"W, 38d39'17.46"N)
Upper Right (   55893.099, 5457301.723) (  0d30' 7.55"E, 43d56'50.54"N)
Lower Right (   55893.099, 4672351.723) (  0d30' 7.55"E, 38d39'17.46"N)
Center      ( -494556.901, 5064826.723) (  4d26'33.65"W, 41d21'17.53"N)
Band 1 Block=22018x1 Type=Float64, ColorInterp=Gray

gdalinfo content for the file generated with the option "-tr 150 150" :

Driver: GTiff/GeoTIFF
Files: out_150_150.tif
Size is 7339, 5233
Coordinate System is:
PROJCS["WGS 84 / Pseudo-Mercator",
    GEOGCS["WGS 84",
        DATUM["WGS_1984",
            SPHEROID["WGS 84",6378137,298.257223563,
                AUTHORITY["EPSG","7030"]],
            AUTHORITY["EPSG","6326"]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433],
        AUTHORITY["EPSG","4326"]],
    PROJECTION["Mercator_1SP"],
    PARAMETER["central_meridian",0],
    PARAMETER["scale_factor",1],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    EXTENSION["PROJ4","+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +
x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext  +no_defs"],
    AUTHORITY["EPSG","3857"]]
Origin = (-1045006.900697760100000,5457301.723357450200000)
Pixel Size = (150.000000000000000,-150.000000000000000)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (-1045006.901, 5457301.723) (  9d23'14.84"W, 43d56'50.54"N)
Lower Left  (-1045006.901, 4672351.723) (  9d23'14.84"W, 38d39'17.46"N)
Upper Right (   55843.099, 5457301.723) (  0d30' 5.93"E, 43d56'50.54"N)
Lower Right (   55843.099, 4672351.723) (  0d30' 5.93"E, 38d39'17.46"N)
Center      ( -494581.901, 5064826.723) (  4d26'34.46"W, 41d21'17.53"N)
Band 1 Block=7339x1 Type=Float64, ColorInterp=Gray

Is my projection well suited ?


texts: Join/relation not an available widget option in the attributes form?

texts:
I cannot figure out why, but, in a tree/planting projects database, I have a fairly simply data model where I have one polygon layer to symbolize the area/project (eg. Maple Oak Park, John Smith School), and a point layer to symbolize the trees planted in that area/project. I created a relation between the two layers using the area/project name as the join field (image #1).
I want to then add this relation to the form for the area/project layer, so that all the trees associated to the area/project appear on the form. See image #3 for an example where this worked for me a few months ago on another set of layers. However, for some reason the relation simply does not appear as an option in the attributes form (image #2).
https://imgur.com/a/ovPhWqw


texts: Validate that Openlayers GeoJSON layer is in the correct format with OpenLayers.Protocol.HTTP call

texts:
I am completely new to JavaScript and OpenLayers, I am trying to figure out a way to validate the received layer with the OpenLayers.Protocol.HTTP call. I want to validate that it is in the correct format (GeoJSON). I have tried callback function, and it doesn't seem to work. Also I have tried through OpenLayers. Request.GET method but I can't really notify the user when the query has given an error. I can see that the GET method only validates if it has established a connection but not if the GeoJSON vector has the correct format.

In the test.php script called in the URL parameter, the output in case of a successful query is a functional GeoJSON layer.

What I am trying to achieve is when the query is not successful, (has returned 0 rows).

None of the alerts I have tried so far have been useful.

        function handler(request) {
        alert(request.responseText);
        alert(request.status);
        alert(request.getAllResponseHeaders());
        var featurecollection = request.responseText;
        var geojson_format = new OpenLayers.Format.GeoJSON();
        var geojson_vector = new OpenLayers.Layer.Vector();
        geojson_vector.addFeatures(geojson_format.read(featurecollection));           alert(geojson_format.isValidType(request.responseText,"FeatureCollection"));
        alert(geojson_format.extract.collection);
        if (geojson_format.read(request.responseText)) {
            alert("CORRECTO. SI SE PUDO CONECTAR A " + geojson_layer.protocol.options.url);
        }
        else{
            alert("ERROR. NO SE PUDO CONECTAR A " + geojson_layer.protocol.options.url);
        }
        }
        function refresh(proy) {
            var request = OpenLayers.Request.GET({
            url: "test.php?proyecto="+proy,
            format: new OpenLayers.Format.GeoJSON(),
            callback: handler
            });
            if (geojson_layer) {
            geojson_layer.protocol.options.url="test.php?proyecto="+proy;
                geojson_layer.refresh({force:true});

            }
            else {
                geojson_layer = new OpenLayers.Layer.Vector("Consulta", {
                    styleMap: new OpenLayers.StyleMap({
                        "default": new OpenLayers.Style({
                        pointRadius: 15,
                        fillColor: "#ffcc66",
                        fillOpacity: 0.2,
                        strokeColor: "#ffff00",
                        strokeWidth: 2,
                        strokeOpacity: 0.8 } ),
                        "select": { fillColor: "#8aeeef",
                        strokeColor: "#32a8a9",
                        labelYOffset:13,
                        label:"${name}"} //Text entspricht feature.attributes.name
                    }),
                    projection: new OpenLayers.Projection("EPSG:4326"),
                    strategies: [new OpenLayers.Strategy.Fixed(), new OpenLayers.Strategy.Refresh({force:true,active:true})],
                    protocol: new OpenLayers.Protocol.HTTP({
                        url: "test.php?proyecto="+proy,
                        format: new OpenLayers.Format.GeoJSON()
                    })
                });
                map.addLayer(geojson_layer);
            }
        }


texts: Question about the TableInfo function in MapBasic

texts:
I have a problem I really can't seem to get my head around, MapInfo crashes when I try to run and if statement.

I have pretty much exactly the same code in another sub to check for other things and that works perfectly fine. It crashes on the zlRow = TableInfo("zerolengths",TAB_INFO_NROWS), I've tried this with zlRow = TableInfo(0,TAB_INFO_NROWS) as there will only be one TAB file open when this rool is running. I just don't know why it works for other subs and queries but not this!!

Sub LenSub

nMID = FrontWindow()

nLayers = MapperInfo(nMID, MAPPER_INFO_LAYERS)

sLayer = LayerInfo(nMID, nLayers, LAYER_INFO_NAME)

Dim sDir as String

sDir = PathToDirectory$(TableInfo(sLayer,TAB_INFO_TABFILE))

    '###Check for blanks

Dim sLen as string

Dialog
    Title "Length"
    Control statictext
        Title "Enter length field"
    Control edittext
        Value sLen
        into sLen
    Control OKButton
    Control CancelButton

if commandinfo (CMD_INFO_DLG_OK) then

Else
End Program
End if

Run Command "select * from " & sLayer & " where " & sLen & " < 1 into zerolengths"

Dim zlRow as smallint

zlRow = TableInfo("zerolengths",TAB_INFO_NROWS)

if zlRow > 0 then
    Export Multipolylines Into sDir & "zerolengths.csv" Type "CSV"
    Delimiter ","
Titles

    Dialog
            Title "Zero lengths"
            Control statictext
                Title "Sections with zero lengths exported into zerolengths.csv"
            Control OKButton
Else
End If

End Sub

Thanks, Alan


texts:
Using QGIS to print a map but I cannot choose a size bigger than A4 in the map composer.

How do I print a map for example 841x1189 mm size ?


texts: how to create a point shapefile with azimuth and elevation?

texts:
I have downloaded a GPS almanac file which gives azimuth and elevation for each GPS satellite at a particular latitude and longitude.

I want to create a point file with azimuth and elevation angle which represents each GPS satellite and drape it on a DEM having 3D model of buildings to carry out line of sight analysis.

How can I define a point shapefile having azimuth and elevation angle in ArcScene?


texts: How can I create a multi-dimensional spectral plot using R?

texts:
I have a raster stack of several layers such as NDVI, Slope, texture, and TCT bands, and would like to create a multi-dimensional spectral signature plot from my ground thruth points representing different vegetation classes.

How can I create a multi-dimensional spectral plot using R based on my raster stack and my ground truth data?

Here is my data:

Stack<-stack(NDVI, NDVI_glcm$glcm_variance, Mean_slope, Sr_TCT)
str(GT_corr)
'data.frame':   38098 obs. of  3 variables:
 $ long  : num  237481 237481 237481 237761 237761 ...
 $ lat   : num  1212239 1212239 1212239 1212240 1212240 ...
 $ M.Info: Factor w/ 23 levels "Bas fond","Bas fond brulee",..: 18 18 18 19 19 19 19 19 20 20 ...

For example, the plot could look like something like this if I would consider only 3 of the layers from my stack:

where each point would represent a ground truth point for a given vegetation class, and a different color would be used for each class.


texts: Can you add a WMTS as a GeoServer Data Source?

texts: From this GeoServer tutorial, (under "Creating a store") I can only see Vector, Database and Raster data types as the only possible data sources. But I was wondering if a Web Mapping service such as a WMTS could be added as a source?


texts: Intersect_analysis in ArcPy "The table was not found."

texts:
Is there a bug in ArcPy?  I have two features I want to intersect in a feature data set using ArcPy (this is a small part of a bigger program).  Let me preface this by saying that when I do this in ArcMap 10.7 with the same two feature classes, it works perfectly.  I need it to work in my program.

The code leading up to the intersect statement is below. All the code leading to the intersect statement works and the file used to intersect the buffer layer is created, the buffer layer is also in the workspace.

import arcpy
import math

filepath = r'C:\Users\aschi\TestPurposes.gdb\TestFeatureDataSet'
arcpy.env.workspace = filepath
sr = arcpy.Describe("MMFSimplifyPt10").spatialReference

# Add field to keep track of turns
arcpy.AddField_management("MMFSimplifyPt10", "turn_count", 'LONG')

arcpy.CreateFeatureclass_management(filepath,'MMFSimplifyPt10_Angles',"POINT",spatial_reference=sr)
arcpy.AddField_management('MMFSimplifyPt10_Angles','Angle',"DOUBLE")
arcpy.AddField_management('MMFSimplifyPt10_Angles','Trip_ID',"DOUBLE")

# Initiate insert cursor to store turn angles
Fields = ['SHAPE@', 'Angle','Trip_ID']
insCur = arcpy.da.InsertCursor("MMFSimplifyPt10_Angles",Fields)

def GetAngle(A, B, C):
    distAB = B.distanceTo(A)
    distBC = C.distanceTo(B)
    distAC = C.distanceTo(A)
    try:
        cosangB = ((distAB * distAB) + (distBC * distBC) - (distAC * distAC)) / (2 * distAB * distBC)
        angB = math.acos(cosangB)
    except:
        angB = math.pi
        print 'Divide by zero error.  Coinciding points, going to next.'
    return angB

# Finds the angle at each vertex and stores those that are "small enough" to be possible turns.
with arcpy.da.SearchCursor("MMFSimplifyPt10", ["SHAPE@", "trip_id"]) as Scursor:
    for row in Scursor:
        for part in row[0]:
            pt_count = 1
            for pnt in part:
                pnt = arcpy.PointGeometry(pnt, sr)
                if pt_count > 1:
                    if pt_count > 2:
                        angle1 = GetAngle(twoBack, oneBack, pnt)
                        print angle1
                        if angle1 < 7*math.pi/9:
                            insCur.insertRow([oneBack, angle1, row[1]])
                    twoBack = oneBack
                oneBack = pnt
                pt_count += 1
del insCur

arcpy.Intersect_analysis(['MMFSimplifyPt10_Angles', 'TriCounty_Junctions_Buffer'], 'MMFSimplifyPt10_PTurns', output_type='POINT')

According to the documentation on intersect this is the correct syntax, the list is my two input features to be intersected, the second parameter is the output feature class to be created by "Intersect" and I want the output to be point data.

When I run the code I get the following error:

arcgisscripting.ExecuteError: ERROR 999999: Error executing function.
The table was not found.
The table was not found. [MMFSimplifyPt10_PTurns]
The table was not found.
The table was not found. [MMFSimplifyPt10_PTurns]
Invalid Topology [Topoengine error.]
Failed to execute (Intersect).

Process finished with exit code 1

Notice that the table that wasn't found is my output table that should be created by the Intersect tool.  I've seen other posts here where people got this error but it was their input features.  Also, I did try creating the output feature with a blank attribute table first and got the error that "MMFSimplifyPt10_PTurns" already exists.

Does anybody know of a way to fix this?  It seems like if it works in ArcMap it should work in Python too, but it doesn't.


texts:
I have one set of polygons and one set of points sharing a common attribute (value 0-503).

When applying style with new random color  ramp (Categorized) the first set of polygons is displayed correctly but when applying the same new random color ramp to the set of points colors are incorrect.

I tried to create an SVG colormap using "i want hue", and SV2ColoR plugin and it worked quite well (except that QGIS seemed to have problem with a large number of colors in the xml so I restricted to 100 in the svg in the end) but still wrong colors for the set of points. It looks likes something totally different is picked to display point.

Any ideas?


texts: Can QGIS generate custom formatted, nice looking maps?

texts:
I just finished a long trip, and with the gpx/kml track, I want to generate a nice looking map of our route. I only want to do some basic things like overlay the route over a basemap with significant locations and borders shown, show different days of the track by color, etc.

I'm new to GIS, but I set up QGIS and imported the gpx track no problem. Before I embark on learning enough GIS to accomplish the above though, I'm wondering if GIS can do this sort of cartography? Every GIS use-case I've heard is for some sort of map analysis, so I'm wondering if just generating a custom formatted and nice looking map is not in its realm of capability.


texts: Converting personal geodatabase .mdb file to .gdb

texts:
Is it possible to convert an MS Access personal geodatabase file to a SDE enterprise geodatabase file format? What kind of python scripting would it entail? I need to convert it in such way to maintain consistent geometric/spatial properties for an application.

I have looked for similar questions but to no avail, please direct me to one if I missed it!


texts:
I have a separate database supermarkets and their suppliers. I have as well database which suppliers provide products to particular supermarket.
I would like to create such a tool/solution where after clicking on the point on map (which is supermarket) I will receive a message with their suppliers. It could be short message and "flash" of the suppliers.

I am new in gis, does anyone have an idea how to solve it?


texts: Show all contour line labels in Layout in QGIS

texts:
I created a map with labeled contour lines. When I am now adding the map to the layout view not all labels appear (e.g. 2000, 2500 and 3500 meters has a label, but not 3000)

[.

In the original map the area with contour lines is a lot bigger, thus all contour lines show up how they should. How can I force QGIS to show the labels at a specific site in the layout without clipping my original file?


texts: Detecting an ArcGIS Split Rule / Ratio Policy

texts:
I have two identical layers in ArcGIS Desktop, one of them with two numeric fields that have a Ratio Policy applied to them.

Is there any way in ArcGIS Desktop to identify a field [or a set of fields] in a layer with a Ratio Policy?

If I get an ArcMap Document with these two layers, how would I know which layer has a field with Ratio Policy?


texts: Understanding UnicodeDecodeError: 'utf8' codec in ArcPy script?

texts:
I work with arcmap 10.3 and python 2.7.8. I have more than 500 shapefiles that located in many folders and subFolders. All Sub Folders are located in one large directory. I try with arcpy to detect all  shapefiles that have in their attribute table ,in field name "YEUD", the value 20. I search all shape files that begin with letters "mig". Finally i tried to print all the shapefiles that had been found with value 20 in it.
When i run this code:

import arcpy,os,fnmatch,unicodedata,codecs

rootPath = r"C:\Project\layers"
pattern = 'mig*.shp'
for root, dirs, files in os.walk(rootPath):
    for filename in fnmatch.filter(files, pattern):
        shp = os.path.join(root, filename)
        if arcpy.ListFields(shp, "YEUD"):
            print("{} has YEUD field".format(shp))
            with arcpy.da.SearchCursor(shp, ["YEUD"]) as rows:
                for row in rows:
                    if row[0] == 52:
                        print("{} has a record with YEUD = wanted row".format(shp))
                        break

i get an error when the python meet files and folders with right to left font:

UnicodeDecodeError: 'utf8' codec can't decode byte 0xe7 in position 23: invalid continuation byte

For completeness, i asked this question in https://geonet.esri.com/message/519769#519769
and marked it as correct answer for files and folders names that written in left to right fonts ,but when i run this code i get an error when the python meet files and folders names with right to left fonts.

In GeoNet i didn't receive helpful answer. I also searched answers in stackOverflow but didn't understand how to unicode the script.


texts:
I've seen this question several times, but I'm not 100% sure what needs to happen for my situation.  Quite often the answer is to install OpenGEO4W & GDAL, but it appears I already have GDAL.
I just did a fresh install of the all in one PostgreSQL and PostGIS windows installer.  PostGIS seems to be working quite fine.  I'm attempting to use the ogr2ogr utility to ultimately migrate geo tables into a SpatiaLite database.  I only found the OGR exectubles in one spot and that is in the pgsql bin folder.
ogr2ogr --version results in this:
GDAL 3.2.0, released 2020/10/26
To simplify this I can generate the error I'm getting by running ogrinfo and attempt to connect to my local pgsql host, like so:
\database\PostgreSQL\13\bin>ogrinfo PG:"host=127.0.0.1 port=5432 user='postgres' password='blablabla' dbname='postgis_30_sample'"

And the instant outcome is this:
FAILURE:
Unable to open datasource `PG:host=127.0.0.1 port=5432 user='postgres' password='blablabla' dbname='postgis_30_sample'' with the following drivers.
  -> ESRIC
  -> PCIDSK
  -> PDS4
  -> VICAR
  -> JP2OpenJPEG
  -> PDF
  -> MBTiles
  -> EEDA
  -> OGCAPI
  -> ESRI Shapefile
  -> MapInfo File
  -> UK .NTF
  -> LVBAG
  -> OGR_SDTS
  -> S57
  -> DGN
  -> OGR_VRT
  -> REC
  -> Memory
  -> BNA
  -> CSV
  -> GML
  -> GPX
  -> KML
  -> GeoJSON
  -> GeoJSONSeq
  -> ESRIJSON
  -> TopoJSON
  -> OGR_GMT
  -> GPKG
  -> SQLite
  -> ODBC
  -> WAsP
  -> PGeo
  -> MSSQLSpatial
  -> OpenFileGDB
  -> XPlane
  -> DXF
  -> CAD
  -> FlatGeobuf
  -> Geoconcept
  -> GeoRSS
  -> GPSTrackMaker
  -> VFK
  -> PGDUMP
  -> OSM
  -> GPSBabel
  -> SUA
  -> OpenAir
  -> OGR_PDS
  -> WFS
  -> OAPIF
  -> HTF
  -> AeronavFAA
  -> Geomedia
  -> EDIGEO
  -> SVG
  -> CouchDB
  -> Cloudant
  -> Idrisi
  -> ARCGEN
  -> SEGUKOOA
  -> SEGY
  -> XLS
  -> ODS
  -> XLSX
  -> Elasticsearch
  -> Walk
  -> Carto
  -> AmigoCloud
  -> SXF
  -> Selafin
  -> JML
  -> PLSCENES
  -> CSW
  -> VDV
  -> MVT
  -> NGW
  -> MapML
  -> TIGER
  -> AVCBin
  -> AVCE00
  -> HTTP

So, Is OpenGEO4W needed?  And what install?  It seems like the standard install is a large suite of tools?  It is confusing to me that the instance packaged with the PgSQL install doesn't have the driver necessary to connect to PgSQL.  Is there a way to simply add the PgSQL driver for ogr?


texts: Syntax to POST query to Soil Data Access Web Service

texts:
I want to query the SSURGO data through the "post.rest" service described here:

http://sdmdataaccess.nrcs.usda.gov/WebServiceHelp.aspx#PostRest

The way I'm reading the documentation, I think this should return a result:

curl -d "SELECT * from component where mukey = '458913' ORDER BY comppct_r DESC" http://SDMDataAccess.sc.egov.usda.gov/Tabular/SDMTabularService/post.rest

However, I get the error:

<?xml version='1.0' encoding="UTF-8" standalone="no" ?>
<ServiceExceptionReport xmlns="http://www.opengis.net/ogc" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" xsi:schemaLocation="http://www.opengis.net/ogc http://schemas.opengis.net/wms/1.1.1/OGC-exception.xsd">
<ServiceException>
Cannot parse ampersand-delimited POST request</ServiceException>
</ServiceExceptionReport>

The error doesn't make sense because my query doesn't have an ampersand.

I have confirmed that the query SELECT * from component where mukey = '458913' ORDER BY comppct_r DESC is valid by running it here:

http://sdmdataaccess.nrcs.usda.gov/Query.aspx

I tried a variety of different header and body variations, but nothing helped.

What am I doing wrong?


texts: Qneat3 / shortest distance calc sensitive to altitude?

texts:
I'm trying to calculate shortest distance using Qneat3 (using both the shortest path and OD-matrices for different purposes). And am wondering if there is any possobility to factor in altitide (it is a very montanous region) to acurately represent the distances travelled? or does the plugin already do this? seem not to find information on this anywhere in the ether..

Frida


texts: How to link data from a CSV table to a shapefile

texts:
I've recently started using QGIS and have been doing alright at self teaching myself. I have read a few similar questions however I have not got this to work for myself.

I have a shapefile which maps a variety of areas, I am trying to link numerical data inputted from a CSV table. The table has unique identifiers which can be linked to the shapefile.

I've tried using the add vector join function to link the data however when applied I can not find the linked data anywhere in the fields section.

I want to show the data from the table on the shapefile to show a gradient difference however the column from the CSV file does not show up on the shapefile properties.

Also, when I leave the layer properties section after joining the table, when I go back inot the properties it is no longer there?

Not sure if there is something simple I've missed?


texts: How to exaggerate the elevation in QGIS 3.6?

texts:
I´m using QGIS 3.6 to create a Thematic Map with a basemap (STRM).

So, to do that, I downloaded the STRM file (*.TIF) in this site: https://earthexplorer.usgs.gov/

File Name: s11_e040_1arc_v3

After the download I imported the two Layers (STRM.tif + DISTRICT.shp) to QGIS, and the result was this:

I´m not happy with the result because I would like to make the STRM more visible, note that the district.shp has 45% of OPACITY...

How can I make the map/STRM more expressive to see the elevations variation?

Do I have to exaggerate the elevation? if, yes how can I do?

NOTE: I´m working in 2D not 3D (I know how to do in 3D view).


texts: I suspect I'm missing something basic but when ever I sort a table with an integer field it ends up sorting 1, 10, 100, not 1, 2, 3.


texts: Installing ArcSDE engine plugin for Eclipse?

texts:
I'm looking for the following jars for the ArcSDE API:

jsde_sdk.jar
jpe_sdk.jar
concurrent.jar
icu4_3_2.jar

thanks to this answer I've found that they can be gotten installing the ArcSDE engine plugin for Eclipse.

But on ESRI's help page it calls for the following directory:

%ARCGISHOME%/java/tools/eclipse_plugin/arcgis_update_site/engine

My installation, instead, only has the following one:

%ARCGISHOME%/java/tools/eclipse_plugin/arcgis_update_site/arcobjects

That is, the arc Objects plugin, instead of the engine plugin.

Where (or how) can I get the ArcSDE engine plugin for Eclipse?


texts: Add a map element with a frame border and background

texts:
I am trying to insert a north arrow and a legend into the map with a white background frame and a border with a 3pt gap. The code below partially works, but it changes the background colour of the dataframe as well (which I don't want) and does not change the north arrow size to 125. The code below is just for the north arrow.

Private Sub NorthArrowAdder2(ByVal pArrow As INorthArrow, ByVal x As Double, ByVal y As Double)
    Dim activeView As IActiveView
    Dim gc As IGraphicsContainer
    gc = m_pMxDoc.ActiveView
    gc.Reset()
    Dim mapSurround As IMapSurround
    mapSurround = pArrow  ' QI
    pArrow.Size = 125

    Dim mapSurroundFrame As IMapSurroundFrame
    Dim mapFrame As IMapFrame = gc.FindFrame(m_pMxDoc.Maps.Item(0))
    Dim pUid As ESRI.ArcGIS.esriSystem.UID = New ESRI.ArcGIS.esriSystem.UID
    pUid.Value = "esricore.MarkerNorthArrow"
    mapSurroundFrame = mapFrame.CreateSurroundFrame(pUid, pArrow)

    '***** set anchor point
    Dim element2 As IElementProperties3 = mapSurroundFrame
    element2.AnchorPoint = esriAnchorPointEnum.esriTopLeftCorner

    '***** white background
    Dim backGround As IBackground
    backGround = mapFrame.Background ' match background to map frame
    backGround.Gap = 3
    Dim symBackground As ISymbolBackground = backGround
    Dim fillSymbol As IFillSymbol = symBackground.FillSymbol
    Dim colour As IRgbColor = New RgbColor
    colour.Red = 255
    colour.Blue = 255
    colour.Green = 255
    fillSymbol.Color = colour
    symBackground.FillSymbol = fillSymbol
    symBackground.Gap = 2
    mapFrame.Background = backGround

    activeView = m_pMxDoc.ActiveView
    Dim pElem As IElement = mapSurroundFrame
    pElem.Geometry = NorthGeometry(x, y)
    gc.AddElement(pElem, 0)
End Sub


texts: Is there a way to fix the polygon-pattern scale to the zoom level in CartoCSS

texts:
I'm using openstreetmap-carto for my own tileserver and added a new landuse to the style. But I can't wrap my head around, how I could fix the scaling of my polygon-pattern to a specific zoom level.
For simplicity, let's say the pattern is a grid. With zoom level 16, the grid is perfectly 100+100 meters. What I want to achieve is that the pattern scales proportionally to the zoom level, so that the grid is still 100x100 meters in other zoom levels.
Do I need to use the polygon-pattern-transform: scale(x) parameter or is there an easier way?


texts: MapServer - Not Rendering All Features With Complex DATA statement

texts:
I have run into a little problem when using MapServer to render complex data statements.  To keep this post short and simple, I think the output of shp2img is the best place to start.  When I use shp2img with a simple DATA statement (ex: GEOMETRY FROM FOO_TABLE USING SRID 4326) AND the -all_debug 5 flag set - an msOracleSpatialLayerNextShape line is written to the console for every single row.

When I use a more complex DATA statement (see below) - shp2img writes out about a dozen of the msOracleSpatialLayerNextShape line - AND THEN the last line is "Fetched 1234 more rows"....  It does not loop through every row like above.  When I run the sql that shp2img uses to retrieve the data - it is right on the money.  Both shp2img and SQLDeveloper report the same # results.

Is there some reason why shp2img doesn't loop all the results?

Also, this only seems to be an issue with LINE geometry type.  The same type of data statements work as expected with POINT.

Complex Data Statement:

DATA "GEOMETRY FROM (SELECT CASE WHEN F.ch_v_avg_a_min IS NULL THEN 0 WHEN F.ch_v_avg_a_min > 242 THEN 1 WHEN F.ch_v_avg_a_min < 242 THEN 2 ELSE 0 END as THRESHOLDCLASS, N.SHAPE AS GEOMETRY FROM CHN_DIS_NODE_MONTH_FACT F INNER JOIN GIS_NETWORK_FEATURE N ON N.OBJECTGUID = F.NODE_FEATUREID WHERE N.NETWORKELEMENTID = 2 and F.MONTHID = TO_DATE('2014-05-01', 'YYYY-MM-DD') AND N.SHAPE IS NOT NULL) USING NONE"


texts:
I would like to know how and where I can download qgis2web Plugin.
I use QGIS with 2.10.1 version.


texts: Offset line, when zooming out line extends at intersections

texts:
[1[]2I'm doing offset in ArcGIS in a cartographic line symbol. Before the offset the line in all scales appears OK. But when zooming out, the offset line extends out at intersections. While when zooming in, the offset line appears OK.

Could some one help me on that issue?


texts: How to create a table from existing tables (applying new schemas) in postgis DB?

texts:
I have a table that is geospatially enabled; that contain world borders. I want to break it down to different schemas (global, europe etc..) . My thought is to use use something that looks like this:

CREATE TABLE europe.borders
  AS (SELECT *
         FROM wo_borders
         WHERE admin_lvl2='eu' ); //just an example!

What's the correct to create a tables from existing tables taking account new schemas and making them geospatially enabled? Do I have to create beforehand the new tables using my default postgis template?


texts:
Using QGIS,I attempted  Africa's watershed delineation taking the following steps:

Downloaded DEM  from earthexplorer USGS  (GTOPO30)
Created a mosaic of the downloaded DEM
Reprojected the mosaic DEM to EPSG:54010 World_Eckert_ VI(For equal area world map projections)
Rescaled map cell size to 1km by 1km(For faster processing time) using the Align Rasters function
Filled depressions using   Wang & Liu XXL algorithm
-Created stream order using strahler order function
Choose  strahler order threshold using  raster calculator function
Created the basins and channels using SAGA channel network and drainage basins function

The output (Picture 1) do not look satisfactory.

How do I improve the results?

The goal is to have a clearly delineated River Niger Basin- the basin with thickened blue border line in picture 2


texts: PyQGIS Creating Numpy Array from QByteArray (QgsRasterBlock.data()) and vice versa (QGIS3)

texts:
I am looking for a way to convert the values contained in a QgsRasterLayer object to a numpy array and vice versa without using any GDAL methods.
To obtain the numpy array I go through QgsRasterDataProvider to QgsRasterBlock and then read the values into a numpy array through the QgsRasterBlock.value(i,j) method. This is not very elegant, but runs in an acceptable time.
Analogously - after working with the numpy array - I change the values in the QgsRasterBlock through the QgsRasterBlock.setValue(i, j, value) method, which takes forever and is not acceptable for my application... so I am looking for an alternative!

I've read that handling the pixel values in the QgsRasterBlock object as QByteArrays and creating numpy arrays directly from the QByteArrays is possible and efficient (as compared to the method described above).

So my first question: Is this true?

My second question: What does the numpy function to create an array from QByteArray look like? And which function can I use to go back again?


texts: How to edit geometry of multiple layers simultaneously in QGIS 3?

texts:
I have several line layers (route, pipes, cables) that share the same geometry and would like to edit them simultaneously.

Is there a way to easily snap the geometry of the pipes and cables layers to the route layer so that I can simply change the geometry of the route layer and don't have to repeat all changes on the other two layers as well?

To clarify, I plan a fibreglass route and have three layers. The first one describes the route that I have to take to pass all homes that have to be connected. The second one contains the pipes in which the cables lie and the third layer is the cable layer. So all features on those layers have to have the same geometry. Once I change the pathway of the route layer, I have to apply the same changes to the other two layers to ensure that they are still aligned with the route layer.
I'm not sure whether I understand your comments correctly, but I use the basic editing tools in QGIS 3.


texts: How to open and save QGIS datasource with/to lower case fields?

texts:
I want to be able to change all field names to lower case from a data source when loading them in QGIS (SHP in particular).

Is it possible? Via the GUI, via Python.
Is it possible to convert field names to lower case when saving a layer with the "Save as" dialog?

Why? I regularly go back and forth between GIS files from OpenData portals and databases like PostGIS or Spatialite. Doing queries with upper case fields require quoting in these database and I want ease the pain. I would like to convert all my SHP field names to lower case.

For info, I've seen it's possible when using Grass e.g the -w option in v.in.ogr.
I've seen using bash with ogrinfo is also possible but remains less friendly IMO.


texts: Google Earth Engine: Maintain area calculation accuracy while reducing resolution and exporting from

texts:
I am struggling to find a best practice for exporting data from Google Earth Engine at a reduced resolution for further analysis outside of GEE, while maintaining similar results  analysis in the code editor.  My principle questions are when and how to define projection to maintain positional and area accuracy.

You have to project in order to reduce resolution to a specified scale, but do you also include projection/scale information on export?  Do you reproject the raw data first, before reducing resolution, then project again to define a larger pixel size?  Do you work with a snap raster before exporting?  I have tried various ways of doing this and am so far unable to maintain similar summary statistics at a reduced resolution.  I know reprojecting and adjusting resolution changes pixel values, but I think I should at least be able to get close. I would be open to javascript of python solutions, as I'm thinking of migrating workflow to the python API.

My workflow currently:

var reduced.image = RawImage.reduceResolution({
                 reducer: ee.Reducer.mean()
                 maxPixels: maxPixels //number of pixels to aggregate
                 }).reproject(crs, null, scale)//crs taken from desired NLCD projection//scale pixel size of intended image

//exporting
Export.image.toDrive({
  image: reduced.image,
  description: 'reduced.image',
  folder: 'reduced.image_folder',
  scale: scale,
  region: region,
  maxPixels:1e13,
  //crs: crs  //do I include crs and transform again here?
  //crsTransform: crsTransform,
  //pyramidingPolicy:{".default":"mode"},
});


texts: How to drop a table in spatialite that starts with a certain string?

texts:
I am using spatialite (v4.2.0) on a Windows 7 (64bit) PC.  I need a way to drop a specific table that always begins with "_old" from my database.  The file may not be in the database in the first place.  I do not know the rest of the table name since it can change (ex:  _old_20140801, or _old_20140715, etc.).  I can always find out if a table that fits my criteria exists with the command:

.tables '_old%'

but how do I save that result (is it storable in a local variable?) and then use it?


texts:
I have installed GeoServer 2.5.5 on Win 7. When I try to view a raster pyramid I see this message:

OutOfMemoryError : Java heap space.

I google it, but it wasn't useful because I don't use Tomcat and it was no wrapper.conf in my PC.
How can I increase Java heap space?


texts:
I try to merge some multiband rasterfiles (Sentinel 2 data) with Rasterio into a single rasterfile but the result is incorrect.
The files are too big to fit in my RAM completely, but I have enough SWAP and I don't get a Error Message.
Any ideas, what could have gone wrong?

My code:

import rasterio
from rasterio.merge import merge
import glob

files = glob.glob("../data_multiband/*.tif")
mosaic = []
for file in files:
    tmp = rasterio.open(file)
    mosaic.append(tmp)
mosaic, out_trans = merge(mosaic)
out_meta = tmp.meta.copy()
out_meta.update({"driver": "GTiff",
                 "height": mosaic.shape[1],
                 "width": mosaic.shape[2],
                 "transform": out_trans})
with rasterio.open("../data_merged/test.tif", "w", **out_meta) as dest:
    dest.write(mosaic)

That is how the result looks:

As you can see (red marks) some rasterfiles occur two times. The one on the left side are wrong.
Thats how the original files look:


texts: QGIS 2.14 Ubuntu 16.04 "Project properties" menu missing

texts:
Could you help me to understand why my QGIS configuration doesn't show me the Project properties sub menu found in Project menu ?

QGIS version : 2.14

Ubuntu distrib : 16.04


texts: How to label the start\ end of the line in QGIS?

texts:
I work with this line layer:

The label is placed in the middle and i seek for a way to place it in the start\ end of the line. I search for that option in the placement tab in the properties but with no results.

I know that when i pan the map the label jump for the best location, but still i need to place the label in the start\ end of the line.

I try to find an option similar to the simple function in the ArcGIS:


texts: Iterate through the attribute tables of several SPDFs in R

texts:
In short I want to know:
What would be the appropriate way to iterate through several SPDFs and do a calculation on every row in their attribute tables?

Some more details: I have some SPDFs and I can iterate through the polygons of each one using something like:

x = readOGR(...)
for (i in 1:length(x)){
    x$Volume_Norm[i] <- x$Volume[i]/max(pre_max)
    x$Volume_Sqrt[i] <- sqrt(x$Volume_Norm[i])
    x$Volume_Thrd[i] <- nthroot(x$Volume_Norm[i], 3)
    x$HAZARD_1[i] <- x$Volume_Thrd[i] * x$HAZARD_NEW[i]
  }

But I don't know how I could work this into a apply-function or something similar, so that it works on all my SPDFs and I don't have to change the names in the for-loop separately.
If I try:

list <- list(x, y, z) # these are my SPDFs
fun <- function(x){
  for (i in 1:length(x)){
    x$Volume_Norm[i] <- x$Volume[i]/max(pre_max)
    x$Volume_Sqrt[i] <- sqrt(x$Volume_Norm[i])
    x$Volume_Thrd[i] <- nthroot(x$Volume_Norm[i], 3)
    x$HAZARD_1[i] <- x$Volume_Thrd[i] * x$HAZARD_NEW[i]
  }
}

lapply(list, fun)

nothing changes and I get :

> lapply(list, fun)
[[1]]
NULL

[[2]]
NULL

[[3]]
NULL

I've asked the same question over at https://stackoverflow.com/q/61549967/11728010, but got no reply so far, so I thought, I could try it here.


texts: How to Zoom and Pan using scroll wheel in ArcObjects?

texts:
is there any way to zoom using scroll wheel ? i have tried AutoMouseWheel of map control, it works well for the first used custom tool, but once i change the tool, it stops working , i would like to know the reason of this weird behavior.

i don't know if it is possible to catch mouse wheel events in custom tool ? that's seems to be an alternative solution.

Any suggestion or help is welcome :)


texts: Can QGIS ask for the preferred user profile on startup?

texts:
I'm using QGIS user profiles for different tasks and would like to choose the right one  at QGIS startup - is that possible?
I know how to set the default profile by editing the profiles.ini in QGIS profile folder, but I would like to choose it on QGIS start.
I also found that solution by starting qgis in a specific profile via command line: Launching QGIS with specific user profile


texts: Generate classification schema in ArcGIS Pro using labeled raster - fail

texts:
Trying to do supervised classification of L1TP Landsat 8 OLI bands using USDA Crop Data Layer .tif file through ArcGIS Pro. The Crop Data Layer file has one band, with labels like "Corn", "Cotton", "Rice." Doing image classification through ArcGIS pro requires a classification schema, with the default being National Land Cover Database. Classification Wizard offers "Generate Classification Schema from Classified Raster," but inputting the Crop Data Layer .tif file into this generates the following error: "Failed to generate classification schema from classified raster." No other explanations attached.

Not sure if there is something inherently off about the CDL file format/way the info is organized? All search results about generating classification schema lead to drawing polygons--should I export the CDL raster as a polygon?
Below is the attribute table for Crop Data Layer:


texts: How do you disable archiving for a table?

texts:
I would expect to have a Disable Archiving... like I do for a featureclass.

But for a table, the Enabled Archiving menu-item is grayed out.

Here's the menu for a table:

And here's the menu for a featureclass, which is what I would expect according to the documentation:

Is there another way to disable archiving on a table? Or is this a bug?

ArcCatalog 10 SP3 client.  Seeing the issue in both SQL Server/Oracle SDE.


texts: How to reproject manually created label positions?

texts:
I changed position to some features' labels and I activated the x,y position in the label properties.
I saved the layer in .shp and with another projection.
Once created a new project with that projection I loaded the .shp.
Now I can see only labels of those features I didn't move the label. So I think that labels' X,Y coordinates saved in the layer aren't reprojected

How could I have these coordinates reprojected?


texts: QGIS: How to modify the appearence of selected attributes

texts:
I am using QGIS version 2.14.1 (Essen). When I open the attributes table of a layer, I select some data which is then highlighted in yellow on the map. I was wondering if it is possible to change the appearance of this yellow (play with the transparency or colour)?

The other idea is to create a single shapefile from all the different attributes but that would take too long as I have lots of data.

Anyone know if this is possible in QGIS?


texts: Lat/Long data to shapefile with many values in CSV dataset

texts:
I have my polygon spatial data in a format of CSV with the column of lat/long delimited by comma. Is there any way to upload this csv and got polygon shapefiles (or any other visualisable format)? My primary GIS tool is QGIS.
There are several ways how to do it just for one WKT at the same time but I cannot find a way how to do it with data with several rows in data.
I have tried to Add delimited text layer with geometry definition of WKT but it does not work (as can be seen in the screenshot. Data are not in the readable WKT format.
I am up to some solution even outside QGIS environment


texts:
Troubleshooting and debugging ArcMap defintion queries that go beyond a basic single statement can be very frustrating. Statements that work on their own easily break when combined. How can this pain be alleviated?

For example say we want to select all records of last 3 years, and then from that set only the month of march.

In a file-gdb this query selects all records 3 years old and newer:

EXTRACT(YEAR FROM "Date") >= (EXTRACT(YEAR FROM CURRENT_DATE) - 2)

While this query selects the month of March:

EXTRACT(MONTH FROM "Date") = 3

However combining them by inserting year_query in place of 'Date' in the month query will fail with "an invalid SQL statement was used". (On multiple lines for readability, it's all in a single line in query builder):

EXTRACT(MONTH FROM
   (EXTRACT(YEAR FROM "Date") >= (EXTRACT(YEAR FROM CURRENT_DATE) - 2))
) = 3

For the example the solution is to use AND instead of nesting:

(EXTRACT(YEAR FROM "Date") >= (EXTRACT(YEAR FROM CURRENT_DATE) - 2))
AND (EXTRACT(MONTH FROM Date ) = 3)

but how to get this answer when you don't know it already?


texts: How to change vector tile source provider in TangramJS

texts:
I rely on the Mapzen API to serve vector tiles to TangramJS which render in Leaflet using a style similar to Bubble Wrap (https://github.com/tangrams/bubble-wrap/blob/gh-pages/bubble-wrap-style.yaml).

Since the unfortunate announcement of Mapzen closing it's services in Feb 2018 I am looking to continue using TangramJS and hopefully just a difference tile source.

Do the other providers (Mapbox, OpenMapTiles, ThunderForest) suggested by Mapzen's blog post (https://medium.com/postzen/mapzen-alternatives-aa6d857e3532) just 'drop in' their own URL and things continue on exactly the same?

Below is the sources section from the Bubble Wrap style.

sources:
    mapzen:
        # EraserMaps the app only supports MVT format tiles
        type:  MVT #TopoJSON
        url:  https://tile.mapzen.com/mapzen/vector/v1/512/all/{z}/{x}/{y}.mvt
        #url:  http://tile.dev.mapzen.com/mapzen/vector/v1/all/{z}/{x}/{y}.mvt
        #url:  //localhost:8080/vector/all/{z}/{x}/{y}.mvt
        # road labels in Tangram JS are broken when overzooming, set max_zoom: 18 to preview fix
        url_params:
            api_key: global.sdk_mapzen_api_key
        tile_size: 512
        max_zoom: 16


texts: How to add data from jTable to JSON file

texts:
I want to create a tools that add my selected item to the JSON file in java netbeanse;

for this I make a ComboBox for selecting my item , and also I have a jTable for display the ArrayList of selected item, and the items in jTable is add/remove by clicking on [add/delete] Button from combobox list item.

So my question is how to add my selected items which is display in jTable into json file. and if I delete item form the jTable it is also delete on JSON file.
finally I need the json file with my selected items,
you can see the schema of my json file in below:

Do you have any idea about this.

Adding:

private void jButton1ActionPerformed(java.awt.event.ActionEvent evt) {

DefaultTableModel model = (DefaultTableModel) jTable1.getModel();
int row1 = jTable1.getModel().getRowCount();
model.addRow(new Object[]{row1+1,jComboBox1.getSelectedItem()});
jTable1.setModel(model);

Deleting:

private void jButton2ActionPerformed(java.awt.event.ActionEvent evt) {

 int row = jTable1.getSelectedRow();
 if(jTable1.getRowCount()>0){
 DefaultTableModel model = (DefaultTableModel) jTable1.getModel();
 model.removeRow(row);
};

JSON file:

{
    "tags": {
        "load_all": true,
        "exclude": [
            "created_by",
            "source"
        ]
    },
    "use_single_id_space": true,
    "tables": {
        "all": {
            "fields": [
                {
                    "type": "id",
                    "name": "osm_id",
                    "key": null
                },
                {
                    "type": "geometry",
                    "name": "geometry",
                    "key": null
                },
                {
                    "type": "hstore_tags",
                    "name": "tags",
                    "key": null
                }
            ],
            "type": "geometry",
            "type_mappings": {
                "points": {
                     "Item 1": ["__any__"],
                     "Item 2": ["__any__"]

                },
                "linestrings": {
                     "Item 1": ["__any__"],
                     "Item 2": ["__any__"]
                },
                "polygons": {
                     "Item 1": ["__any__"],
                     "Item 2": ["__any__"]
                }
            }
        }
    }
}


texts: Values extracted from NetCDF not matching original raw data

texts:
I am working with monthly rainfall data from University of Delaware which is available as a NetCDF from NOAA and raw text from UDel's homepage. I would prefer to use the NetCDF format but when I extract values, these don't match what the ones in the raw data files.

Can you suggest why the comparison might not be working?

# Extract data
from netCDF4 import Dataset
ncfile =  "C:/precip.mon.total.v401.nc"
fnc = Dataset(ncfile, 'r')
print fnc.variables
lat = fnc.variables['lat']
lon = fnc.variables['lon']
monthyear = fnc.variables['time']
rain = fnc.variables['precip']
# Manually check values / example one
print lat[37] # 71.250
print lon[0]-180 # -179.750 (need to subtract 180 because the lon values in the nc file range from 0 to 360)
print rain[0:12, 37, 0]*10
# Wrong output: [-- -- -- -- -- -- -- -- -- -- -- --]
# Another example
print lat[46] # 66.750
print lon[54]-180 # -152.750 (need to subtract 180 because the lon values in the nc file range from 0 to 360)
print rain[0:12, 46, 54]*10
# Wrong output: [26.6   30.5   23.6   27.1   29.7   26.1   101.9   44.4   31.   79.8   49.9   29.7]

Acc to the raw data the values corresponding to lon -179.750 and lat 71.250 are:

0.0     5.4     0.0     4.0     7.3     3.9    33.5    63.2     9.4     4.0     8.0    23.6

In case of lon -152.75 and lat 66.75, the extracted data should have been:

11.6    11.4    26.4    28.5    14.3    31.7    30.5    67.0    59.6    22.9     9.1    33.3

(The netCDF may be downloaded from https://www.dropbox.com/s/6zi0imr5i1qlqgs/precip.mon.total.v401.nc?dl=0 and the raw data file is available here https://www.dropbox.com/s/6klwmj37kvwj46w/precip.1900?dl=0)


texts: Time of the gps fix in GPGGA sentence is always rounded

texts:
I am using QSTARZ 818XT device for my app.
I am using GPGGA sentences to get the GPS fixes, but I have noticed that timings are always rounded like xxxxxx.000 - for 1Hz GPS frequency and xxxxxx.x00 for 10 Hz GPS frequency.

Why is this happening?
Is that something device specific?
Can I trust these timings or should I make some corrections?
My app does not function in real time, but it really needs those milliseconds.


texts: How can I convert a MapDB file to other formats?

texts:
I'm currently running a webapp locally that stores/caches data in MapDB format.
I would like to play with the data in PostgreSQL. Do there exist tools to convert the .db files to other formats?


texts: How to visualize traffic data in a Mapbox map?

texts:
I am trying to visualize traffic data which I have estimated using FCD data. The traffic data format is as follow :
osm_start_node | osm_end_node | speed(Kmh)
  113054533    |  113096757   |   54
  170190194    | 4525170049   |   93
  172637810    |  172637809   |   70

I want to know how should I format this data to pass it to the Mapbox traffic layer and visualize it on the same Mapbox map using MapBox GL-JS and to have a result similar to this :


texts: Plotting points conditionally using cluster number in qgis

texts: I have a set of locations such that there are partitioned into five clusters, I used the categorized property under styles to differentiate points belonging to different clusters using different colors but I want to run queries such that I can select and display the points on the map according to cluster numbers using the queries in qgis.


texts:
I have been able to batch clip rasters in one directory using python in ArcGIS.

However, I am not sure how to do the same in QGIS using the Python console.

Is this possible?


texts: Select OSM admin boundary which links with pelias geocoder result

texts:
We use pelias geocoder environment in our project. There is a case "selection a polygon, which links with pelias result." For example we choose "Muenich" in search field -> on client we see this:

How to realise it? Geocoder answers have only bbox and node coordinate (if i use &source=osm too), how to links them with osm admin boundary?


texts: Is there software that will do reverse geocoding to convert photographs' GPS locations to addresses?

texts:
I am using a standalone GPS tracker that produces GPX files and then clock-sync photos with it using Microsoft Pro Photo Tools. This works pretty well. But I was also using a feature of Pro Photo Tools called "Get location text" which uses the latitude/longitude information from EXIF and translates it using Windows Live services to actuall addresses. At least it used to, but stopped recently, probalby because of some change of the online API or whatever, and a lack of update for the software.

So, the questions is: do you know of any other software that can bulk translate location info into addresses and write that back into the photo's metadata?

I expect the solution to use Google Maps or Bing Maps APIs - both provide (reverse) geocoding world-wide and publicly, and work quite well through their websites. I find it hard to believe that no one has created an app yet that would do this in bulk for photos.


texts: mod_tile renders OpenStreetMap without some water

texts:
My rendered map:

The same view from openstreetmap.org:

What went wrong?

I downloaded europe_latest.osm.pbf and nigeria_latest.osm.pbf and imported these files using --append. Now, I'm trying to import it again by merging both .pbf files - will see how it goes.


texts:
I'm referring to a question I found here (How can I move diagrams?) about pie charts in QGIS. I too would like to move some pie charts because they are overlapping each other. I tried to follow the instructions on the link, but didn't understand the approach properly.
I was able to create a dbf with the centroids of the polygons, however, i am not sure if the centroids were of the polygons of the layer OR the pie charts. Also, I did not know how to move the diagrams with the dfb-file I had created.

I'm using QGIS 2.18.3. (Las Palmas).


texts:
I was installing TileStache and and followed all the necessary steps, Upon running the tilestache-server.py file, It starts the server at http://127.0.0.1:8080/. And opening this into my browser gives the message 'TileStache bellows hello' which is okay. But when I try to run the url for osm basemap layer which is http://localhost:8080/osm/0/0/0.png it always gives me an error message saying "Internal Server Error".

I have tried installing different TileStache versions but the issue is same with every version I install. I am not sure about the reasons that are causing this issue.

Screenshot of the error shared below:


texts: From a Grid map to a dataframe of presence/absence of each cell

texts:
Working in R, I have a grid map of the world (as SpatialPolygonsDataFrame object). I have also many indipendent shapefiles (formed by several polygons) expressing species ranges. I grouped all these shapefiles representing distributions in a single list.

species <- list(sp1, sp2, sp3, sp4, sp5)

What I want to obtain as output is a dataframe where for each cell of my grid (put in the rows of the dataframe) I have in the columns the presence/absence (expressed by the values 1 and 0, respectively) of all species of my list. See the example I want to obtain:


texts: How much sense does it make to average (lat,lon) samples in order to increase 2D accuracy of a GPS location?

texts:
Several GPS applications, like this one, or this one, acquire multiple (lat,lon) samples of a given location, assuming that the GPS unit is not moving, and then take an average of the samples in order to compute a "more precise" 2D location.

(We do not care about the elevation/altitude position here!)

The second app (GPS Averaging) uses the accuracy value associated to each sample as a weight for the current location, and then calculates the weighted average accordingly. It also provides an estimate of the accuracy of the averaged location.

Questions:

1) While common sense pushes us to believe that averaging should lead to increased accuracy, how much sense does it make for handheld units like phones (i.e. simple devices that do not use differential GPS)?

2) Would you recommend another method than GPS Averaging's method to compute the average location?

3) How to compute an estimate of the accuracy of the averaged location?

4) Is there a way different than averaging to get a better 2D positioning by acquiring multiple (lat,lon) samples of a given location?

UPDATE 1: results of my preliminary study with 2 handheld GPS units (Sony phone models ST15i and ST17i) acquiring 3m accuracy fixes at the same position during 4.5 hours gave the following data:

=> It is quite interesting to note that even though the supposed accuracy of the fixes was 3 meters, the ST17i model had lots of points further than 3 meters from the median/average.

=> Also remarkable is the monotone drift of the longitude on the ST15i model.

(Note that the ST15i seems to have a more sensitive antenna than the ST17i as I could analyse it used on average 3 more satellites for its fixes than the ST15i!)

UPDATE 2: some further stats and numbers, still from the same datasets

=> The data is definitely not normal

=> I also computed the distance between the median location of the ST15i and the median location of the ST17i: it is 3 meters, as if the study was playing with us, since all fixes used had an accuracy of 3 meters or better. This definitely validates the suggestion below of using a known reference in order to derive meaningful conclusions about the accuracy of each GPS unit!


texts: No clue of the Georeferencing tool in gvSIG

texts:
I have kept trying to locate the georeferencer in gvSIG but I could not find it. Some tutorials say it is under Tools from the main menu!
I installed several versions of gvSIG, both on Windows and Linux but the georeferencer does not appear anywhere. I also tried to install a georeferencing plugin from http://downloads.gvsig.org/download/gvsig-desktop/pool/org.gvsig.georeferencing/ but after I installed it, the gvSIG won't start.
Any idea?


texts:
This has kind of been asked before in different ways, but I'm struggling to adapt the answers to my specific scenario.

Say I have a ground station and a plane in the sky. I know the GPS lat, lon, and altitude of both. I'm trying to derive the XYZ position of the plane relative to the ground station on a new grid system where the ground station sits at (0,0,0), the Z-axis points out of the ground station into the sky, the Y and X axis are tangent to the altitude of the ground station, and the Y axis projects over the north pole (points north when zoomed in close).

Here's a picture to try and convey this:

The orange dot is the ground station, the brown dot is the plane, and the red axes are how I'd like to define my new relative grid.

I tried doing this: Convert the lat/lon/alt to earth XYZ coordinates, then rotating about the earth Z axis by ground-station longitude, then rotating about the earth Y axis by ground-station latitude, then translating along earth X axis for ground-station altitude (plug earth radius). I then relabeled the axis: earth-X = Z2, earth-Y = X2, and earth-Z = Y2. I got a result that placed the ground station at (0,0,0), but the plane results don't look right: the altitudes shown on the Z-axis are out of order compared to the plane's GPS altitude readings. I'm testing with plane coordinates that are close enough to the ground station that the curvature of the earth should not affect the Z-axis very much (maybe a couple hundred feet, while these planes are at 5000ft to 48000 ft).


texts: Using Network Links updates successfully with ArcGIS Explorer?

texts:
I've been trying to do updates of kml through a network link without success.

Here is an example from the KML 2.1 tutorial.

Load the following file (Point-load.kml) into AGX:

<?xml version="1.0" encoding="UTF-8"?>
<kml xmlns="http://earth.google.com/kml/2.1">
<NetworkLink>
  <name>Loads Point.kml</name>
  <Link>
    <href>http://code.google.com/apis/kml/documentation/Point.kml</href>
  </Link>
</NetworkLink>
</kml>

This should create two placemarks.

Then load the following file (Update-load.kml):

<?xml version="1.0" encoding="UTF-8"?>
<kml xmlns="http://earth.google.com/kml/2.1">
<NetworkLink>
  <name>Update</name>
  <Link>
    <href>http://code.google.com/apis/kml/documentation/NetworkLinkControl-Update.kml</href>
  </Link>
</NetworkLink>
</kml>

This should update one of the Placemarks to a different name and does in Google Earth, but doesn't seem to work in ArcGIS Explorer.


texts: Segmentation fault (core dumped) QgsApplication

texts:
import os, os.path,sys
from    qgis.core   import  *
from    qgis.gui    import  *
from    PyQt4.QtGui import  *
from    PyQt4.QtCore    import  *

class   MapExplorer(QMainWindow):
    def __init__(self):

        QMainWindow.__init__(self)
        self.setWindowTitle("Landmark Explorer")
        self.resize(800,    400)

def main():

    QgsApplication.setPrefixPath(os.environ['QGIS_PREFIX'], True)
    QgsApplication.initQgis()
    app =   QApplication(sys.argv)
    window  =   MapExplorer()
    window.show()
    window.raise_()
    app.exec_()
    app.deleteLater()
    QgsApplication.exitQgis()

if  __name__    =="__main__":
    main()

Then run this code and get:

Segmentation fault (core dumped)

What did I do wrong?

OpenSuse 42.2


texts: How to make QGIS 3.18 display Chinese characters in labels

texts: I am using vector tiles from Maptiler in QGIS 3.18. I have set the language of the tiles to Chinese on the maptiler website. The labels display in QGIS in Pinyin (romanization of Chinese characters.) I need to see the actual characters. Actually the ideal situation would be English and Chinese side by side (this is the expected behavior based on the maptiler preview map) but I would settle for just Chinese characters.


texts: Display points with lowest elevation within a 1m2 grid - QGIS 2.18.2

texts:
I have a few thousand points measured with totalstation from an archaeological site. The points are measured quite unevenly which makes TIN triangulation and subsequent DEM generation a mess. Some areas has a lot of tightly measured points, others almost none.

I have created a 1m2 grid overlaying the points containing XYZ.

What I would like to achieve: Display one point for each 1m2 and that this one is the one with the lowest Z elevation value.

PS: Yes, the Z-values are negative, because the excavation was conducted on old seabed.


texts: Add and fill a field in Python Toolbox (pyt) in ArcMap

texts:
I am using python toolbox in ArcMap, and I want to add a new field in a shapefile and after that I want to fill this field with the name of the shapefile.

I achieve to add a new field using arcpy.AddField_management function, but then I cannot fill this field. I am using arcpy.CalculateField_management function to fill it, however I think that this is a wrong way.

Any suggest is welcome.

The code is:

import arcpy, os, sys
from arcpy import env

class Toolbox(object):
    def __init__(self):
        """Define the toolbox (the name of the toolbox is the name of the
        .pyt file)."""
        self.label = "Toolbox"
        self.alias = ""

        # List of tool classes associated with this toolbox
        self.tools = [Tool]
class Tool(object):
    def __init__(self):
        """Define the tool (tool name is the name of the class)."""
        self.label = "FIELD"
        self.description = ""
        self.canRunInBackground = False

def getParameterInfo(self):
    """Define parameter definitions"""

# First parameter

    param0 = arcpy.Parameter(
        displayName="Shape entrada",
        name="in_features",
        datatype="DEShapefile",
            parameterType="Required",
        direction="Input")

# Second parameter

param1 = arcpy.Parameter(
        displayName="Nuevo campo",
            name="field",
        datatype="Field",
            parameterType="Required",
        direction="Input")

    parameters = [param0, param1]
    return parameters

def isLicensed(self):
    """Set whether tool is licensed to execute."""
    return True

def updateParameters(self, parameters):return

def updateMessages(self, parameters):return

def execute(self, parameters, messages):
    """The source code of the tool."""
    a= parameters[0].valueAsText
b= parameters[1].valueAsText

a=arcpy.AddField_management(a, b, "TEXT", "", "", "", "", "NULLABLE", "NON_REQUIRED", "")

a=arcpy.CalculateField_management(a, b, "hola" ,"VB")

    return a


texts: Performing raster cleanup in ArcObjects?

texts:
I am working on Raster cleanup tool and using below but is the output rectangle or block?

http://edndoc.esri.com/arcobjects/9.2/ComponentHelp/esriArcScan/IConnectedCells_LocateByDiagonalOfEnvelope.htm

Private Sub LocateByDiag()

    Dim pMxDoc As IMxDocument
    Dim pConCell As IConnectedCells
    Dim pFoundBlock As INotifySelectionBlock

    Set pConCell = Application.FindExtensionByName("ESRI ArcScan Tools")
    Set pMxDoc = ThisDocument
    Set pFoundBlock = New FoundBlock

    pConCell.LocateByDiagonalOfEnvelope 100, True, True, True, pMxDoc.ActiveView.Extent, pFoundBlock, Nothing, Nothing

End Sub

Place the following code into a class and name it FoundBlock
Implements INotifySelectionBlock

Private Sub INotifySelectionBlock_SelectionBlock(rectangle As tagRECT, ByVal block As Variant, ByVal cancel As ITrackCancel)
    Debug.Print "Found block"
End Sub


texts: Google Earth Engine: pop-up when clicking on geometry/Feature

texts:
This may be more of a JavaScript question than a GIS question.
I'm adding feature points to a GEE map (https://code.earthengine.google.com/2d6c3146e5468cb9f14e10039ca22b59). My goal is to map point features from an asset table (which has associated data), and view those data when clicking on the features on the map.

Google Map API has an infowindow function (https://developers.google.com/maps/documentation/javascript/examples/infowindow-simple) but I'm not sure if or how it can be incorporated into a GEE script.

Is it possible to have a pop-up bubble (info window) show up when a feature is clicked? I have looked through GEE guides but couldn't find any function already existing.

Here is an example script:

// Import DEM
var collection = 'WWF/HydroSHEDS/03VFDEM';
var band = 'b1';
var dataset = ee.Image(collection).select(band);
var elevationVis = {
  min: -50.0,
  max: 3000.0,
  gamma: 2.0,
};
// create point geometry (will call an asset in the real script)
var point = ee.Geometry.Point([-61.334, 15.416]);

Map.setCenter(-61.334, 15.416, 10);
Map.addLayer(dataset, elevationVis, collection);
Map.addLayer(point, {}, "acou");

// example from Google Map API
// should be placed earlier in the script but causes GEE to stop runing
var contentString = '
'+
            '
Hello'+
            '
';
var infowindow = new google.maps.InfoWindow({
          content: contentString
        });
var marker = new google.maps.Marker({
          position: point,
          map: map,
          title: 'Uluru (Ayers Rock)' });
        marker.addListener('click', function() {
          infowindow.open(map, marker);
        });


texts: Sub-dividing a polygon into smaller polygons

texts: In OpenStreetMap data, some large areas are tagged with a residential=yes tag. Actually within this large polygon there are multiple smaller (rectangular) houses (which are actually not drawn as yet). I would like to 'slice' the large polygon into smaller rectangular polygons. Or looking at it another way, I would like to or auto-generate as many fixed size rectangles (houses) into the larger area. I would like to achieve this either using QGIS or MapWindow, as I do not have ArcGIS. Is there any way to do this?


texts: Change rendering order of objects in MapInfo table?

texts:
i have one big MapInfo file where i did a lot of digitizing work. Now i have the problem of disapearing points under polygons (or polygons under polygons, lines under polygons etc.).

How can i manipulate the order objects are rendered?

My solution so far is to create a couple of more physical layers and add objects to them when i want them "above" the main layer. I have to send the digitized map to a customer, i don't want to send her 4 files just to have the outlines of a park bench not dissapearing under the green "grass" polygon it is situated in.


texts: Creating filter checkboxes for Google Maps Javascript Web Map?

texts:
I'm currently using google fusion tables to store research data and visualize it in the Google Maps Javascript API. I'm new to coding and am having trouble developing some functional components of the web map. I was hoping someone here could help me out with a filtering by checkboxes component.

I am currently working off the IN query Google Fusion Tables Layer example found here: https://developers.google.com/fusiontables/docs/samples/in

Basically what I'm trying to develop is a series of checkboxes so that users can filter the dataset by year. Below is the code I have so far, some of which I tried to steal and alter from the example. I think there is an error in my Javascript logic but like I said, I'm a novice at best.

Below is my code:

https://gist.github.com/cv26/1452706d9513215cb5e90ef9b13d070d


texts: Inserting hyperlink to PDFs for each Polygon using QGIS?

texts: I am currently using QGIS 2.18.16 and would like to know if you can you add a hyperlink to each line in the attribute table, I am wanting to hyperlink to PDF reports for each polygon.


texts: ArcGIS: Layer to KML (with image attachments) Error: Could not modify Schema?

texts:
I would like to convert to KMZ with attachments so that colleagues can view the location of photos in Google Earth. I found this Layer to KML tool to use:
https://www.arcgis.com/home/item.html?id=5d8704c938ea4715b59eebabcd96c1d9

I am using ArcGIS Desktop 10.7
I downloaded my data from ArcGIS Online as an FGDB to get the images and point locations together. Each point location has at least 1 photo, some points 3-4 photos. I viewed my ObjectID, verified the points were non-sequential and ensured that the Allow Unique ID field was checked on. When I run the tool, I get the below error.

Do I need to convert to a different file format (currently FGDB from ArcGIS Online)?

In the comments for the tool, the creator specifically states that the tool only requires a basic licence to run, so I am confused by the error message indicating the licensing is not valid.

Failed script Layer2KMLAttach...

Traceback (most recent call last):
File "C:\Users\Katie\Downloads\Toolbox\Layer_to_KML_attachment.py", line 177, in
arcpy.AddField_management(inputFeatures, "tempIDField", "TEXT")
File "c:\program files (x86)\arcgis\desktop10.7\arcpy\arcpy\management.py", line 3429, in AddField
raise e
ExecuteError: ERROR 999999: Error executing function.
The application is not licensed to create or modify schema for this type of data
Failed to execute (AddField).

Failed to execute (Layer2KMLAttach).


texts: Rendering a shapefile with a custom projection in GeoTools

texts:
I have an image, a shapefile, and a method which can translate an x/y pixel location from the image into a latitude/longitude point or vice versa. I'm trying to use Geotools to render the shapefile on the image, but I can't seem to get it to work.

All of the examples and tutorials I can find use factories to retrieve pre-existing transforms and mappings, but how can I use my custom projection?

Edit: The image is the view from a satellite, which means that it cannot see the entire earth. Is Geotools capable of handling this?

Edit2: I don't think that this can be done easily in GeoTools, I ended up just writing what I needed from scratch.


texts:
Not sure if this has been posted yet. But I couldn't find any software that can do this yet.

I am given a large set of random locations and a set of set office location. I would like to map the random locations to the set of office locations based on travel time and efficiency od the route and mapping. So basically, given a time or distance that each office can cover(which can change base on the office), is there anyway to create a most efficient routing?

Any help is greatly appreciated.


texts:
Anyone able to convert ESRI PBF into anything "human readable"? For example GeoJSON?
The ArcGIS dev doc is not very concrete here:

The PBF format is binary and is thus smaller and faster for clients to parse than JSON. It is used by newer versions of some Esri APIs and applications to improve efficiency. The PBF format is supported by some feature service operations in ArcGIS Enterprise 10.7 as well as hosted feature layers in ArcGIS Online.

From the esri Community discussion it seems the format is based on MapBox specification but nothing concrete here. Just a promise they'll work on documenting (since late 2019).
This is an example file: https://services.arcgis.com/V6ZHFr6zdgNZuVG0/arcgis/rest/services/NYC_Footprints_fgdb/FeatureServer/0/query?f=pbf...
Tried various tools but without success.

EDIT: sure there is an option to use f=geojson, f=json or f=pjson but the question is about dealing with pbf format, not the ArcGIS query.


texts: QGIS gives dll errors after installation

texts:
I'm having issues installing QGIS 3.10. I tried uninstalling and reinstalling, with the same errors as a result.

Below you can read the setup.log.

"2019/11/05 16:17:48 Starting OSGeo4W install, version 2.579
2019/11/05 16:17:48 using locales .1252
2019/11/05 16:17:48 io_stream_cygfile: fopen(/etc/setup/setup.rc) failed 2 No such file or directory
2019/11/05 16:17:48 Current Directory: C:\Users\schelsma\AppData\Local\Temp
2019/11/05 16:17:48 Root Directory: C:\OSGeo4W64
2019/11/05 16:17:48 Could not open service McShield for query, start and stop. McAfee may not be installed, or we don't have access.
2019/11/05 16:17:55 source: network install
2019/11/05 16:17:56 root: C:\OSGeo4W64
2019/11/05 16:17:58 Selected local directory: C:\Users\schelsma\AppData\Local\Temp
2019/11/05 16:17:58 Menu name: OSGeo4W
2019/11/05 16:17:59 net: Direct
Cached mirror list unavailable
get_url_to_membuf http://download.osgeo.org/osgeo4w/mirrors.lst
getUrlToStream http://download.osgeo.org/osgeo4w/mirrors.lst
getUrlToStream failed!
get_url_to_membuf failed!
Defaulting to empty mirror list
2019/11/05 16:18:52 source: network install
2019/11/05 16:18:57 source: network install
2019/11/05 16:19:08 root: C:\OSGeo4W64
2019/11/05 16:19:15 Selected local directory: C:\Users\schelsma\AppData\Local\Temp
2019/11/05 16:19:15 Menu name: OSGeo4W
2019/11/05 16:19:16 net: Direct
2019/11/05 16:19:21 io_stream_cygfile: fopen(/etc/setup/setup.rc) failed 2 No such file or directory
2019/11/05 16:19:21 Ending OSGeo4W install"

When starting up QGIS, a series of dll errors are shown. Ending in this last error.

Any help?


texts: Adding feature id column that autoincrements using PyQGIS?

texts:
Is there anyway to add the feature id column which increases automatically in attribute table programmatically? Here is a code to add a column:

fields = []
fields.append(QgsField('SOIL_ID', QVariant.Double))
vLayer.dataProvider().addAttributes(fields)
vLayer.updateFields()

But I want to add a column whose name is 'id' and it should increase in order without my touch.


texts:
I have a Python add-in with a good amount of buttons and menus. I'm trying to remove a menu using the addin_assistant, but the menu isn't being removed from config.xml. I've tried removing the add-in in ArcMap and re-adding it, but that didn't help.

The process I'm doing:

Open the Python Add-in Wizard, browse to my add-in folder.
Find the menu. Right click and Remove.
Save. Click OK for the backup warning message.
The config.xml document date is modified automatically by saving.
Open the XML and see that the menu is not removed.

Quick thoughts:
I've tried restarting my computer.
I've removed buttons successfully.
I can't remake the add-in every time I want to remove a menu.
Customize Mode -> Delete for the menu (in ArcMap) isn't feasible because this add-in will be shared with hundreds of users.

Any thoughts? Please let me know if you need more info.

Edit:

I've copied the contents of my add in folder to a different folder. I tried removing the menu the way I previously described and it didn't work.

I edited the XML in Notepad++ to manually remove the menu and ran the makeaddin.py file, then removed and re-added the add-in in ArcMap. This didn't remove the menu.

Double clicking and installing the esri addin file worked to remove the menu since I've edited the XML, but now when I try to access the add-in directory using the add-in assistant I get an error window that says:

Error initializing Add-In

add-in name.menu_remove  (which was the name of the menu ID).

It seems that the addin_assistant is expecting to find menu, but since it was removed in XML it's throwing an error.

Edit 2

I didn't fully delete all of the required lines in XML. The error thrown above no longer appears. I had to delete both the Menu caption for the menu category and the Menu refID that referenced my menus ID (menu_remove).

Editing the XML works well enough to consider this answered. Esri believes the issue is a bug or limitation with the python addin wizard. I'll add that info when I receive it.


texts: Is it possible to add a raster file as layer from an URL?

texts: I have several geoTIFFs stored on a server that I can access through http. I would like to add them as raster layers analogous to adding a local raster file through the "Add Raster Layer" button, but with an URL as source, without storing the image as a layer beforehand (in PostGIS, for example). Is that possible?


texts: Given a shapefile, how can I return the polygon containing a lat/lon location?

texts:
I have a shapefile generated by the ESRI product suite. I do not have any ESRI products and my budget does not support their purchase. Aside from fiscal issues, there may be technical limitations with deploying ESRI products but that's a moot point since we have no money.

Using .NET and, if necessary, open source libraries, how can I read my shape file data and return the polygon that contains a latitude and longitude coordinate? I see this being built as service code and given a coordinate, it returns the polygon. Conceptually not too difficult but I'm not sure which libraries out there will help with this.

Thoughts?

edit #1 - I'd settle for returning the polygon's ID or description or some other string property. If I want more info, I will fetch that using the polygon's ID or name or ???

edit #2 - The solution should be something I can call from a service. For example, a library I can call from my .NET code would work just fine. An application that must be installed and data entered through its UI will not work since my application will be deployed.


texts:
With single images, addBands works, but what about an imageCollection?

Code for the problem i am talking about

In this code, I have generated drought indices with variable names, VCI, TCI, PCI, VHI, SMCI, and Arid.
All these variables hold imageCollections with single bands. Now I want to compile them into a single imageCollection, but have no clue on how to achieve this goal.


texts: Labeling Selected Features Only using ArcGIS Engine?

texts:
I've been labeling using the LabelEngine and simple expressions and it is working for every layer. I want to offer the option of labeling only selected features and thought it would be this simple:

annotateLayerProperties.LabelWhichFeatures = esriLabelWhichFeatures.esriSelectedFeatures;

Unfortunately that doesn't do it. Every polygon is still labeled. I verified that the features were selected (I have a routine that just returns the selected keys.) I also tried replacing the Engine labeling with a:

annotateLayerProperties.WhereClause = BuildWhere( keyName, selected );

using both my keys and the features OIDs - that does not label at all.

Does anyone know how to label only the selected features?


texts: I got a mesh as Dolfin XML and would like to know if there is anybody that had some experiences with the web visualisation of it. Is it possible to convert in X3D?


texts: Optimizing script with in_memory selection

texts:
I had a problem solved here.

I had 2 features: polygons and points, both containing similar attributes(example: street name, number). Each polygon must have equivalent point added(same street name, same number), and I had to detect which polygons didn't have a point feature(or if they had one and the values are different, because of typos and such). The problem was solved by doing a spacial join and then filtering the ones with different values.

My new problem is that I need to run the script on an selection of features. Here
you can find something that i plan to use. As I see it , before I run  arcpy.MakeFeatureLayer_management I should use the arcpy.CopyFeatures_management  feature, with arcpy.env.workspace = "in_memory" defined.

But what should I grab for the "selected_features" ? In that code, that does "pipes" define? Were does "in_memory" kick in?

import arcpy arcpy.env.overwriteOutput = True

from arcpy import env

arcpy.env.workspace = "in_memory"

# *make 'in memory' selection from featurelayer  arcpy.CopyFeatures_management (in_feat, out_feat)

selected_features = "D:/ALTE PROIECTE/script/cladiri adrese/Adrese_Output.shp"

sel_adrese = "D:/ALTE PROIECTE/script/cladiri adrese/selectadrese.shp"

arcpy.CopyFeatures_management ("D:/ALTE PROIECTE/script/cladiriadrese/Adrese_Output.shp", sel_adrese)

selected_features = "D:/ALTE PROIECTE/script/cladiri adrese/cladiri_Output.shp"

sel_poligon = "D:/ALTE PROIECTE/script/cladiri adrese/selectpoligon.shp"

arcpy.CopyFeatures_management ("D:/ALTEPROIECTE/script/cladiriadrese/cladiri_Output.shp", sel_poligon)

# *Make Feature Layers from memory selection

arcpy.MakeFeatureLayer_management(sel_adrese, "lyr_Adresa")

arcpy.MakeFeatureLayer_management(sel_poligon, "lyr_Poligon")

# *Do a spatial join

inputFC = "lyr_Poligon" joinFC = "lyr_Adresa" outputFC = "D:/ALTE PROIECTE/script/cladiri adrese/JoinPoligonAdresa.shp"

arcpy.SpatialJoin_analysis(inputFC, joinFC, outputFC)

# *Make a new Feature Layer that just contains the ones that need fixed

arcpy.MakeFeatureLayer_management("D:/ALTE PROIECTE/script/cladiriadrese/JoinPoligonAdresa.shp", "lyr_MySpatialJoinLayer", "(\"Strada\" <> \"DenumireAr\") >OR (\"Nr_postal\" <> \"Nr_Posta_1\")")

# *Create a new shapefile (or FGDB) with just the polygons that need fixed

arcpy.CopyFeatures_management("lyr_MySpatialJoinLayer", "D:/ALTE PROIECTE/script/cladiriadrese/Erori.shp")

That code will do the spacialjoin, but not on the features I selected. what I am missing out? I am using arcgis 10.2


texts: Avoiding unexpected raster search with Python script?

texts:
I am trying to write a python script that will find all raster files in my server including 'GRID'.

Here is my sample script:

import os
import arcpy
path = r"F:\Geodaten"
for root, dirs, files in os.walk(path):
    #print ("we are in " + dirs)
    #Accessing to directories...................................

    #accessing files...........................
    for fileName in files:
        basename, extension = os.path.splitext(fileName)
        if (extension == ".jpg") or (extension == ".tif"):
            fullPath = os.path.join(root, fileName)
            #mxd = arcpy.mapping.MapDocument(fullPath)
            #print ("we are in " + dirs)
            print "bilder: " + fullPath
walk = arcpy.da.Walk(path, type="GRID")

for root, dirs, files in walk:
    for filename in files:

        if "." not in filename:
            fullPath=os.path.join(root, filename)
            print "Grid: " + fullPath

print("......searching Done...............")

my result is showing some unexpected grid files that I could not access
like this:

  Grid: F:\Geodaten\Eigene\GPS\Laukesa2008\shp\20080529.shp\20080529
  Grid: F:\Geodaten\Eigene\GPS\Laukesa2008\shp\20080530.shp\20080530
  Grid: F:\Geodaten\Eigene\GPS\Laukesa2008\shp\20080531.shp\20080531
  Grid: F:\Geodaten\Eigene\GPS\Laukesa2008\shp\20080602.shp\20080602
  Grid: F:\Geodaten\Eigene\GPS\Laukesa2008\shp\20080603.shp\20080603
  Grid: F:\Geodaten\Eigene\GPS\Laukesa2008\shp\20080604.shp\20080604
  Grid: F:\Geodaten\Eigene\GPS\Laukesa2008\shp\20080605.shp\20080605

How could I avoid unexpected results?


texts: Postgis error "unknown spatial reference system" on ST_GeomFromGML

texts:
I'm trying to create a Postgis geometry from GML using ST_GeomFromGML. It results in an error message "unknown spatial reference system".

The GML has been created by Oracle function sdo_util.to_gml311geometry and looks like this:

<gml:Polygon srsName="SDO:28992" xmlns:gml="http://www.opengis.net/gml">
   <gml:exterior>
      <gml:LinearRing>
         <gml:posList srsDimension="2">...</gml:posList>
      </gml:LinearRing>
   </gml:exterior>
</gml:Polygon>

srsName SDO:28992 refers to the Oracle SRID 28992. What should I put in the GML so that PostGis understands this?


texts: What tool breaks up raster cells into individual values based upon contiguity?

texts:
The green values have a COUNT attribute of 13418.
I will want to use a tool that lumps the COUNT attribute and reclassifies based on values. But first it appears I will need to break up the raster cells into individual values using a tool? I am learning raster, in the past I only used vector.


texts: Getting polylines back from leaflet when using esri-leaflet featureLayer

texts:
I have a leaflet map that uses esri-leaflet to display a feature layer created for arcgis that is a number of paths. These paths are created correctly and I can see them on the map, but they aren't the visual I need. Instead I'm interested in getting the LatLng of the centroid of the path.
Since I know leaflet polylines have a getCenter method, I just need to know how to access these objects after they are drawn to the map. Is this possible?
Code for this follows:
let mymap = L.map("mapid").setView([lat, lng], 18);
L.esri.basemapLayer("Gray").addTo(mymap);
var subtleUA = L.esri
  .tiledMapLayer({
    url:
      "https://mymapserver.../"
  })
  .addTo(mymap);
let roomFeatures = L.esri.featureLayer({url:"https://mymapserver/../interiors/16"}).addTo(mymap)


texts:
I am not familiar with gis or writing sql, so my understanding in this field is lacking. I have a hundreds of lat/lng geographic coordinates stored in a single table called Places like this:

ID | Coordinate
-----------------------------
1  |  POINT(192.938 -28.348)
2  |  POINT(57.349 -88.482)

They are in 4326, and I'm using geography type. I want to find all nearby points that are within 10 miles of a given point.

I want to use st_dwithin(), but I don't understand how to use it. I came up with this, but it doesn't seem quite right:

SELECT * FROM places WHERE ST_DWithin(places.Coordinate, ST_GeographyFromText('POINT(195.888 -28.348)'), 16093.4);

16093.4 metres is 10 miles.

How do I get all points within 10 miles of a given point (in this example, from POINT(195.888 -28.348))?


texts: Getting average line of two non-parallel lines using ArcGIS Desktop?

texts:
I have two nonparallel lines and I need to find the average line between them.  This is proving to be more difficult than expected.

Is there a fairly simple way to do this?

I thought I could use a series of points for each line, then find the mid point of each set of points then created a line using said points, however this didn't work has I had hoped.  The lines are not the same length fyi.  I do have an ArcGIS Advanced license.


texts: Creating Minimum Bounding Polygon in QGIS?

texts:
Does anybody know how to create a "minimum bounding polygon" from a set of polygons in QGIS or elsewhere?

For example, in the figure below I can create a Convex Hull per the purple outline, but I'm after a tighter fit that matches the inner reddish shape.  I can do a buffer with a small distance, but there are interior holes that are tedious to delete manually.


texts:
I have a raster stack of 108 elements

The stack 'images' is created using 4 bands (RGB+NIR) of 27 satellite images. The sequence of the rasters in the stack is:

Position [[1]]: blue band image 1

Position [[2]]: green band image 1

Position [[3]]: red band image 1

Position [[4]]: NIR band image 1

Position [[5]]: blue band image 2

...

Position [[108]]: NIR band image 27

I want to derive two basic vegetation indices for each image (27 NDVI and 27 NDWI2 indices): NDVI - (NIR-R)/(NIR+R) and NDWI2 (G-NIR)/(G+NIR)

So far I have doing it manually. For example:

NDVI_1-> (images[[4]] - images[[3]]) / (images[[4]] + images[[3]])

NDVI_2-> (images[[8]] - images[[7]]) / (images[[8]] + images[[7]])

and

NDWI2_1-> (images[[2]] - images[[4]]) / (images[[2]] + images[[4]])

NDWI2_2-> (images[[6]] - images[[8]]) / (images[[6]] + images[[8]])

Since my raster stack as a pattern, I am looking for an automatic way to derive the indices. Something like a loop that calculates the indeces, save it with the appropriate name and then moves to the next 'position' to repeat the process.

Any help on that?


texts: How to save a Leaflet or OpenLayers map for opening later?

texts:
Can a webmap state be saved to a text file or something for restoring later?
Like if I set up a map with a particular zoom and choose my layers and draw a polygon or two, can it be saved to a file and reopened in the same state later on?

Thanks


texts: White squares in WMTS for OS Raster in QGIS

texts:
A WMTS connection I have had in QGIS for some time now is currently not working as expected. I have not had this issue before and I have not changed the URL connection.

The connection in question is for OS Open Raster from ESRI's Vector Tile Services and the connection URL is:

  http://tiles.arcgis.com/tiles/qHLhLQrcvEnxjtPr/arcgis/rest/services/OS_Open_Raster/MapServer/WMTS/1.0.0/WMTSCapabilities.xml?cacheKey=99417257241d59ec

Here is A GIF to highlight the issue: