texts:
I am currently trying to calculate the angle between two lines or rather vectors (which are not touching) like in the picture shown below:

In order to calculate the angle α, I think I will have to lengthen v1 somehow, so I get an intersection point of both vectors, like this:

Is there a way to do this within a python script (and how?) or can I do this with the field calculator of QGIS?

As I have two vector layers, I need this calculation more than one time and I will have to store the result into an attribute.


texts: Purely client-side geocoding options for a webmap

texts:
I am creating a simple webmap that should be able to:

Display external maptiles (any of the following: Cloudmade, Mapquest, Google).
Overlay additional static layer in Geojson or KML format, that can show a popup box on a mouse-click.
Use a search box for geocoding that is limited only to towns and cities.  And then zoom in to the geocoded location at a specific zoom level (z8-12 in webmap speak).

I'm looking at several webmapping framework like Openlayers, Leaflet, Polymaps.  So far, I was able to do 1 and 2.

I'm having difficulty integrating the third feature.  There are several geocoding APIs available (OSM's Nominatim, Geonames and Google's geocoder).  However, I can't find any that limits geocoding to towns/cities.

A major challenge is to do this purely on the browser/javascript, no server side processing since the webserver is only limited to static htmls.


texts: Projection - Some of my data are in odd position

texts:
Some of my data has this position lon: 181 / lat: 91
So when I try to draw them on the Google map:

var geographic = new OpenLayers.Projection("EPSG:4326");
var mercator = new OpenLayers.Projection("EPSG:3857");

var wfsProtocol = new OpenLayers.Protocol.WFS.v1_1_0({

               url: GeoServerURL,
               featureType: "Vector_Layer",
       featureNS: "http://mapmap.org",
       srsName: "EPSG:3857",
               version: "1.1.0",
               extractAttributes: true,
               isBaseLayer: false,
               visibility: true
           });

wfsProtocol.read({
    filter: filter,
    callback: processTheQuery,
    scope: [new OpenLayers.Strategy.BBOX({
                   ratio: 1
               })],
  })

Layer = new OpenLayers.Layer.Vector("test", {
maxExtent: new OpenLayers.Bounds(-180, -85, 180, 85).transform(geographic,mercator),
displayOutsideMaxExtent: true,
displayInLayerSwitcher: false
         })

function processTheQuery(request) {
if (request.data && request.data.bbox) {
var b = request.data.bbox;
sExt = new OpenLayers.Bounds(b[0],b[1],b[2],b[3]);
Layer.destroyFeatures();
Layer.addFeatures(request.features);
map.zoomToExtent(sExt)

I got this error:

<?xml version="1.0" encoding="UTF-8"?>
<ows:ExceptionReport version="1.0.0"
  xsi:schemaLocation="http://www.opengis.net/ows http://myDomain/geoserver/schemas/ows/1.0.0/owsExceptionReport.xsd"
  xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" xmlns:ows="http://www.opengis.net/ows">
  <ows:Exception exceptionCode="NoApplicableCode">
    <ows:ExceptionText>java.lang.IllegalStateException: A transformation exception occurred while reprojecting data on the fly
A transformation exception occurred while reprojecting data on the fly
Latitude 91°00.0'N is too close to a pole.</ows:ExceptionText>
  </ows:Exception>
</ows:ExceptionReport>


texts:
I am migrating my Flex eSearch widget to the CMV Search widget.  I need to be able to pass a query through an url within the widget.  Currently my widget has 8 URLs that are in the widget.
Flex Code:


texts: Error when using UPDATE to populate a geometry column

texts:
I'm trying to update a column with new geometry but am getting

ERROR:  column "current_geom" is of type geometry but expression is of type boolean.

UPDATE this_table t1
    SET "current_geom" = ST_Intersects(t2.point_geom,t2.geom)
FROM other_table t2
WHERE t1.current_geom is null

I'm obviously doing something wrong but not sure why it doesn't accept it?

Using: PostgreSQL 9.6 and pgAdmin 4.2


texts: Error while creating spatial index on geogarphy column in sql server

texts:
I used the following code to create geography column using latitude, longitude

alter table [locate].[dbo].[BFLGeocodedResult] add geographyColumn as geography::STGeomFromText
('POINT('+convert(varchar(20),[Final Longitude])+' '+convert(varchar(20),[Final Latitude])+')',4326)

Also created Clustered Index and then tried to create spatial index using the below query

       CREATE SPATIAL INDEX [BFLGeocodedResult_sidx] ON [locate].[dbo].[BFLGeocodedResult]
(
    geographyColumn
)USING  GEOGRAPHY_GRID
WITH (
BOUNDING_BOX =(-90, -180, 90, 180))

when I run the above query I get the error

Cannot create primary xml, selective xml or spatial index 'BFLGeocodedResult_sidx' on table 'locate.dbo.BFLGeocodedResult', column 'geographyColumn', because the column is computed.


texts: Allow user to add notes to an OpenLayers map

texts: I am looking for any ideas here. I want a user to be able to interact with a map. I know I can use HTML5 canvas to allow the user to draw over a map. But what about if I want something simpler? For instance, the user can click on the map, which will add a textbox and pinpoint - the user can then write notes for any locations they place this point. Is there a library which will allow me to do this? Thanks!


texts:
Using QGIS 3.8.0. Despite similar questions out there I can't reproject a vector layer to a geopackage while specifying the original table name in a qgis python processing script. FYI all tablenames have underscores, so no spaces or invalid characters.

I see two options:

processing.run('native:reprojectlayer', alg_params....
QgsVectorFileWriter.writeAsVectorFormat(vlayer, geopP,'utf-8', crsObject, "gpkg")

Option 1: the native:reprojectlayer, I can't see how to specify the table name in the alg_params:

alg_params = {
            'INPUT': parameters['inputlayer'],
            'TARGET_CRS': parameters['selectcoordinatesystem'],
            'OUTPUT': outData}

When I specify the tablename for the geopackage in 'outData', like c:/temp/geopackage.gpkg|layername=buildingData
it throws an error:

Could not create layer c:/temp/MyExport.gpkg|layername=buildingData.gpkg: Creation of data source failed (OGR error: sqlite3_open(c:/temp/MyExport.gpkg|layername=buildingData.gpkg) failed: unable to open database file)

So it looks like the 'native:reprojectlayer' output is looking for a destination file ending '.gpkg' and has added that on the end of my table name. Does this algorithm support writing specific tablenames?

The only thing that works is only passing in the gpkg without the table and that repojects the dataset but exports the table with the same name as the geopackage file (c:/temp/MyExport.gpkg|layername=MyExport).

Option 2: Trying QgsVectorFileWriter I'm using vector layer objects taken from the selected layer in the QGIS project:

def initAlgorithm(..

self.addParameter(QgsProcessingParameterVectorLayer(self.selectdataV, 'selectdataV', types=[QgsProcessing.TypeVector], defaultValue=None))

def processAlgorithm(...

layerV = self.parameterAsVectorLayer(parameters, self.selectdataV, context)

layname = layerV.name()

When using SaveVectorOptions() I want to specify layername.
options = QgsVectorFileWriter.SaveVectorOptions()
options.actionOnExistingFile = QgsVectorFileWriter.CreateOrOverwriteLayer
options.layerName = layname

writer = QgsVectorFileWriter.writeAsVectorFormat(layer=layObj, fileName=geopP, layerOptions=options, fileEncoding='utf-8', destCRS=crsObject, driverName="gpkg")

I get the error:
TypeError: QgsVectorFileWriter.writeAsVectorFormat(): arguments did not match any overloaded call:
overload 1: argument 'layerOptions' has unexpected type 'SaveVectorOptions'
overload 2: 'layerOptions' is not a valid keyword argument

If I take out layerOptions it will just create the same table name as the Geopackage name like in Option 1.

What am I doing wrong? baring in mind I need to specify the new CRS and table name?


texts: Extracting pixel values by overlaying polygons

texts:
I am trying to extract pixel values by overlaying polygons. I use a code after Patrick Grey (http://patrickgray.me/open-geo-tutorial/chapter_5_classification.html). When I masked the image with the shape features I got out_image. Then the next step is removing 0, which totally mess up the array as values are not present according bands. I have been coding only for a month and I tried many ways to remove 0 and keep the order of band values according class. In R I can do it without any problem and when I export the data as CSV and train the algorithm everything works fine in Python environment.
How can I extract pixel values and keep the numbers band and class-wise?
    X = np.array([], dtype=np.int8).reshape(0,8) # pixels for training
    y = np.array([], dtype=np.string_) # labels for training

   with rasterio.open(img_fp) as src:
        band_count = src.count
        for index, geom in enumerate(geoms):
            feature = [mapping(geom)]

    # the mask function returns an array of the raster pixels within this feature
    out_image, out_transform = mask(src, feature, crop=True)
    # eliminate all the pixels with 0 values for all 8 bands - AKA not actually part of the shapefile
    out_image_trimmed = out_image[:,~np.all(out_image == 0, axis=0)]
    # eliminate all the pixels with 255 values for all 8 bands - AKA not actually part of the shapefile
    out_image_trimmed = out_image_trimmed[:,~np.all(out_image_trimmed == 255, axis=0)]
    # reshape the array to [pixel count, bands]
    out_image_reshaped = out_image_trimmed.reshape(-1, band_count)
    # append the labels to the y array
    y = np.append(y,[shapefile["Classname"][index]] * out_image_reshaped.shape[0])
    # stack the pizels onto the pixel array
    X = np.vstack((X,out_image_reshaped))


texts: Pixel values change after clipping a raster layer

texts:
I'm trying to evaluate mean temperatures for the distribution area of European fish species. Therefore, I downloaded the E-OBS data set for mean annual temperature in Europe (0.25 deg. reg., "TG"), loaded it into QGIS and clipped the raster using a polygon shape as a mask.
I did not alter the projection, but I added the command
--config GDALWARP_IGNORE_BAD_CUTLINE YES
to the console:
After that, I used the value tool to extract pixel values for random locations of the newly clipped layer and noticed, that strongly deviate from the original data. Unfortunately, I can't post more that 2 links here, but the original data is in °C values and the clipped layer has values in the magnitude of plus-minus several hundreds.

Is there an obvious flaw in my workflow I'm not seeing?

edit: Text fraction of the Layer obtained from Layer properties->Metadata->Properties:

  D:/GoogleDrive/Sonstiges/GIS/tg_0.25deg_reg_1950-1964_v15.0.nc/tg_0.25deg_reg_1950-1964_v15.0.nc
  latitude#long_name=Latitude values
  latitude#standard_name=latitude
  latitude#units=degrees_north
  longitude#long_name=Longitude values
  longitude#standard_name=longitude
  longitude#units=degrees_east
  NC_GLOBAL#Conventions=CF-1.4
  NC_GLOBAL#Ensembles_ECAD=15.0
  NC_GLOBAL#history=Wed May 31 15:20:56 2017: ncks -a -d time,0,5478 tg_0.25deg_regular_1.nc tg_0.25deg_reg_1950-1964_v15.0.nc Wed May 31 15:19:16 2017: ncks -a --mk_rec_dmn time tg_0.25deg_regular.nc tg_0.25deg_regular_1.nc
  NC_GLOBAL#NCO="4.6.2"
  NC_GLOBAL#References=http://www.ecad.eu\nhttp://www.ecad.eu/download/ensembles/ensembles.php\nhttp://www.ecad.eu/download/ensembles/Haylock_et_al_2008.pdf
  NETCDF_DIM_EXTRA={time}
  NETCDF_DIM_time_DEF={5479,4}
  NETCDF_DIM_time_VALUES={0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,43,44,45,[...]


texts: I have GPS data - point cloud in a shapefile. Now I want to average this data and get on output a single GPS track. Here I found images but can't find algorithm. Maybe gdal or postgis can do it? Any advise?


texts: Running zonal statistics on multiple raster tiles using QGIS?

texts:
I'm new to QGIS (2.18.11). I've got OS Terrain 50 elevation data, in around 150 ASC tiles (10 km x 10 km) covering my whole area of interest (north-east Scotland). There is just a single band. I've got a grid of 1 km squares for my area of interest (made using TomBio plugin).

I'd like to extract mean elevation values per 1 km square (I may look for further stats later but mean is fine for starters). I've used Zonal Statistics to do this on one tile at a time and it adds a field to the Table of Attributes of the TomBio grid. That's all fine, and just what I need - except, is there a quicker way to do it for all 150 tiles at the same time?


texts: QGIS WMS URL parameters for zoom / scale?

texts:
I have a link for WMS in QGIS Server like this:

http://localhost/qgis/qgis_mapserv.fcgi.exe?SERVICE=WMS&VERSION=1.3.0&CRS=EPSG:21781&REQUEST=GetMap&map=C:/OSgeo4W/SERVER/SQL.qgs&BBOX=608350,264875,615466,268677&WIDTH=1200&HEIGHT=1300&LAYERS=TS&FILTER=TS:%22ID%22+=+3&FORMAT=image/png&

Is there any possibility to give an additional parameter for zooming in URL like:

zoom-order=n& or scale=n& (not for printing)


texts:
This is the first real world python program am trying to write. I am trying to create a scrip that takes a shapefile that has n polygons with an integer formatted  date (e.g. yyyymmdd or 6/24/2014 is 20140624) attribute and returns a shapefile which has a single polygon for each date, and that polygon is the summation of all polygons with that date or earlier. Additionally I need the polygon to retain the the attributes of the most recent polygon.

My initial thought was to manually sort the attribute table form most recent to oldest and have something like this to do the rest:

import arcpy
form arcpy import env
import shapefile

w = shapefile.Writer()
# Create a list of all off the features in the shapefile.
def generate_working_list(a):
    env.workspace = "a"
    wl = arcpy.ListFeatureClasses()

# copy the oldest polygon, merge it with the next oldest, copy the newly created polygon, repeat
def aggregate(wl):
    twl = []
    fl = []
    wl[-1].append(twl)
    wl[-1].append(fl)
    wl.pop()
    while len(wl) > 1:
         wl[-1].append(twl)
         '''Dissolve/ merge twl'''
         twl[0].append(fl)
         wl.pop()
    swl[0].append(twl)
    arcpy.Append_management(twl, fl) '''I think that this is the right tool to use'''
    twl[0].append(fl)
    return fl
    w.save("Merged")


texts: Stumped on python looping and data dictionaries

texts:
I have been stumped getting my python code to work for a while now. I've been using python for the last several years on and off for simple things, but I am far from an expert in python. I am trying to create a loop that cycles through specific folders and reprojects all rasters in each particular folder into a new specified empty folder in a different location, then I need the program to do the same for the next specified folder and so on. Below is my stripped down code with my specified input and output folder locations. I was thinking a data dictionary would be best here but I couldn't get the dictionary to run with the loop. Any help would be greatly appreciated!

Thanks so much in advance!

# Import arcpy module
import arcpy
from arcpy import env

# Load required toolboxes...
arcpy.CheckOutExtension("3D")
arcpy.CheckOutExtension("spatial")

arcpy.env.workspace = "z:\\z11_10m\\"

input_wksp = ["z:\\z11_10m\\s41120",
              "z:\\z11_10m\\s42120",
              "z:\\z11_10m\\s42121",
              "z:\\z11_10m\\s43120",
              "z:\\z11_10m\\s43121",
              "z:\\z11_10m\\s44120",
              "z:\\z11_10m\\s44121"]

out_wksp = {"z:\\z11_10m\\project\\s41120\\",
            "z:\\z11_10m\\project\\s42120\\",
            "z:\\z11_10m\\project\\s42121\\",
            "z:\\z11_10m\\project\\s43120\\",
            "z:\\z11_10m\\project\\s43121\\",
            "z:\\z11_10m\\project\\s44120\\",
            "z:\\z11_10m\\project\\s44121\\"}

for rast in arcpy.ListRasters("*"):
    arcpy.ProjectRaster_management(rast, out_wksp + rast, "PROJCS['NAD_1983_UTM_Zone_11N',GEOGCS['GCS_North_American_1983',DATUM['D_North_American_1983',SPHEROID['GRS_1980',6378137.0,298.257222101]],PRIMEM['Greenwich',0.0],UNIT['Degree',0.0174532925199433]],PROJECTION['Transverse_Mercator'],PARAMETER['False_Easting',500000.0],PARAMETER['False_Northing',0.0],PARAMETER['Central_Meridian',-117.0],PARAMETER['Scale_Factor',0.9996],PARAMETER['Latitude_Of_Origin',0.0],UNIT['Meter',1.0]]", "NEAREST", "10", "", "", "PROJCS['NAD_1983_UTM_Zone_10N',GEOGCS['GCS_North_American_1983',DATUM['D_North_American_1983',SPHEROID['GRS_1980',6378137.0,298.257222101]],PRIMEM['Greenwich',0.0],UNIT['Degree',0.0174532925199433]],PROJECTION['Transverse_Mercator'],PARAMETER['False_Easting',500000.0],PARAMETER['False_Northing',0.0],PARAMETER['Central_Meridian',-123.0],PARAMETER['Scale_Factor',0.9996],PARAMETER['Latitude_Of_Origin',0.0],UNIT['Meter',1.0]]")
    print rast, "projected and extracted"

print "finished"


texts: Disconnecting QGIS from internet via settings

texts: Recently QGIS is always very slow (or better freezing) when I am connected to the internet. Some processes take several minutes (even when opening or clicking on a dropdown menu from the taskbar). As soon as I disconnect my internet everything works fine and as usual. Is there an option in the settings to deny QGIS access to the internet. I have already deactivated the display of news and the search for updates. I use QGIS 3.16.4.


texts:
I am trying to come up with a way to export an attribute table as a formatted excel table in ArcMap 10.2. I am aware that there is a Table To Excel tool in the conversion toolbox, but I am not sure if there is more functionality with Excel that can be accessed using Python.

Right now I have a macro in Excel that reads in another Excel table and sorts and formats it, but I would like to be able to do this all in one moodel or script within the ArcGIS Desktop environment.

The first screenshot is the exported excel table and the second is the formatted table after running the macro.


texts: Download ALL data for an area from OpenStreetMap

texts:
I am trying to download ALL the OpenStreetMap data for a certain region.  Meaning, I want the roads, the stores, the rivers, the parks, the traffic signals, the bus stops, ... everything.  I'm expecting that results to be big, so I plan do it for a specific administrative area at a time, and sweep the relevant admin areas.  However, if there is a way to retrieve and download all the data for a larger region, I'm interested in that too. I know how to get all the data of one type by using the Overpass-Turbo wizard; for example, to get all the "highway" data in Shinjuku, this worked:
[out:json][timeout:25];
{{geocodeArea:Shinjuku}}->.searchArea;
// gather results
(
  // query part for: “highway=*”
  node["highway"](area.searchArea);
  way["highway"](area.searchArea);
  relation["highway"](area.searchArea);
);
out body;
>;
out skel qt;

Ideally I'd be able to enter "*" in place of "highway" and get everything, but that didn't work.  I could sweep the types, but I can't find a list of all available element types, so I can't do a "highway or station or ..." type command manually either.

How do I select all kinds of elements?

My current plan is to incorporate that into a Python-based query to get, format, and save the data using a script like the one below.  It's basically the same as the Overpass Turbo wizard query, but note the last line.
tokyo_areacode = 3601543125
query = f"""
    (
        // query part for: “railway=station”
        node["railway"="station"](area:{tokyo_areacode});
        way["railway"="station"](area:{tokyo_areacode});
        relation["railway"="station"](area:{tokyo_areacode});
    );
    // print results
    out body;
    >;
    out skel qt;
    """
    api = overpass.API(timeout=900)
    result = api.get(query, responseformat='csv(::id,::lat,::lon,name)')

Because I'm getting all the data, I want to save all the data too, and I think JSON is best for me in this case (I'm going to load it into NetworkX).

How do I write the responseformat to save all the tags for all the elements as a JSON?

It may seem like two questions, but it's just one question about how to get ALL the data for an area, and I highlighted the two steps in that process where I can't find any solutions.


texts: Adjusting buffer between arrows and both begin and end points separately

texts:
I'm currently making a map of the Netherlands in which the flow of student migrations upon or shortly after their graduation is depicted. Here's what I've made so far:

I'm getting a bit more pleased with the result, but I'm not quite there yet. Something I'd like to adjust in particular is the buffering space of the arrows. In particular, I'd like to change the buffering between the end points of the arrows and their destinations only.
So for the the starting points, I'd like to keep everything the same. I'm okay with how it looks at the Southernmost arrow, for instance:

When we look the the destinations of the arrows, however, things get more messy. They're lumped together at the exact same location:

In order to create some more space between the arrows at the end points, I tried adjusting the I've found in this blog post by Anita Graser. The code is used in the Geometry Generator of the virtual layer of the arrows:
difference(
    difference(
        make_line(start_point($geometry), centroid(offset_curve($geometry, length($geometry)/(-5 + log(3.5,intensity)/1.9))), end_point($geometry)),
        buffer(start_point($geometry), 0.01)
        ),
    buffer(end_point($geometry),0.5)
    )

Adjusting the values of the buffer at the start_point and end_point does not seem to have any effect on the actual buffering space of the arrows on the map, though. So my question is: how to change the buffering space between the arrowheads and their destination locations to a suitable value?


texts: QGIS does not recognize the CRS when the data is provided by SQLServer Database with Esri SDE

texts:
We have a MS SQL Server: Microsoft SQL Server 2016 database with SDE from Esri on Windows Server 2012 R2. Loading data into QGIS results in the QGIS unkown CRS dialog. If we don't select one CRS the layer is displayed with the questionmark for the unknown CRS. But in the layer properties the correct CRS is displayed anyway (QGIS: Version 3.12.2 and 3.16). Does someone knows a solution or workaround?

and


texts: Installing qgis2web in qgis3 (2.99) on ubuntu 16.04

texts:
When installing the plugin qgis2web through the plugin manager on the qgis 2.99 nightly release I get the following error message. Does anyone know what I need to do to get it working. I've not seen an obvious corresponding issue on the github repo.

Couldn't load plugin 'qgis2web' due to an error when calling its classFactory() method

ImportError: No module named 'PyQt5.QtWebKitWidgets'
Traceback (most recent call last):
  File "/usr/lib/python3/dist-packages/qgis/utils.py", line 336, in startPlugin
    plugins[packageName] = package.classFactory(iface)
  File "/home/will/.local/share/QGIS/QGIS3/profiles/default/python/plugins/qgis2web/__init__.py", line 28, in classFactory
    from .qgis2web import Qgis2Web
  File "/usr/lib/python3/dist-packages/qgis/utils.py", line 664, in _import
    mod = _builtin_import(name, globals, locals, fromlist, level)
  File "/home/will/.local/share/QGIS/QGIS3/profiles/default/python/plugins/qgis2web/qgis2web.py", line 30, in
    from qgis2web.maindialog import MainDialog
  File "/usr/lib/python3/dist-packages/qgis/utils.py", line 664, in _import
    mod = _builtin_import(name, globals, locals, fromlist, level)
  File "/home/will/.local/share/QGIS/QGIS3/profiles/default/python/plugins/qgis2web/maindialog.py", line 55, in
    from qgis.PyQt.QtWebKitWidgets import QWebView, QWebInspector, QWebPage
  File "/usr/lib/python3/dist-packages/qgis/utils.py", line 664, in _import
    mod = _builtin_import(name, globals, locals, fromlist, level)
  File "/usr/lib/python3/dist-packages/qgis/PyQt/QtWebKitWidgets.py", line 26, in
    from PyQt5.QtWebKitWidgets import *
  File "/usr/lib/python3/dist-packages/qgis/utils.py", line 664, in _import
    mod = _builtin_import(name, globals, locals, fromlist, level)
ImportError: No module named 'PyQt5.QtWebKitWidgets'

Python version: 3.5.2 (default, Nov 23 2017, 16:37:01) [GCC 5.4.0 20160609]
QGIS version: 2.99.0-Master Master, 585a4d3

Python Path:
/home/will/.local/share/QGIS/QGIS3/profiles/default/python/plugins/gml_application_schema_toolbox/extlibs
/usr/share/qgis/python
/home/will/.local/share/QGIS/QGIS3/profiles/default/python
/home/will/.local/share/QGIS/QGIS3/profiles/default/python/plugins
/usr/share/qgis/python/plugins
/usr/lib/python35.zip
/usr/lib/python3.5
/usr/lib/python3.5/plat-x86_64-linux-gnu
/usr/lib/python3.5/lib-dynload
/home/will/.local/lib/python3.5/site-packages
/usr/local/lib/python3.5/dist-packages
/usr/lib/python3/dist-packages
/home/will/.local/share/QGIS/QGIS3/profiles/default/python


texts:
I have created 3 custom tools, using vb.net, to run in ArcGIS. They all runs fine.They are inside a same project, in three different tool classes . Now i want them all to be added in a single toolbar, so that it will be easier to access them.

To do so, i am using the basetoolbar Item, inside Extending ArcObjects

here is little bit of  whats inside the class

Public Sub New()

'
'TODO: Define your toolbar here by adding items
'
 AddItem("CustomTool_Openlayer")
 AddItem("CustomTool_selectfeature")
 AddItem("CustomTool_splitline")

'BeginGroup() 'Separator
'AddItem("{FBF8C3FB-0480-11D2-8D21-080009EE4E51}", 1) 'undo command
'AddItem(New Guid("FBF8C3FB-0480-11D2-8D21-080009EE4E51"), 2) 'redo command
End Sub

Public Overrides ReadOnly Property Caption() As String
    Get
        'TODO: Replace bar caption
        Return "My VB.Net Toolbar"
    End Get
End Property

after i run the program, i get the MY VB.net Tool bar in the arcmap, but it looks something like this, and none of the tools work.

this is the small description of the tool

MyBase.m_category = "mytools" 'Localizable text.
        MyBase.m_caption = "Open layer And Select Feature" 'Localizable text.
        MyBase.m_message = "Open layer And Select Feature" 'Localizable text.
        MyBase.m_toolTip = "Open layer And Select Feature" 'Localizable text.
        MyBase.m_name = "CustomTool_Openlayer"

someone tell , what to write in the additem() function, i aslo tried with the GUIDs of the tools

Public Const ClassId As String = "4f76ddec-f1b4-472e-90e2-97257c9c23a2"
    Public Const InterfaceId As String = "abdd6064-2ed1-47d4-8af5-fe4836515927"
    Public Const EventsId As String = "ac56be15-367c-4082-966f-129871b28438"

should i use this basetoolbar item, or is there some other way. i just read that in basetoolbar commands are added only not tools. ?? please correct me.


texts: How to treat bias? Calculate zonal statistics proportionally?

texts:
Using ArcGis 10.4

I don't know if this is even possible, but I want to give it a try. The following problem: I have subcatchments (built from a DEM) and I have a river shape. This river shape is divided into 100m stretches with every stretch having many values (environmental data). Important: Within these 100m stretches there is no subdivision. If there are for example 10 found species, we dont know where. And that is fine, this is the accuracy we get.

Now, I need this data in my small catchments (sometimes the sum, sometimes STD or whatever), so I am using zonal statistics. But obviously the 100m stretches can be within more than one subcatchment. If I just transform it to a raster, then I get very wrong results, because every single pixel gets the whole values of the 100m stretch (hence, the sum values in the catchments get too large). My idea was to use "feature to point" to only have this value at ONE point. But then the problem occurs, that the middle point of a line is not necessarily in the catchment with the largest share of the 100m. (See picture, using the feature to point tool and then zonal statistics would mean that all values are assigned to the blue marked catchment, whereas in reality its very very unlikely all of the recorded data was found on that tiny stretch, but rather in the catchments left and right).

A little gerneralization is absolutely fine, since we can't really know where on that 100m stretches values are distributed anyway. But this way it seems very biased. I hope it gets clear what I mean. Any ideas how to improve that??


texts:
I'm having problems trying to compile (r.landscape.evol and another) on the 7.6.1 version of GRASS from the g.extension command. It does not found the file in the default repository. I am running Windows 10. It´s seems the repository does not work.

g.extension extension=r.landscape.evol
Downloading precompiled GRASS Addons <r.landscape.evol>...
ERROR: Extension <r.landscape.evol> not found


texts: Leaflet EasyPrint plugin only prints left half of screen

texts:
I'm trying to print maps from my Leaflet app using the Leaflet EasyPrint. Installation is easy, I made it run but it only prints the left half of the screen, not the entire window.

I see that one can only configure: title, position of button and select the html elements to hide.

Does anyone have experience with the EasyPrint plugin and can advise how to configure it to print the entire screen?


texts: How to reset QGIS installed on Windows to 'first run' state?

texts:
The title sums it up, I'm trying to reset QGIS to its first run state under Windows 7 and XP. Deleting the .qgis folder in my user directory goes a long way to achieving this, but certain things persist such as:

The load state of core plugins
Recent Projects list
Options

And presumably others besides this list

Where are these bits of information stored?

Using QGIS 1.7.4 (411aff6)


texts: Converting raster pixels to polygons with GDAL python

texts:
I want to convert raster pixels to vector polygons as below image that is output of "Processing Toolbox--> Vector Creation--> Raster pixels to polygons" in QGIS 3.8:

I tried to do this vectorize with below code:

from osgeo import gdal, ogr, osr
import sys
import os

#gdal.UseExceptions()
os.chdir("H:/aa")
fileName = "H:/aa/image.tif"
src_ds = gdal.Open(fileName)
srs = osr.SpatialReference()
srs.ImportFromWkt(src_ds.GetProjection())
if src_ds is None:
    print('Unable to open %s' % src_fileName)
    sys.exit(1)
srcband = src_ds.GetRasterBand(1)
print (srcband)
dst_layername = "H:/aa/aapoly"
drv = ogr.GetDriverByName("ESRI Shapefile")
dst_ds = drv.CreateDataSource(dst_layername + ".shp")
dst_layer = dst_ds.CreateLayer(dst_layername , srs = srs, geom_type=ogr.wkbMultiPolygon)
newField = ogr.FieldDefn('DN', ogr.OFTReal)
dst_layer.CreateField(newField)
gdal.FPolygonize(srcband, None, dst_layer, 0, [], callback=None)
dst_ds.SyncToDisk()
dst_ds=None

But output of this code is:

In second image, pixels with same value that are neighborhood were merged in vector output, but I don't want to merge them. First image is desirable but second image as output of the code is not desirable.
It is important to do it with GDAL code and not another tools such as in Convert Raster to vector - create polygons based on each pixel.

With another class, "Polygonize", its result is more different from what I expected (Based on @RafDouglas suggestion in answer section):

from osgeo import gdal, ogr, osr
import sys
import os

#gdal.UseExceptions()
os.chdir("H:/aa")
fileName = "H:/aa/image.tif"
src_ds = gdal.Open(fileName)
srs = osr.SpatialReference()
srs.ImportFromWkt(src_ds.GetProjection())
if src_ds is None:
    print('Unable to open %s' % src_fileName)
    sys.exit(1)
srcband = src_ds.GetRasterBand(1)
print (srcband)
dst_layername = "H:/aa/aapoly"
drv = ogr.GetDriverByName("ESRI Shapefile")
dst_ds = drv.CreateDataSource(dst_layername + ".shp")
dst_layer = dst_ds.CreateLayer(dst_layername , srs = srs, geom_type=ogr.wkbMultiPolygon)
newField = ogr.FieldDefn('DN', ogr.OFTReal)
dst_layer.CreateField(newField)
gdal.Polygonize(srcband, None, dst_layer, 0, [], callback=None)
dst_ds.SyncToDisk()
dst_ds=None

Unfavorable Result:


texts: GeoServer - How to display Ortho ortho (JP2 file) when georeference of each file are located in jgw file

texts:
I wouldlike to create a layer with all my tiles of ortho (jpeg2000). The problem here is the georeference of each jpeg2000 file are in .jgw file. When I try to add this on geoserver, I have always the same error:

Caused by: org.opengis.referencing.operation.OperationNotFoundException: Aucune transformation du système "EngineeringCRS[Wildcard 2D cartesian plane in metric unit]" vers "GeographicCRS[WGS84(DD)]" n'est disponible.


texts: QGIS: Symbology: Create Tercil for all values > 0 and tercil for all values< 0

texts:
I have a layer with multiple polygons each with a number of the commuter balance. You can see the histogram below (1st pic). I now want to symbolize this data.
I want a tercile calculation for all values < 0 and have 3 different red colors for it and a separate tercil calculation for all values > 0 and have 3 differenz green colors for it.
With the given algorithms it always includes 0 somewhere in the middle of the class breaks.
So preferably something like in the second image (the breaks here are just done by hand)
Is there any possibility beside split the layer apart and do two separate calculations?


texts:
I am deploying QGIS 3.6 on Windows 10.

Using the Raster > Projections > Warp (Reproject) tool, I need to reproject a single-band, float32, GeoTIFF DEM from EPSG:3413 to EPSG:3857.

The output raster must snap to a pre-existing GeoTIFF, also in EPSG:3857.

PROBLEM: The tool does not include a snapping option. The tool's help page is silent on snapping, as is a Google search.

There is a curious input section entitled "Additional Command Line Parameters".  But again, no guidance.

How to include snapping when reprojecting?


texts: How to create a network from nodes and boundaries?

texts:
I have a table of nodes (points) and i wish to route between the points. Its my understanding that i first need to build a pgrouting compatible topology (graph ?) that connects all the nodes using "source" and "target" fields.

I was just wondering what tools i can use to generate this topology ? assign_vertex_id doesn't seem to work.

I'm using postgis on a postgresql database.

edit: The nodes are basically points in free space (water), so there aren't any lanes as such. The only constraints are land mass boundaries.


texts:
I have a RGB raster layer (3 bands) in a File Geodatabase in ArcGIS. I export this layer to TIFF and then import that in QGIS. The problem is that QGIS shows a different (darker) view than ArcGIS. I checked raster style and all bands are correctly set to Red-Green-Blue.
What can i do in QGIS to solve this?

What i see in ArcGIS:
!()

what i see in QGIS:


texts: Colours for positive and negative values in a raster layer in QGIS

texts:
I've got a raster layer in QGIS 3.10. It's statistics are the following (hopefully understandable although in Spanish):
Valor mínimo: -0.36928805708885193
Valor máximo: 0.6040433645248413
Intervalo: 0.9733314216136932
Suma: 14433.653211801786
Valor medio: 0.12369652924773997
Desviación estándar: 0.09271679256133081
Suma de cuadrados: 1003.0713567335182

I would like to represent it with a Red-to-Black colour ramp, but centered to zero, so values < 0 appear in the red palette and those > 0 in the grey one. Logically, my attempts with an automatic classification always set the central white color into the positive values. Is it possible to do what I want and how?


texts: PyQGIS data defined property for label orientation

texts:
I have a QgsVectorLayer object with field Direction which determines orientation of the symbol of the layer.  I use setDataDefinedproperty method to adjust the orientation of the symbol like:
symbol = vLayer.renderer().symbol()
symbol.symbolLayer(0).setDataDefinedProperty(QgsSymbolLayer.PropertyAngle,QgsProperty.fromField("Direction"))

Now I want to add label for Name field with same orientation as symbols. I tried :
label_settings = QgsPalLayerSettings()
label_settings.fieldName = "Name"
label_settings.enabled = True

label_property_collection = QgsPropertyCollection()
label_property_collection.setProperty(QgsPalLayerSettings.TextOrientation,QgsProperty.fromField("Direction"))
label_settings.setDataDefinedProperties(label_property_collection)

labels = QgsVectorLayerSimpleLabeling(label_settings)
vLayer.setLabeling(labels)

But there is no change for label orientation. How can I achieve that?


texts: C# library for bounds based geometry search

texts:
I am searching for C# library which has code that builds spatial index in memory based on geometry bounds and enables fast searching for all geometry inside given bounds.

I prefer free library.


texts: Query a JOINED ArcGIS Online resultset (via an HTTP request)

texts:
In a related question, we determined that it's possible to query an ArcGIS Online FC via an HTTP request.
(I'll refer to the feature class as FC1.)
In addition to the above, I want to dynamically join from FC1 to a second feature class called FC2.

The join could be described as: Left join FC1 to FC2 on FC1.ID = FC2.ID.

Can such a join be performed in ArcGIS Online -- in a way that allows the resultset to be served up to an external system? (as a JSON object that is requested via HTTP)


texts: I am trying to create a transect profile and extract information for each specific bands of Landsat which will be further used to create a transect graph. Is there any direct way of achieving this using QGIS?


texts: How to convert a png with known bounds from Plate Carree to Mercator with GDAL

texts:
I have a png image that shows the region defined by the geo bounds UpperLat, LeftLon, BottomLat, RightLon (in degrees & -ve for West Lon and South Lat) in Plate Carree EPSG:4326. I want to convert it to Mercator EPSG:3857 using GDAL.

When I tried doing it this way:

Step1) gdal_translate -of Gtiff -a_srs EPSG:4326 pc_map_piece.png pc_map_piece.tiff

Step2) gdalwarp -s_srs EPSG:4326 -t_srs EPSG:3857 pc_map_piece.tiff mer_map_piece.tiff

I got this error:

ERROR 1: Unable to compute a transformation between pixel/line and georeferenced coordinates for pc_map_piece.tiff. There is no affine transformation and no GCPs.

Then I tried doing it this way:

Step1) gdal_translate -of Gtiff -a_ullr UpperLat LeftLon BottomLat RightLon -a_srs EPSG:4326 pc_map_piece.png pc_map_piece.tiff

Step2) gdalwarp -s_srs EPSG:4326 -t_srs EPSG:3857 pc_map_piece.tiff mer_map_piece.tiff

I didn't get the above error but the output image is not converted to mercator projection, I simply got the same image upside down as output.

Any ideas to get around this problem and convert the image to Mercator? Thanks.


texts: Where do I store Layer level information in QGIS?

texts:
I am using a python script to load the content of some files to layers in QGIS. These file usually contain header information that I am skipping (like the Number of header lines to discard parameter in the Delimited Text dialog box).
I would like however to preserve this information as I may have to output my data in the original format.
Where can I store these header ? Is Metadata the way to go ? If so, how do I store my information there and retrieve it later ?


texts: Openlayers 4.6.5 Drawing not populating on map

texts:
The code is as follows

@Input( 'coreMap' ) coreMap :CoreMapComponent
  vectorLayer : VectorLayer;
  modifyLayer : Modify;
  draw: Draw;
  vectorSource: VectorSource;
  style: Style;
  snap: Snap;

  ngOnInit(){
    this.vectorSource = new VectorSource({wrapX: false});
    this.style = new Style({
      geometry: 'circle',
      fill: new Fill({
        color: '#ffcc33'
      }),
      stroke: new Stroke({
        color: '#ffcc33',
        width: 2
      }),
      image: new CircleStyle({
        radius: 7,
        fill: new Fill({
          color: '#ffcc33'
        })
      })
    })
      this.vectorLayer = new VectorLayer({
        source: this.vectorSource,
        style: this.style
      });
    this.modifyLayer = new Modify({source: this.vectorSource});
    }

  addLayer(){
    this.coreMap.layerManager.addLayer(this.vectorLayer)
  }
  getLayers(){
    console.log(this.coreMap.map.getLayers().getArray());
  }

  addInteraction() {
    this.coreMap.map.addInteraction(this.modifyLayer);
    this.draw = new Draw({
      source: this.vectorSource,
      type: 'Circle'
    });
    this.coreMap.map.addInteraction(this.draw);
    this.snap = new Snap({source: this.vectorSource});
    this.coreMap.map.addInteraction(this.snap);
  }

  changeInteraction(){
    this.coreMap.map.removeInteraction(this.draw);
    this.coreMap.map.removeInteraction(this.snap);
    this.addInteraction();
  }

This will allow me to bring up the circle drawing tool and set the size ect but when I double click to end the drawing it doesn't save anything on the map, I thought maybe my new layer is being put at the bottom of the array but ol documentation says it should be at the top of the array.

*** Update

If I remove the styles from the new layer everything seems to work, why is this? As I am planning on making the styles customisable.


texts:
I'd like to give QField a try but I can not install it from Google Play Store. I do not want to create a Google account for privacy concerns.

Is there a possibility to install QField without using Play Store and a Google account, for example by using F-Droid?


texts: Getting the largest date value for each record

texts:
I have a layer in QGIS that has four date fields - date1, date2, date3, date4. Using pyqgis I wish to find the latest date (ignoring null values) for every record from these four columns and copy that value into a field called MaxDate.

If the four date fields are null, then MaxDate is set at the value of last_date

It feels like there should be an easy solution, but I have yet to crack it!

The fields are all 'QDate' format (yyyy-MM-dd)


texts: Automatically calculating geometry area when adding polygon in ArcMap?

texts: In ArcMap 10, I realize I can manually do this after the polygon feature is calculated by right-clicking the column in the attribute tables to calculate geometry or create a python script to do this but neither of these seem to allow me to automate this for new features without having to take those manual steps after the polygon is created and added to the table.  I want this to happen automatically once the feature is created and shows up in the attribute table, much like the auto-generated shape_area column works.  Maybe I can use those values to convert them to acres in a new column?


texts: Finding name of field that contains values that meet certain criteria using ArcPy?

texts:
I am trying to write a script that will loop through a table, and based on the values that are contained in the field of a table of a feature class in a geodatabase, it will print the name of the field that contains values that meet a certain crireria. Here is what I have so far:

fields = [field.name for field in fieldlist]
cursor = arcpy.da.SearchCursor(fc, fields)
 for row in cursor:
    for col in row:
       if isinstance(col, ( int, long ) ):
          pass
       elif type(col) == unicode:
           Strcol = unicodedata.normalize('NFKD', col).encode('ascii','ignore')
           if int(len(Strcol)) < 10 and Strcol.isdigit():
                if Strcol.startswith("55") or Strcol.startswith("56") or Strcol.startswith("57"):
                    print Strcol
                    print field.name

What this script does for me currently is it print out all of the values in all of my tables that match the criteria, but, with the field.name print statement, it automatically defaults to the last field in my field list. Is there a way to print the name of a field automatically when it finds values of the field that meet a desired criteria. (eg a field that contains values that are unicode strings, less than 10 characters, that are digits, which start with 55, 56, or 57)


texts: QGIS: RGB Image display results in odd colors (MODIS)

texts:
I am trying to make a rgb map of MODIS data (Dataset09 Radiances) with QGIS.

First I converted it to tif with gdal (and python). Since I did not really understand all of GDALs funcionality I used a lot of temp. results:

gdal_merge.py -v -o band1.tif HDF4_EOS:EOS_GRID:MOD09GA.A2001244.h20v11.005.2008275152329.hdf:MODIS_Grid_500m_2D:sur_refl_b01_1 HDF4_EOS:EOS_GRID:MOD09GA.A2001244.h20v12.005.2008275162824.hdf:MODIS_Grid_500m_2D:sur_refl_b01_1
gdal_merge.py -v -o band3.tif HDF4_EOS:EOS_GRID:MOD09GA.A2001244.h20v11.005.2008275152329.hdf:MODIS_Grid_500m_2D:sur_refl_b03_1 HDF4_EOS:EOS_GRID:MOD09GA.A2001244.h20v12.005.2008275162824.hdf:MODIS_Grid_500m_2D:sur_refl_b03_1
gdal_merge.py -v -o band4.tif HDF4_EOS:EOS_GRID:MOD09GA.A2001244.h20v11.005.2008275152329.hdf:MODIS_Grid_500m_2D:sur_refl_b04_1 HDF4_EOS:EOS_GRID:MOD09GA.A2001244.h20v12.005.2008275162824.hdf:MODIS_Grid_500m_2D:sur_refl_b04_1

gdal_merge.py -separate band1.tif band4.tif band3.tif -o MOD09GA.A2001244-allBands-noRepro.tif

I also did some reprojection:

gdalwarp -of GTIFF -t_srs EPSG:4326 MOD09GA.A2001244-allBands-noRepro.tif warpRepro.tif

If I load the result into QGIS every pixel has a different color, but similar values do not have similar colors. (The picture reminds me of the "ant-wars" or "blizzard" or "static noise" in old TVs, when there was no signal)
I googled a bit and read that QGIS can only display Rasters with values between 0 and 255. I tried scaling, similar to:

gdal_translate -scale -ot Byte band1.tif band1rescale.tif

If I load the new scaled and merged tif in QGIS, the result is a rgb image but the colors look too weak (brightness and contrast seems off. i think in a photosoftware I would solve the problem with an histogramm normalization (equalization???) )

Maybe somebody here knows a solution?

I also wrote my own normalization in python, but with the same result. The snow on the scene is super bright and all the other ground is dark and undistinguishable.


texts: Extract by location in QGIS using a Python script

texts:
I am using QGIS version: 3.10.3-A Coruña.
I have two shapefile layers, the first is polygons for world countries and the second is a 1-degree spaced polygon worldwide.
I would like to extract which polygon(s) from the second layer intersecting each country from the first layer.
I would do it easily with "Extract by Location", BUT my goal is to rename each output following the country's name and not to do it individually.

Here is the log file for one country when I do it individually.

Algorithm 'Extract by location' starting…
Input parameters:
{ 'INPUT' : 'C:/Users/Familie Hammoud/Desktop/ne_50m_admin_0_countries/World 1deg grid.shp', 'INTERSECT' : QgsProcessingFeatureSourceDefinition('ne_50m_admin_0_countries', True), 'OUTPUT' : 'C:/Users/Familie Hammoud/Desktop/Neuer Ordner/Germany.shp', 'PREDICATE' : [0] }

I am thinking if there is Python script to do this mission but sadly I don't know how.


texts: Earth engine reproject: why does reprojecting a pixel at 30m scale gives me pixels with an area of ~550?

texts:
Trying to understand how reproject() in Google Earth Engine works. To find an example, I am using ee.Image.pixelArea(), which gives an image containing the area of each pixel in square meters.

Now reprojecting this at the 30m scale, I expect that I will end up with pixels having an area of 900.

However, the data shows that my pixels have actually an area of ~557... what is hapenning? What is reproject() doing in that case?

var projection = ee.Projection('EPSG:3857').atScale(30)
var image_area = ee.Image.pixelArea()
var image_area_30 = ee.Image.pixelArea().reproject(projection)

Use the inspector to see the value of each individuals pixels, or:

Map.addLayer(image_area_30, {min:500, max:550}, "image_area_30")

EDIT/ANSWER

Bad question: the native projection of earth engine is a Mercator one, which does not preserve area. So indeed, need to go close to Equator, even if there seems to be around 893. Choosing instead an Albers equal area shows values for each pixel at 899.88 far away from Equator.


texts: ArcMap doesn't recognize geodatabase files under "Add Data"

texts:
I have a folder filled with geodatabase files of the following formats: ATX, FREELIST, GDBINDEXES, GDBTABLE, GDBTABLX and SPX. I would like to convert these to shapefiles using ArcMap. I've been looking around StackExchange for similar problems, but none are quite the same. For instance, see the solutions listed at Converting File Geodatabase feature classes to shapefiles in ArcGIS for Desktop?.

I tried to export the gdb to a shapefile after adding it as a layer. But, under "+ Add Data", the gdb files are not listed. I can't add anything as a layer to TOC unless it shows up in the "Add Data" window.

I tried the "Feature Class to Shapefile" tool, but for input files, the files are not listed here either. I also tried playing around with Catalog, but no success with the conversion.

Any idea how I can work something out of this?


texts:
I am trying to interpolate my soil data (metal concentration) using IDW in ArcMap10. My question is: is there a way to set the interpolation extent to include only the immediate area where the sampling points are located? I have a number of areas containing  sampling points (on one map) and I would like to interpolate them within their respective polygons on the same map. The distances between my separate areas are hundreds of meters and I would like the areas between them to remain blank (showing just my background map). Can I use a mask in IDW interpolation? If so, how do I go about it? I was unable to solve my problem so far.
I'm a beginner in GIS and any help would be greatly appreciated.


texts: What happened to the landsat add-ons i.landsat.dehaze and i.landsat.trim in grass7?  I keep seeing references to them but cannot find them in the repository.


texts: Creating coordinate points based on polygon

texts:
I am trying to create coordinate points based on a polygon that I've created in Google Maps and imported into R as an sf object. Specifically I have sample points of physical samples I've taken, but want the coordinates of every vine in the vineyard I sampled. The vines are spaced regularly in an 8ft * 8ft square over a rectangle with 26 rows and 59 columns. How can I create a dataframe of points?

Here's what I've done thus far to create regularly spaced points within the polygon extent.

i <- grep("California", epsg$note, ignore.case=TRUE)
   # first three
epsg[i[1:15], ]
# 717                +proj=lcc +lat_1=41.66666666666666 +lat_2=40       +lat_0=39.33333333333334 +lon_0=-122 +x_0=2000000.0001016 +y_0=500000.0001016001 +datum=NAD83 +units=us-ft +no_defs
# 718 +proj=lcc +lat_1=39.83333333333334 +lat_2=38.33333333333334 +lat_0=37.66666666666666 +lon_0=-122 +x_0=2000000.0001016 +y_0=500000.0001016001 +datum=NAD83 +units=us-ft +no_defs
# 719            +proj=lcc +lat_1=38.43333333333333 +lat_2=37.06666666666667 +lat_0=36.5 +lon_0=-120.5 +x_0=2000000.0001016 +y_0=500000.0001016001 +datum=NAD83 +units=us-ft +no_defs
library(raster)
library(rgdal)
f <- readOGR("dem_shape.kml")
f
#-122.8037, -122.8019, 38.6641, 38.66549  (xmin, xmax, ymin, ymax)
crs(f)
newcrs <- CRS("+proj=lcc +lat_1=39.83333333333334 +lat_2=38.33333333333334 +lat_0=37.66666666666666 +lon_0=-122 +x_0=2000000 +y_0=500000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs")
rob <- spTransform(f, newcrs)
plot(rob)
#
library(tidyverse)
library(sf)
RR<-st_as_sf(rob)
st_area(RR)/5.94725
#=9900.83 m^2
#=1664
bbox(f)
plot(rob) ; points(spsample(rob, n=1664.649, type='regular'), col='red', pch=3, cex=0.5)```

The issue is that, if each cross in the plotted image were a vine, they are not spaced in a square pattern in realtion to the polygon extent.


texts:
I'd like to get some expert opinions on how to get started in web mapping. First, a bit about my skills, then what I'd like to do.

I've got many years of experience with GIS, mapping, and cartography. And I'm pretty good with PHP, mySQL, HTML, and CSS. But I now need to bring all of that together into a web mapping project (I can't avoid it any longer!).

What I want to build is a mapping site where I can have something like a Google or Bing map base, then on top of that bring in my own township grid, and many thousands of points, with several pieces of data for each one. These points will represent wells or sampling points and will be updated monthly, and some of them will be updated daily. I'd even really like to have it so that users can add their own information to these points, and I'd really really like to have users able to draw in their own lines and polygons, and text tabels.

Is there a general type of approach to something like this?

I'm looking at Manifold [which I already use] and mapserving.com, or MapserverPro and one of the mapping packages they provide.

What about OpenLayers and embedding OL maps right into my own pages?


texts:
I'm struggling with deleting an mxd in a for loop. My code is:

    import arcpy, glob, os, datetime, time
    from time import strftime

    for filename3 in glob.glob(os.path.join(folderPath, "*.mxd")):

    fullpath = os.path.join(folderPath, filename3)
    basename, filename2 = os.path.split(fullpath)

    mxd = arcpy.mapping.MapDocument(filename3)

    if os.path.isfile(fullpath):
        basename, filename2 = os.path.split(fullpath)

        (shortname, extension) = os.path.splitext(filename2)

    mxd.saveACopy(draftloc + "\\" + "SentToDraftOn(" + strftime("%d%b%y %Hh%Mm%Ss") + ")"+ "_" + shortname + extension)
    arcpy.Delete_management(mxd)

I'm trying to do some processing on each mxd in a directory and then I want to saveACopy to a new location and then delete the old mxd. I keep getting the following traceback error:

  Object: Error in executing tool

I think it's an indentation thing, but I can't seem to figure it out. Can anyone suggest what I am missing?

Thanks,
Mike


texts: Filtering for slope using kernel size in ArcGIS Desktop?

texts: I have problem in smoothing the slope derived from digital elevation model. I want to do filter for the slope <5 degree in 7 X 7 kernel size, 5-20 degree in 5X5 kernel and >20 degree slope in 3X3 kernel size. But I am not able to put these kernel size condition in filter tools in ArcGIS Desktop 9.3.


texts: Is having rectangular cell sizes in a raster bad (after resampling)?

texts:
I used rasterio to scale a raster I have to a new resolution. It was at (5.0, 5.0) resolution, and I wanted to make it (15.0, 15.0).

However, the resulting resolution is not square anymore.

Here's my code (it's pretty standard, pretty much just from the rasterio documentation):

with rio.open(f_path) as ds:

        orig_res = ds.res[0] # 5.0
        des_res = 15.0

        # Calculate scale factor for this specific input layer
        upscale_factor = orig_res / des_res
        print("upscale factor: ", upscale_factor)

        # Scale while reading in array
        scaled_arr = ds.read(
            out_shape=(
                ds.count,
                int(ds.height * upscale_factor),
                int(ds.width * upscale_factor)
            ),
            resampling=Resampling.bilinear
        )

        # Also have to update the transform to be scaled
        scaled_transform = ds.transform * ds.transform.scale(
            (ds.height / scaled_arr.shape[1]),
            (ds.width / scaled_arr.shape[2])
        )

        scaled_profile = ds.profile
        scaled_profile.update(transform=scaled_transform, height=scaled_arr.shape[1], width=scaled_arr.shape[2])

        with rio.open(f_out_path, 'w', **scaled_profile) as scaled_ds:
            scaled_ds.write(scaled_arr)

            print(f"Orig arr shape: ({ds.count}, {ds.shape[0]}, {ds.shape[1]})")

            print(f"Scaled ds shape: ({scaled_ds.count}, {scaled_ds.shape[0]}, {scaled_ds.shape[1]})")

            print("Orig ds res: ", ds.res)
            print("Scaled ds res: ", scaled_ds.res)

Here is the output printed:

upscale factor:  0.3333333333333333
Orig arr shape: (1, 9240, 2671)
Scaled ds shape: (1, 3080, 890)
Orig ds res:  (5.0, 5.0)
Scaled ds res:  (15.0, 15.00561797752809)

So is it bad that my new raster does not have square pixels anymore? (15.0, 15.00561797752809)

Will this affect anything with how I use the raster, or will it be okay? I also have other rasters that I'm up/down sampling to all be (15.0, 15.0), and they are also not exactly 15.0, so will it cause issues if my rasters have very slightly different resolutions?

Or is there a way to remedy this?


texts: Ordering of Band Import into QGIS Raster Calculator?

texts:
I am using QGIS 3.4,3

I want to use Raster Calculator to create NDVIs using Planetscope AnalyticMS_SR.tif

I know that the band order for Planetscope is BGRN
When I import the multiband SR into Raster calculator I see the bands numbered 1 to 4

When I look at the multiband SR in Layer Properties/Histogram the bands are numbered and colored indicating that B1 is R, B2 is G, B3 is B and B4 is purple and must be NIR. This is not the order of bands that Planetscope provides

When I look at the SR file in S-A CP Band Set I see a list of bands 1-4, but I cannot find other information about them because Planetscope is not included in the list for quick wavelength settings

How can I be sure that I have Red and NIR bands correctly identified in Raster Calculator to create the NDVI?


texts: Convertind a sketch up model to GOCAD 3ds .ml file type

texts:
I use a program called GOCAD. It exports the models created in .ml file format. I was wondering if there is a way to take a simple terrain block created in Sketchup and convert it to .ml?
Any help or insight?


texts: How to calculate the scale of current Zoom Level with parameters like extent, resolution?


texts: Sync Enabled Feature Service at ArcGIS for Server Version 10.1

texts: Is sync enabled Feature Services available for ArcGIS server at version 10.1?  If not at what version was this introduced.  I'm having trouble finding a definitive answer on the web.


texts: How to enable pop-up window for more than 1000 features in ArcGIS Online

texts:
Is there any way to enable pop-up window in ArcGIS Online for more than 1000 features for layers' information?

When I try to add a shape file in arcgis.com to the map I get the message saying:

  The maximum  features allowed are 1000.


texts:
I'm new to ArcGIS/Esri.  I have a working sandbox example using SketchViewModel, but need to change that to:

I want to customize the drag and drop/resizing points on the surrounded graphic
I want to be able to do a 'snap to' on a resize of a graphic in one dimension (Y for example), I only use rectangles.
Ideally I'd like to branch in Git off of SketchViewModel.  Alas the un-minified source is not evident in Git, leaving me to recode SVM.


texts:
I'm trying to convert the census tract data with demographic profiles to a raster so that I can do a site suitability analysis with proximity to high youth population areas.

I downloaded the 2010 Census Demographic Profile 1 Census Tract shapefile from https://www.census.gov/geo/maps-data/data/tiger-data.html and used select by location to extract all polygons within the county boundary.

I then added a field and calculated the combined ages of population under 19. I want to use this field to convert to raster. I've tried polygon to raster, feature to raster, I've tried merging the data first and each time the file comes up blank. ArcMap will create the file but it will not contain any data. I'm not sure what I'm doing wrong or if I'm missing a step here.


texts: Getting ST_ClosestPoint for Geography type?

texts:
Using PostGIS 2.0.0 I'd like to find the point on a LINESTRING that is closest to a given point. The LINESTRING represents a great circle line (ie. geography type). ST_ClosestPoint appears to do exactly what I want, however, I find that the returned point does not lie on the great circle line:

WITH points AS (
    SELECT ST_GeomFromEWKT('SRID=4326;POINT(41 12)') AS point,
        ST_GeomFromEWKT('SRID=4326;LINESTRING(40 5,41 15)') AS border
)
SELECT ST_AsText(ST_ClosestPoint(border, point)), -- POINT(40.7029702970297 12.029702970297)
    ST_DWithin(ST_ClosestPoint(border, point)::geography, border::geography, 700) -- false
FROM points;

It's actually over 700 metres from the line. I believe this happens because ST_ClosestPoint operates on a plane. How would I do the same thing on the spheroid (or at least on the sphere)?

I tried @MerseyViking's function (ported to plpgsql) and it mostly works very well, but in some cases returns a point that is not on the line. A simple example:

WITH points AS (
    SELECT ST_GeomFromEWKT('SRID=4326;POINT(10 3)') AS point,
    ST_GeomFromEWKT('SRID=4326;LINESTRING(10 5,10 15)') AS border
),
calc AS
(
    SELECT ST_ClosestPoint(border, point)::geography AS plane_closest,
        closest_point_to_line(border, point) AS sphere_closest
    FROM points
)
SELECT ST_AsText(plane_closest) AS plane_closest, -- POINT(10 5)
    ST_DWithin(plane_closest, border::geography, 1), -- true
    ST_AsText(sphere_closest) AS sphere_closest, -- POINT(10 3)
    ST_DWithin(sphere_closest, border::geography, 1), -- false
FROM points, calc;


texts: Extract specific amenity Osmosis with Java api

texts:
I need to extract restaurants from a planet.osm.pbf file. I assume the best way to do this is with Osmosis, however I need to save this data in my own database format. (Not something standard like Postgres or something else with existing plugins). Therefore, I need to write a Java application to do this, yet the Osmosis documentation seems very confusing to me. What I need to do is this:

Pseudocode:

while(file.hasNext("amenity=restaurant"))
{
     Amenity restaurant = this.amenity;
     String street = restaurant.street;
     String city = restaurant.city;
     String region = restaurant.region;
     String name = restaurant.name;
     String cuisine = restaurant.cuisine;
     // I do my own thing with this information here...
}

I can read the file in with FileInputStream, etc... I just don't understand how to parse and extract data specifically.

This information would be extremely valuable, allowing people to extract this type of data with ease, so they don't have to scour the internet as I did.


texts: Is it safe to store coordinates on NFC tags always in one CRS and rely on conversion for others

texts:
I'm involved in NFC (Near Field Communication) standardization of a record structure for geographic point locations. That will allow to put coordinates on NFC memory tags which can be read with a smartphone or to exchange between phones and the like.

A first proposal made was to always include the coordinate tuple (2D or 3D) and the CRS, so a coordinate tuple for a specific CRS could be stored. To this, an alternative opinion claimed that the CRS could set implicitly (for example WGS-84) and coordinates for any other CRS be computed, of course only for the same dimensions.

It may be worth to note that a compact encoding format is important because NFC tags, especially the cheaper ones, do not provide a lot of memory. So any byte less is good if functionality isn't reduced (one envisioned use case is that NFC tags could be affixed to landmarks and allow the location be read electronically - this was thought to sometimes require other CRS than WGS-84).

So the question is, if a coordinate tuple is given in one CRS is it safe to assume that it can be converted to any other CRS of practical relevance?

To make this more open-ended, would there be any drawbacks if a data format for geographical coordinates is fixed to only one CRS and requires conversion to and from other CRSs if needed.

And as further explanation, the intention is not to derive a format for a specific application but for general use by anyone who wants to store or read a location from an NFC tag/card or share with another smartphone when touched (like to show the location on a map application). Similar to what RFC 5870 (Geo URI) provides but more compact.


texts: Error in reading point cloud data (.las) in R?

texts:
I want to generate a DEM using a point cloud (.las) for which I tried reading it in R using rlas, lidR and rLiDAR : readLAS command and none of it seems to work.

The LASfile is not being read at all. The error is :

Error in readBin(con, "raw", size = 1, n = 4, endian = "little") :
  can only read from a binary connection
In addition: Warning message:
In file(LASfile, open = "rb") :
  file("") only supports open = "w+" and open = "w+b": using the former

And the code is:

library(lidR)
library(rLiDAR)
library(rlas)
setwd("D:/Thing/Path")
LASfile <- system.file("extdata", "ABCD.las", package="rLiDAR")
rLAS <- readLAS(LASfile, short = TRUE)


texts: Cannot draw straight line polygon in QGIS 3.16

texts: I tried to draw polygon by using add polygon. It should be a straight polyline and then become a polygon but it turns to a curve line like arc (circular line) so that I couldn't draw a polygon. It's very strange. I download QGIS 3.18 and the result is the same. I don't know how to solve this problem.


texts:
I have a raster that I want to reclassify and I'm using GRASS r.recode from within QGIS v.2.4.0. The rules file I am using has this content:

0:149:1
150:229:2
230:379:3
380:459:4
460:*:5

The raster I have has a maximum value of 515 and according to the rules file, it should be reclassified as 5 but the range in the output raster from r.recode is [1,4].

Any ideas as to what is wrong with my rules file?


texts:
I have a form / survey created for Survey123. It has been published and I can access the data which has been captured through the ArcGIS Online Portal.

I would like to read the attribute information from this form in Excel. Similarly as I have connections to my SQL Databases and I can read the attribute information stored in the databases connected to my Collector for ArcGIS maps.

I know how to download the CSV / SHP / FGDB information, but I would like to connect to this information without the need to have a user logon to the AGOL and download it. Some users will not have ArcGIS Installed on their machines, so using Python is an option, but I cannot use arcpy.

Is there a workflow that will allow me to have a live connection between Excel and the attribute information stored in my Survey123 Form?


texts:
I am struggling to understand the projections using pyproj.

for now my question is to understand the results of projection operations

I have following coordinates that I project on x,y

lon  = [3.383789, 5.822754]
lat = [48.920575, 53.72185]

# with Basemap
from mpl_toolkits.basemap import Basemap
M = Basemap(projection='merc',ellps = 'WGS84')
q1, q2 = M(lon, lat)
print zip([x for x in q1], [x for x in q2])
# this gives me
#[(20414190.011221122, 65799915.85232639), (20685694.35308374, 66653928.94762396)]

# with pyproj
from pyproj import Proj
p = Proj(proj='merc', ellps='WGS84',errcheck = True)
print p(3.383789, 48.920575), p(5.822754, 53.72185)
#this gives me
#(376681.6684318804, 6229168.979819128) (648186.0102944968, 7083182.075116195)

why are the results different while I use the same projections parameters

as a newbie in geospatial data processing I apologize in advance for a question which may be trivial

P.S: I edited the questions following comments


texts: Extract a feature from a shapefile layer, rasterize it and save to a numpy array

texts:
I have the following problem:

The following sequence of operations:

loading a shp file using a gdal package,
extracting a layer out of it (every feature in this layer is a polygon),

leads me to a place where rasterizing every feature in its own extent to a separate numpy.array seems to be an impossible task.

Is there any way to achieve that without hacking a osg package in a really cumbersome way?


texts: Changing default projection in ArcPad 10

texts:
Is there a way to set a default projection in ArcPad 10?
The only way I know is by creating a quickProject based on a project that contains a shapefile with the same projection.


texts: How do the Identify layers Project settings work?

texts:
In QGIS 2.14 project settings one tab says "Identify Layers". Here you can select or unselect Layers. From my understanding this should define if a layer reacts if it´s clicked with the identify tool. But even if I uncheck all layers I can still use the tool and the Object Information Dialog is shown.

Do I missunderstand this setting?
What´s the purpose of it?


texts: Extract one class from global LULC in google earth engine

texts:
I am trying to extract one LULC class (e.g. cropland) for the entire globe from the MODIS global lulc. How it can be done?
Code I am using to display global lulc is:

    var cover = ee.Image('MODIS/051/MCD12Q1/2012_01_01').select('Land_Cover_Type_1');
    var cover_palette = [
    'aec3d4', // water
    '152106', '225129', '369b47', '30eb5b', '387242', // forest
    '6a2325', 'c3aa69', 'b76031', 'd9903d', '91af40',  // shrub, grass
    '111149', // wetlands
    'cdb33b', // croplands
    'cc0013', // urban
    '33280d', // crop mosaic
    'd7cdcc', // snow and ice
    'f7e084', // barren
    '6f6f6f'  // tundra
     ];

    Map.addLayer(cover, {min: 0, max: 17, palette: cover_palette },
                'Land Cover classification');


texts: Changing data (attributes) in .csv of data vector layer and QGIS automatically adopting this?

texts:
I have a .csv with data for 300 points: Each column on the .csv is a particular object and the cells define the amount that is available in each point.

If I need to change this amount, is it possible to change this in the .csv file in order that QGIS takes these edits automatically into its attribute table?


texts:
I think I am a bit confused with CartoDB when doing maps at runtime.

Let's say I have a Cluster Map that includes all my data. However, I want to filter by a field. I thought one could simply set the new SQL (SELECT * FROM mytable WHERE field = 'foo') to the sublayer,as it is done inside the platform, to select a subset and the visualization was to to be applied.

This does not appear to be the case, the SQL is a complex expression that includes all the calculation and has a lot of parameters making very difficult to reuse it.

Does anyone have a solution for this case? I dont think you MUST create a new visualization.

regards,


texts: Print pixel band value for a point to the console in Earth Engine

texts:
I want to display the pixel value of Landsat B4 reflectance for a specific image date and coordinate. The following code shows how I select the pixel of interest. How do I extract the value and display it to the console?

// The pixel's coordinate
var lng = -91.2363;
var lat = 30.0919;
var r = ee.Geometry.Point(lng, lat); // r is the pixel of interest

// Load Landsat 7 imagery and filter it
var collection = ee.ImageCollection('LANDSAT/LE07/C01/T1_SR')
  .filterDate('2004-05-14','2004-05-15') // the date, on 14th May, 2004
  .filterBounds(r);


texts: Pattern did not match any band error for Sentinel 2, 2015 year image in Google Earth Engine

texts:
I am trying to calculate MNDWI index for Sentinel 2 image of an area for 2015 year data in Google Earth Engine.
var image2 = ee.ImageCollection('COPERNICUS/S2')
                  .filterDate('2015-06-23', '2015-12-31')
                  // Pre-filtering to get less cloudy granules.
                  .filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', 20))
                  .filterBounds(point)
                  .select(['B3', 'B11'])
                  .reduce(ee.Reducer.percentile([20]))
                  .rename(['B3', 'B11'])
var mndwi2 = image2.expression('(green - swir) / (green + swir)' ,{
  'green':image2.select('B3'),
  'swir':image2.select('B11'),
}).rename('MNDWI2');
//adding the layer
Map.addLayer(mndwi2.clip(roi), imageVisParam, 'MNDWI_2015');

//making water zones from 2020 mndwi image
var zones_2015 = mndwi2.gt(0.1).unmask(0);
Map.addLayer(zones_2015.clip(roi), VIS_WATER_MASK, 'water zones for 2015');

It is giving error :

MNDWI_2015: Layer error: Image.select: Pattern 'B11' did not match any
bands. water zones for 2015: Layer error: Image.select: Pattern 'B3'
did not match any bands.

But as I am changing my year from 2015 to any other year it is working fine. Is there some issue with particular 2015 data or something wrong in the code? Some days back I apply the entire procedure for the same 2015 year for other area and everything works well but this time it is giving error for trying for other area.
Script link: https://code.earthengine.google.com/47dd2ed2fc7ca17b90c769e1c91ac971


texts:
Should be really straight forward, but I can't understand what I'm doing wrong.  I'm trying to measure the distance between the centre of a polygon (suburb) and an arbitrary point.
SELECT
    ST_AsText(ST_Centroid(geom)) as poly_centroid
    , ST_AsText(ST_SetSrid(ST_MakePoint(151.182995, -33.890395), 4283)) as pin
    , ST_Distance(
        ST_Transform(ST_SetSrid(ST_MakePoint(151.182995, -33.890395), 4283), 3112),
        ST_Transform(ST_Centroid(geom), 3112)
    ) as dist
FROM public.suburbs

I'm using 3112, Australia Lambert projection (which is in metres) - and entering the arbitrary pin point in long/lat (EPSG 4283) - but the resulting distance is clearly wrong:

poly_centroid
pin
dist

POINT(150.9300421895883 -33.72537351810487)
POINT(151.182995 -33.890395)
4229563.228668396

POINT(151.03083602374133 -33.75769279240827)
POINT(151.182995 -33.890395)
4229563.161116713

POINT(151.15847466229258 -33.75004174996872)
POINT(151.182995 -33.890395)
4229563.120652905

POINT(150.93371519943338 -33.73946132331832)
POINT(151.182995 -33.890395)
4229563.214226549

POINT(150.97498052582384 -34.12744863603382)
POINT(151.182995 -33.890395)
4229562.838806998

These should be 10s of Km apart, not 1000's.
Note that the polygons are all loaded in EPSG 3112 projection, so the ST_Transform isn't necessary.


texts: Getting Surface Area on 3D Objects (ArcGIS Pro)

texts:
TLDR: I have a 3D multipatch feature class and need the surface area of the facade of the building.

I brought in some data from Revit into ArcGIS Pro, processed it (project it, georef), converted the data into a GDB feature class (Layer 3D to Feature Class). I used the Add Z Information tool to add the surface area of each 3D object, which is close to what I want, but not quite. I also tried field calculating geometry, using area and sq. ft as the unit type, but it only accurately captures geometries that are horizontal.

As (hopefully) shown in the screenshot, the surface area will calculate the facade plus the interior surface area, and the surface area of the 3D object wall thickness. Does anyone know of a way to just get the facade area? I've been banging my head against the wall on this one.

I have also seen this answer in StackExchange but I think this will just give me the same result as the surface area from Add Z Information.
 Calculating 3d polygons' surface area using ArcGIS Desktop?


texts: Finding least cloudy Sentinel-2 granule for non-rectangular ROI and given season in Google Earth Engine

texts:
We need to download cloudless Sentinel-2 data for a rather large number of regions of interest (ROIs), which are provided in the form of geojson polygons. Unfortunately, relying on the cloud cover metainfo tag will not help us, as the clouds can appear in a corner of the granule, which does not intersect with our ROI - we might up rejecting images that would be good for our ROI or keeping images that are good outside, but bad inside the ROI.

Thus, I thought, it could be a good idea to use Google Earth Engine (GEE) to search for granules, which are not affected in our ROI, exploiting the cloud band. Then, I would export a list of image IDs for subsequent downloading outside GEE.

In the end, the data is subject to a land cover classification task, and we want to have one (mosaicked, if necessary) image per season per ROI.

Any idea how I could get this list of image IDs for a bunch of ROIs? I am not looking for a full-fledged processing chain but rather for hints towards how to formulate the data query. The challenge is the required ROI-based point of view in contrast to the many answers with a pixel-based or granule/image-based point of view...


texts: How can the amount of columns/rows in a raster be altered?

texts:
I want to work with 5 raster sets in R. 3 of them are results of an Euclidian distance analysis (R1-3), 1 is from a kernel density analysis (R4) and the last one is from a viewshed analysis (R5). My goal is to stack them in R. To do this, all of the raster sets need to have the exact same extent (incl. CRC, cell size and columns/rows).
R1-4 are not the problem, they are properly stackable with the following extent/resolution:
X        10
Y        10
Columns  2826
Rows     2911
xmax     5661589.216435
xmin     5632486.563462
ymax     32371381.150053
ymin     32343124.735493
CRC      ETRS89/UTM32N

R5 has the following extent/resolution:
X        10
Y        10
Columns  2827
Rows     2911
xmax     5661589.216435
xmin     5632486.563462
ymax     32371381.150053
ymin     32343124.735493
CRC      ETRS89/UTM32N

So the whole difference between R1-4 and R5 is one single column. I don't really understand this because the cell size and the extent are the same. As far as I understand, the number of columns cannot logically be anything else than 2826.
What I now did was to export R5 in ArcGIS with a number of columns of 2826 which obviously changed the cell size to 9.998... each. I then tried to resample R5 to a cell size of 10. However, the number of columns remained unchanged in the newly exported raster.
My question would now be how the number of columns can be changed or which other ways there are to align the raster extents.


texts: In ArcMap, what is the difference between Copy Parallel and Buffer?

texts:
To build an offset of a linear feature in ArcMap, what is the difference between Copy Parallel and Buffer?  In terms of the resulting features and their linear differences from the source linear feature?

I am building an easement of a centerline feature.  I want to build a feature that is 25feet either side of my linear feature.  Where the centerline is going without direction changes there seems to be no difference.  Where the centerline changes direction the differences start to appear.

Thanks.


texts: Clustering Point for StatePlane WKID 2236

texts:
Does anyone have an example to point me in the right direction to utilize Clustering with the ESRI JS API working with a State Plane wkid: 2236? I want to provide a solution for an internal application just like the demo on ESRI JS website. Though, all our data and the application I'm supposed to update is defined wkid:2236. I've tried modding the demo code though either I'm overlooking something or it's not possible. If this isn't possible any suggestions?

I had overlooked the "spatialReference" option within the constructor. This gets me to the point of displaying the cluster symbols on the map. Though, when I click the symbol I received the following error and I included the call stack. Any ideas?

CHANGES TO ESRI POINT CLUSTERING EXAMPLE

var initExtent = new esri.geometry.Extent({ "xmin": **432163.89125**, "ymin": **1451074.4432**, "xmax": **710435.747749999**, "ymax": **1626999.3528**, "spatialReference": { "wkid": **2236** } });

...

var basemap = new esri.layers.ArcGISTiledMapServiceLayer("**http://gisserver/ArcGIS/rest/services/BaseMap/MapServer/**");

...

var photos = esri.request({

                    "url": "**js/testdata.json**",

                    "handleAs": "json"

                });

...

 clusterLayer = new extras.ClusterLayer({

                "data": photoInfo.data,

                "distance": 100,

                "id": "clusters",

                "labelColor": "#fff",

                "labelOffset": 10,

                **"spatialReference":"2236"**,

                "resolution": map.extent.getWidth() / map.width,

                "singleColor": "#888",

                "singleTemplate": popupTemplate

            });

ERROR
Unhandled exception at line 34, column 167548 in http://serverapi.arcgisonline.com/jsapi/arcgis/?v=3.1

0x800a01b6 - Microsoft JScript runtime error: Object doesn't support property or method 'toJson'

CALL STACK

  toJson [arcgis] Line 34 Script
      updateHighlight [PopupBase.js] Line 19  Script
      onSelectionChange [Popup.js] Line 19    Script
      select [PopupBase.js] Line 19   Script
      _updateFeatures [PopupBase.js] Line 19  Script
      setFeatures [PopupBase.js] Line 19  Script
      _2b5 [arcgis] Line 15   Script
      setFeatures [Popup.js] Line 19  Script
      _addSingles [eval code] Line 275    Script
      onClick [eval code] Line 126    Script
      advice [arcgis] Line 15 Script
      _504 [arcgis] Line 15   Script
      _onClickHandler [arcgis] Line 34    Script
      Anonymous Function [arcgis] Line 15 Script


texts: Mean reducer changes scale in Google Earth Engine?

texts:
I am working with NOAA CFSV2 data trying to get an Image that is a one-year average of an ImageCollection. However, in the process the scale of an image changes:

var uv = ee.ImageCollection('NOAA/CFSV2/FOR6H')
  .select('Downward_Long-Wave_Radp_Flux_surface_6_Hour_Average')
  .filterDate('2000-01-01', '2000-12-31');

// reduce the image collection to one image by summing all rasters within a given year:
var muv = uv.reduce(ee.Reducer.sum());

//Checking the scale:
var scale_uv = uv.first().projection().nominalScale();
print('Scale of uv', scale_uv);
//Scale of uv: 34771.29356168821
var scale_muv = muv.projection().nominalScale();
print('Scale of muv', scale_muv);
//Scale of muv: 111319.49079327357

Does anyone know why this happens and what is the way to avoid that?


texts: Measure length of line to any point along line

texts:
I am trying to create linear referencing to point. I have a line to which I have snapped to points to (without breaking it).
Now I want to find the distance from the start point of my lines to points snapped to it, this I struggle with.

Point 1. 45m from start of line
Point 2. 126 from start of line
etc.

Examples I have seen seem to focus on calculate distance between vertices but I just want to calculate length from start of line to point, for every point.

Does someone have an idea how this could be done?


texts: ArcGIS Flex application not accessing a server - no errors

texts:
I'm working on an application that has to access data from a server. The flex application uses Esri's World_Imagery basemap and is able to access ESRI_Census_USA from the sample server.

Here are the two active lines of my code (almost directly copy and pasted from their tutorial.

<esri:ArcGISTiledMapServiceLayer url="http://server.arcgisonline.com/ArcGIS/rest/services/World_Imagery/MapServer" />
<esri:FeatureLayer graphicAdd="featurelayer1_graphicAddHandler(event)" outFields="[*]" url="http://sampleserver1.arcgisonline.com/ArcGIS/rest/services/Demographics/ESRI_Census_USA/MapServer/0" />

This works just fine as long as I want to show census data, but if I try to access the ArcGIS Server software running on my machine (the same machine), nothing appears. To do this I change the code in the second line to this:

https://localhost:6443/arcgis/rest/services/Pennwell/MapServer/0

Even the SampleWorldCities layer on the server won't load, so it's either the Flex application not accessing the server properly or it's the server itself not being set up right (there might be something else going on, too, but these are the most likely as far as I can tell).

Screenshot of the ArcGIS for Server Manager

I am able to get the SampleWorldCities (and my other data) to open in the online viewer without any problems.

How do I go about getting my Flex application to talk to the server?


texts: Function Evaluator doesnt work in network analyst

texts:
I made a function in my evaluator (% of asphalt in my route), but for some reason, it doesnt work.
The idea of the script is to se the % of asphalt in my route, so basically the equation is

(km of asphalt)/(km of my total route)

It is a really simple equation, but when I try to insert this evaluator, the network does not solve, and a message like this appears:

  Error: The evaluator failed to return a value. [Attribute: Distance,
  Source: TF_4589, OID = 42, EID = 46897]

Anyone know how to solve that?


texts: QGIS: "Invalid geometry" even after using validity checker and v.clean

texts:
I'm trying to intersect two layers in QGIS, similar to this unanswered question. I've performed Vector > Geometry Tools > Check Geometry Validity and I've run GRASS's v.clean, per these questions.

I'm still getting this response:

GEOS geoprocessing error: One or more input features have invalid geometry.

Is there a better way to check geometry?
Is there a way QGIS can tell me where the problem in the geometry is?
Can anyone see what I'm doing wrong using my files?

Here is a link to the underlying set of Voronoi polygons I want to intersect:
https://www.dropbox.com/s/edev53n7nn7b6wy/Voronoi1.zip?dl=0

Here is a link to the problematic outline of Liberia that's giving me an "invalid geometry" error when I try intersecting it with the above (already run through Check Validity and v.clean):
https://www.dropbox.com/s/4oa7ln7rk2n6e23/Liberia%20Outline%20Valid%20Cleaned.zip?dl=0


texts: Angle between the magnetic field vector and the vector joining a satellite to a receiver station

texts:
I was wondering if someone could help me in finding the angle between the  magnetic field vector at a location and the vector joining a satellite to a receiver station. This angle is also known as aspect angle.
I know that we could find the angle between two vectors in 3-D space using the dot product of the two vectors but in order to employ that idea I must have the vectors in the first place.
How can I find the vector notations for the magnetic field vector at a given location and the vector notation for the line joining a satellite to a receiver station?


texts: Create GeoTIFF from raster image and shapefile with extent

texts:
I have a couple of hundred raster images (TIFF) with a specific file name, like 242_13_1_0.tif. I also have a shapefile with polygons that corresponds to the corners of one of these files. All polygons have an attribute id that could be used to link a specific polygon to a specific file.

I have previously, when I have done something similar, exported the shapefile to some text format and then created a Python script and with Python/GDAL mapped the image to the correct extent and then create a GeoTIFF. But it's a bit time consuming to have to write new specific scripts whenever this problem occurs, so my question is if this could be done more automatically? Like in QGIS?


texts: Ipyleaflet map object doesn't display in Jupyter Notebook but it gets created

texts:
When running the code from the example (https://ipyleaflet.readthedocs.io/en/latest/api_reference/map.html#example), no map get displayed.

from ipyleaflet import Map, basemaps, basemap_to_tiles

m = Map(
    layers=(basemap_to_tiles(basemaps.NASAGIBS.ModisTerraTrueColorCR, "2017-04-08"), ),
    center=(52.204793, 360.121558),
    zoom=4
)

m

No error message pops up.
The map object gets created:

type(m)
>>> ipyleaflet.leaflet.Map

It just doesn't get rendered.
Any hints?

After a while, this message appears :
"A Jupyter widget could not
be displayed because the widget state could not be found. This could
happen if the kernel storing the widget is no longer available, or if
the widget state was not saved in the notebook. You may be able to
create the widget by running the appropriate cells."

It suggests that something with the widget is not working, but I can't understand what.

Versions:
ipyleaflet: 0.9.1
ipywidgets: 7.4.2


texts: Mapbasic Error - Alter Object cannot operate on an unitialized object variable

texts:
I get this error when running an mbx in Mapbasic

   Include "MapBasic.def"
   Include "Menu.def"

'assign brush/pen symbology for all but easements
        dim b1 as Brush
        b1 = makebrush (1,16777136,-1)  'Black colour, no fill
    dim p1 as Pen
        p1 = makepen (1,2,0)            'Black line
    dim x as Integer
    dim oUpdateObj as Object

    'Change Fill and Pen
    Open Table StrPathOt+"Easements.TAB" Interactive

    'Easements
    Select * from Easements where Str$(obj) = "Region" into NewStyle
        Fetch first from NewStyle
    x = 1
    Do
        oUpdateObj = NewStyle.obj
            Alter Object oUpdateObj Info OBJ_INFO_BRUSH, b1
            Alter Object oUpdateObj Info OBJ_INFO_PEN, p1

            Update NewStyle set Obj = oUpdateObj where RowID = x

            x = x + 1

            Fetch next from NewStyle

    Loop while x <= TableInfo(NewStyle,TAB_INFO_NROWS)
    Commit Table Easements

The issue reported is on this line

    Alter Object oUpdateObj Info OBJ_INFO_BRUSH, b1


texts: Load layer From Oracle using QgsDataSourceUri in QGIS3 always return Invalid Layer

texts:
I'm tring to load spatial layer from Oracle but when i called function QgsVectorLayer(uri.uri(), "LAYER" , 'oracle') it still inValid layer

This is my code

> def LoadFromOracle(tablename):
>     GeometryColumn = Get_GeometryColumn(tablename)
>     print("GeoMetry Column is")
>     print(GeometryColumn)
>     oracleSetting = settingdb.GetSettingOracle()
>     uri = QgsDataSourceUri()
>     uri.setConnection(oracleSetting[3], oracleSetting[4], oracleSetting[0], oracleSetting[1], oracleSetting[2])
>     uri.setWkbType(QgsWkbTypes.Polygon)
>     uri.setDataSource('', tablename, GeometryColumn)
>     uri.setKeyColumn('ID1')
>     vlayer = QgsVectorLayer(uri.uri(), tablename, 'oracle')
>     print(uri.uri())
>     print ("uri wkbtype == %s layer geometry type == %s wkbType == %s"%(uri.wkbType(),vlayer.geometryType(),vlayer.wkbType()))
>     crs = QgsCoordinateReferenceSystem("EPSG:24048")
>     vlayer.setCrs(crs)
>
>     return vlayer

I'm tring to use uri.setWkbType function but it is not working.
When I try to check valid status of QgsVectorLayer with this code.it still "False"

layer = self.LoadFromOracle('PROVINCE')
print ('Layer is valid = %s'%self.layer.isValid())

result of the code :

> dbname='XE' host=127.0.0.1 port=1521 user='SYSTEM' password='SYSTEM' key='ID1' estimatedmetadata=true type=Polygon table="PROVINCE" (GEOMETRY) sql=
> uri wkbtype == 3 layer geometry type == 3 wkbType == 0
> Layer is valid = False


texts: QGIS categorized styling: show number of objects for a categroy

texts:
When I have a vector layer with several objects and apply a categorized styling to the layer, ArcGIS Shows me the number ob objects that fall into each category. The following image is from QGIS but in ArcGIS it looks pretty much the same. In ArcGIS next to the legend is for each row a number with objects that fall into the category.

Is there a similar option in QGIS I'm not seeing? Is there an easy way to query the numbers?


texts: How to remove black pixels in Landsat image

texts: Does anyone know how I can remove these black pixels around the edge of this Landsat image? They appeared when I layer-stacked the bands in ERDAS IMAGINE - they were not there in the separate band images.


texts:
On the line after the field info I continue to get syntax errors and I'm unable to solve this problem myself, please help me.
The purpose of the code is to move a given point within a set distance in order to protect privacy of said point

# Step 11 - SELECT By Attributes -> where FID == FID_1, export to new Feature Class
# ------------------------------------------------------------------------------------------------------------------
# Make a layer from the feature class
arcpy.MakeFeatureLayer_management(in_features="RDiff_AreaIntersect_Dissolve", out_layer="lyr_eaten",
where_clause="",
workspace="",
field_info= "OBJECTID OBJECTID VISIBLE NONE", "Shape Shape VISIBLE NONE",     "ORIG_FID ORIG_FID VISIBLE NONE", "FIRST_Join_Count FIRST_Join_Count VISIBLE NONE"
FIRST_OBJECTID_1 FIRST_OBJECTID_1 VISIBLE None,
FIRST_gps_latitu FIRST_gps_latitu VISIBLE None,
FIRST_gps_longit FIRST_gps_longit VISIBLE None,
FIRST_OBJECTID FIRST_OBJECTID VISIBLE None,
FIRST_R1 FIRST_R1 VISIBLE None,
FIRST_ORIG_FID FIRST_ORIG_FID VISIBLE None,
Shape_Length Shape_Length VISIBLE None,
Shape_Area Shape_Area VISIBLE None)


texts: Running MaxEnt modelling in "dismo" using "loop" and storing the results in seperate objects and folders

texts:
I am doing a MaxEnt modelling using "dismo". I am using a loop to give multiple runs. I wanted to save results from individual runs in separate objects and at the same time the results from individual runs in separate directories.

The loop I have created runs without error. But when I am trying potential options to save the results, I am getting errors. The loop is as follows -

run = 2
for(i in 1:run){
group.k1 <- kfold(loxA.4, 5)
loxTrain.L1 = loxData.L1[group.k1 != 1,]
loxTest.L1 = loxData.L1[group.k1 == 1,]
background = randomPoints(predictor0, 10000)
pseudo_absence_values = extract(predictor0, background)
pseudo_absence_values <- as.data.frame(pseudo_absence_values)
pseudo_absence_values$africa_ecoregion <- as.factor(pseudo_absence_values$africa_ecoregion)
train_y = c(rep(1,nrow(loxTrain.L1)), rep(0,nrow(pseudo_absence_values)))
train_sdm_data = cbind(pa = train_y, rbind(loxTrain.L1, pseudo_absence_values))
Model.MT1 = maxent(train_sdm_data[,-1], p = train_y, path = "maxent/set_2/output"
show(Model.MT1)
Loop.r = predict(predictor0,Model.MT1, progress='text')
}


texts: How to get layer extent in OpenLayers or zoom to WMS layer extent?

texts:
I have added WMS layers using follow:

basemap = new OpenLayers.Layer.WMS( "Layer Name1",
"http://localhost:8080/geoserver/wms/workspacename", {layers: 'workspacename:layername'} );
utility = new OpenLayers.Layer.WMS( "Layer Name2",
"http://localhost:8080/geoserver/wms/workspacename", {layers: 'workspacename:layername', transparent:"true"}, {isBaseLayer:false} );

Now I want to get zoom to my utility layer, How can I?

I have tried as following:

map.zoomToExtent(utility.getDataExtent());


texts: Alternatives to the Vincenty geodesic functions in OpenLayers 3?

texts:
I'm migrating an application from OpenLayers 2 which uses OpenLayers.Util.destinationVincenty() to compute the projected latitude and longitude of a point from an origin given an angle and distance in metres.

destinationVincenty() and the other related geodesic utility functions are absent in OpenLayers 3.  I guess that they have been replaced by some other facility, but I'm struggling to find anything that helps.  I can't believe that I will have to re-code the Vincenty logic in my application, so...

Does anyone know of an embedded OL3 solution for geodesic arithmetic?


texts: Dropdown option in QGIS 3 processing modeler input?

texts:
I have an Excel spreadsheet with multiple sheets of different XY database. I made a processing model in QGIS 3.10, in which I have used 'make points from table' tool to create a point shapefile.

But in my Excel spreadsheet, there are multiple sheets so I can only select the first sheet as an 'input'. So, how can I make a dropdown option to select any sheet from input Excel?


texts:
I am using org.geotools version 2.7.0 and I am trying to get some shapefiles to GML 3.1.1. Sadly I got an exception:

java.lang.RuntimeException: Encode failed for _Feature. Cause: INVALID_CHARACTER_ERR: An invalid or illegal XML character is specified.
at org.geotools.xml.impl.ElementEncodeExecutor.visit(ElementEncodeExecutor.java:132)
at org.geotools.xml.impl.BindingWalker$BindingExecutionChain.execute(BindingWalker.java:215)
at org.geotools.xml.impl.BindingWalker.walk(BindingWalker.java:181)
at org.geotools.xml.impl.BindingVisitorDispatch.walk(BindingVisitorDispatch.java:44)
at org.geotools.xml.impl.ElementEncoder.encode(ElementEncoder.java:83)
at org.geotools.xml.Encoder.encode(Encoder.java:1100)
at org.geotools.xml.Encoder.encode(Encoder.java:836)
at org.geotools.xml.Encoder.encode(Encoder.java:564)

Maybe it's a space or an accented char (shape's attributes are full of accented chars) but I am very new to geotools and GIS, this error is completely blocking me. I can't find a way to specify a charset to read data from shape (by set up DataStore) or to write data into GML file (by set up Encoder). What am I missing?

I am trying to remove not allowed symbols from attributes and element names by using

StringBuilder sb = new StringBuilder();
if (!s.substring(0, 1).matches("[a-zA-Z_]"))
    throw new IllegalArgumentException();
sb.append(s.substring(0, 1));
for (int i = 1; i < s.length(); i++) {
    if (s.substring(i, i + 1).matches("[a-zA-Z0-9_.-]")) {
    sb.append(s.substring(i, i + 1));
    }
}
return sb.toString();

to sanitize shape feature attribute's values. Now attributes looks good but I am still getting the same exception.

I am new to GIS and org.geotools.


texts: MapInfo Thematic mapping - line width based on a field

texts: I want to apply thematic mapping to a layer in MapInfo with varying line thicknesses based on a float column. Is this possible?


texts: Issue to convert from PostgreSQL input to GPKG using Python GDAL API function gdal.VectorTranslate

texts:
When I run the following code, it fails with message

The error is ERROR 1: SELECT from table mytable failed, no such table/featureclass.

from osgeo import gdal
gdal.UseExceptions()

gdal.VectorTranslate(
    '/tmp/out.gpkg',
    "PG:host=localhost dbname=mydatabase user=username password=mypassword port=5432",
    SQLStatement='SELECT * FROM mytable',
    layerName='tablenameout'
)

The following alternate syntax failed too

gdal.VectorTranslate(
    '/tmp/out.gpkg',
    "PG:host=localhost dbname=mydatabase user=username password=mypassword port=5432",
    options='-sql "SELECT * FROM mytable" -nln layernameout'
)

With ogr2ogr and providing the same input table name and database connexion, the conversion happened correctly.
ogr2ogr -f "GPKG" /tmp/out.gpkg \
"PG:host=localhost user=username dbname=mydatabase password=mypassword port=5432" \
-sql "SELECT * FROM mytable" \
-nln tablenameout

What do I miss in my code syntax to make the Python API gdal.VectorTranslate works to open the database with my query without error?


texts: Integrate copyright label in qgis2web export?

texts:
I want to integrate an always visible link into a qgis2web leaflet map.
My ideas / I tried this:

view... decorations... copyright label = not show there
change the qgis project title, also with html = not the best solution, false position
try to change the index html = no success

Is it possible to take the copyright label within qgis/decoration?
It would be good, because I can use html for this field and it is always present after every map-export.

I use QGIS 2.18.12 (windows 10)


texts: Copying points fall on line to irregular polygon?

texts: I have points fall on  a centre line and I need to copy these points to the surrounding border, the polygon falls in irregular distance from the centre line?