texts: Two layers, same coordinate system, wont show up

texts:
I have vector two sets of vector layers that won't show up together. One set shows up together and the other set together etc but they are of the same area.

The layers I was given did not have a defined spatial reference which I manually projected as the same as the data frame I specified so all layers now have same CRS...yet not showing up together. I have noticed they appear to have different extents... What should I do?


texts: GeoJSON in Leaflet generated by gdal2tiles

texts:
I have generated a leaflet map with gdal2tiles out of my rasters and it works wonderful. Now I want to add my shapes in form of GeoJSONs and I have my difficulties.

To the gdal2tiles generated script, first I assign a variable to my geojson and show the location of it.

// GeoJSON
var Pad = L.geoJSON('./Padarchol/Padarchol.geoJSON');

Than I gather all my geoJSONs under one var, like it is done with the basemaps and the overlaymaps.

var shapes = {"Padarchol": Pad,}

In the end I add all of those layers

L.control.layers(basemaps, overlaymaps, shapes {collapsed: false}).addTo(map);

Result is a blank page. Without  these addings leaflet presents me a wonderful map...


texts: Leaflet and Tangram conctenated into a single script file

texts:
I am trying to concat leaflet and tangram source codes into a single script file using grunt to include it in an api which I am trying to make based on both leaflet and tangram. If I include the leaflet and tangram js files using different script tags in the html, the code works perfectly. But If i concat them using grunt and use a single script file I get the following error:

  ReferenceError: window is not defined

The example code which I tried to run is as follows:

var map = L.map('map').setView([12.96,77.58],14);
var layer = Tangram.leafletLayer({
scene: '../scene/scene.yaml',
      });
layer.addTo(map);
var marker = L.marker([12.96,77.58]).addTo(map);

This is a simple code to setup the map with a scene.yaml file and then add a marker to it. This code runs perfectly if separate script tags are used for both leaflet and tangram in that order. But the same code gives the above mentioned error when a single script tag is used with the concatenated script file with leaflet on the top and tangram following it.


texts: Reading raster values in points gives back a generator object instead of actual value

texts:
Given a raster TIF file and a point feature class in a geodatabase, I am trying to extract the values of the raster on the points and add them as a new column on the point feature class. I am using the following piece of code:

import ....
gdb_path = r'.......gdb'
pnts = gpd.read_file(gdb_path, driver='FileGDB', layer='...')
raster_path = r'.........tif'

#building list of pairs of coordinates
temp = [pnts['geometry'].x, pnts['geometry'].y]
coords = list(map(list, zip(*temp)))

with rasterio.open(raster_path) as src:
    pnts['values'] = src.sample(coords)

The code above creates a new column on the points feature class where the value for each point is a generator object sample_gen at 0x000..... instead of the actual raster value.

I could loop point by point instead, like this:

raster = rasterio.open(raster_path)
for id, row in pnts.iterrows():
    for val in raster.sample([(row['geometry'].x, row['geometry'].y)]):
        pnts.at[id, 'value'] = val[0]

This seems to work ok, but it takes a really long time (I have 4 millions of points). So, this solution is not valid for my purposes.

Could anyone shed some light on how can I make the first piece of code produce what I need? Or any other approach?


texts: I am setting up a hectare grid in QGIS 3.0 using Vector > research tools > create grid > grid type:point and then defining boundaries by a preexisting polygon layer. I would like to get coordinates from these points but the x and y axis created seem to be defined not by the CRS but by seemingly arbitrary numbers. Is there a way I can get coordinates into the attribute table? When I identify feature, I can find the actual lat/long coordinates but I cannot figure out how to import these into the attribute table.


texts: How do I delete an individual label without converting the entire data frame into annotation?


texts:
I have one shapefile that covers an entire city. I also have a list of 100 shapefiles which are buffers in different places in the city. Now I want to clip for each buffer from that city using ArcMap and python. I tried using arcpy in python but the code is not working.

import arcpy
from arcpy import env
from arcpy.sa import *
env.workspace = "U:\Park and Residential Area\Test\SBA park_res_buffer_5\SBA.gdb"
infeature= "U:\Park and Residential Area\Test\park_res_merge.shp"
clipfeatture = arcpy.ListFeatureClasses("*", "polygon")
for i in clipfeatture:
    outclipfeatture = arcpy.Clip_analysis(infeature,i)
    outclipfeatture.save("U:\Park and Residential Area\Test\SBA park_res_buffer_5/"*i)

This is the error i get:

Traceback (most recent call last):
File "<pyshell#10>", line 3, in <module> outclipfeatture.save("U:\Park and Residential Area\Test\SBA park_res_buffer_5/"*i)
AttributeError: 'Result' object has no attribute 'save'


texts: Seeking tutorial on Time Slider in Carto?

texts:
I want to add a time slider to a map. It looks as if I have to do that via CartoCSS but I really don't know where to start. Unfortunately, I don't understand what all the other posts regarding this issue say.

Does anyone know of a super basic tutorial on how to add a slider to a map?


texts: Clustering - Need to group points as I zoom out

texts:
I do not have much experience with ArcGIS Online but there is a work that must be done.

I have a groups of points and they represent addresses where a certain work has been or not performed. I need to develop a map to show, per area as the user zooms out, how many points have the work done and how many need to do it.

Showing it as a pie chart will be the ideal, but also with just numbers is OK (see example). So far I have tried analysis - summarize data - aggregate points, but the result is not what I need.

I cannot find a way to do the clustering without having to modify the code (because I am null programming).

Is there any way to do it?


texts: Creating route that passes from road only once using ArcGIS Network Analyst?

texts:
I want to calculate a route with ArcGIS Network Analyst that only passes once from every road instead of the proposed routes that network analyst produces now which return to the same road more than once.

Does anyone know the evaluators and needed expressions to calculate this?

I am thinking to add a new attribute in the network dataset properties, i.e called Uturns that will have cost as usage and will make it more costly for the calculated route to make a uturn instead of using the hole road till the next junction.

I came with an idea that if I use a script that won't allow the next road to be the same with the previous I might have the desired result. The script is below.

restricted = False

A= OBJECTID //current road object id

B= OBJECTID+1 //next road object id

C=OBJECTID-1//previous road object id

Select Case NoUturn([OBJECTID])

Case B=C : restricted = True

End Select

The question is how to assign values to the variables A.B.C, because I dont think that OBJECTID+1 is correct.


texts: Zonal average weighted by distance from point

texts:
I am using ArcGIS 10.2. I have a vector layer with points, which represent houses and a raster layer with values of 0 or 1, according to a certain threshold. In a buffer of 5 km around each of these houses I want to calculate the average number of raster points above a certain thresholds (value=1), weighted according to the squared distance from the house.

Without the weighting, this would be quite simple: I would itteratively make buffers around the house (to have buffers which do not dissolve) and then use the 'zonal statistics as table' command. Now, the weighting is complicating it.

I tried the following procedure (itterating house by house in a model, see below): first rasterize the houses (point to raster) then calculate the euclidian distance, then use raster calculator to take the square of the distance and to divide my raster with thresholds by the raster with squared distances. In the meantime I make buffers around the houses, rasterize these buffers and then finally run a zonal statistics as table command to get the statistics of the output of my second rasterization within the 5000 km. Finally I use the append command to make one single table for all the houses.

Now, I have two questions:

1) Is this a good approach or am I fundamentally wrong somewhere?

2) The analysis seems to work fine, but I get strange results for some houses (eg a maximum value of 5 for one house). Moreover, if I follow each of the previous steps without a model (manually) for some houses, I get different results. I don't understand why.

Could anyone help me?


texts:
Two friends and I are developing a series of web mapping applications, only with some basic tasks (measure tools, identify polygons, and others). We are looking for suggestions on what to implement to make our applications better. So, my question is:

¿What tasks do you add to your maps to add something extra?

As an example, we recently added a tool that get the distance between the selected points in a polygon, the selection is made by drawing a polygon in the map.

We are using ArcGIS Server 9.3.1 with .NET, but the tools listed in the answers could be implemented in any framework.

EDIT.

This is not only for us, everyone here could benefit of each idea/answer.


texts: What OSM information to use for low zoom levels to get a style similar to MapQuest

texts:
I am building an outdoor map of Europe. Everything looks pretty nice when zoomed in. But when I am zoomed out, my map looks pretty boring. Most parts are white.
In low zoom levels I display polygons only at a certain size. So there aren't many polygons that I could f.e. render as a green forest or grass-lands.

Does anybody know which information f.e. the MapQuest Map uses for zoom levels form 4 - 8?

Which polygons or which information is used to render most parts of Europe and America green, or the northern part of Africa brown, the poles white and so on?

I am using a Natural-Earth raster myself for the first few zoom levels, but after I can use those.

Which OSM information could I use?

Or are there some open shapefiles I could use?


texts:
I'm using 10.2.2 arcgis

I have the following script

import arcpy
mxd = arcpy.mapping.MapDocument ("CURRENT")
df = arcpy.mapping.ListDataFrames(mxd, "Layers")[0]
addLayer = arcpy.mapping.Layer(r"C:\Users\Document\Geomarketing\Datos\20140826\rest\PG_Restaurantes_GT.shp")
arcpy.mapping.AddLayer(df, addLayer, "BOTTOM")
del mxd, addLayer

The function of this script is to add a specific layer to my dataframe.  I want to add this script into a model, my purpose is to have a checkbox to have the option to add or not to add the layer.

I am new to this and do not know what changes should I make to my script to let me do this in modelbuider.

Do I do this in my script?


texts: Requirement to use a spatial database - Why use a spatial database?

texts:
The background to this question comes from a specific requirement I have in trying manage borehole (point) data.

Currently I have an MSAccess database which holds borehole information, the key data being borehole name, x, y, z, and total depth. My desire at the moment is to be able to have the boreholes in the database link directly to QGIS so I don't have to export weekly CSV's from Access to update the current borehole shapefile, or create a new shapefile.

I would just like the link to exist so that when I add a new borehole in the database, it will add into my mapping workspace also.

As I have searched for ways to 'live' link to my MSAccess datasheet however I keep finding references to using a spatial database for my data. What I don't understand is how a spatial database, like PostGIS which I see mentioned a lot, practically varies from my current database.

Does a spatial database such as PostGIS hold table data with z,y,z, and other attributes, or does it hold the actual spatial files, i.e. shp files.

If it does hold the actual files then how does that work in with the folder structure I currently use to store my GIS data?

I have spent the last couple of weeks searching and reading but I fear there is some logic gap I am just not making to put together how these pieces all relate and are used effectively.


texts:
I have elevation rasters which loaded just fine in a GRASS session. They represent adjacent areas, so I want to meld them into one raster.

I ran r.patch using the GRASS GUI, and it would only accept one input raster at a time in the input files field. So I used the command line in a terminal window. The result was an error message saying that the rasters: cannot be found -- the directory does not exist.

To run r.patch, do the input files have to reside in the Location folder I'm using for the GRASS session? That doesn't seem right. GRASS had no problem finding the same files to simply load them into the map window.

I'm running GRASS on Ubuntu machines (11.10 Oneiric Ocelot). Thank you.


texts:
I'm newbie to QGIS and to this forum.

I want to develop application using Qt on QGIS to real time updation/mapping of a symbol using GPS co-ordinates.

Initially i have installed QGIS on Windows 32 bit machine.

I'm looking into how to make vector maps in QGIS and the same map to be called on Qt Application to pin/map the location on map with a symbol reference by GPS Co-ordinates.


texts: Selecting a feature using expression if attribute field is "NULL"

texts:
I am new to QGIS and am having, what I expect to be, an easy issue. I am working with a multi-polygon shapefile of taxable property parcels in QGIS 3.16.
All it is I am trying to do is select all features whose "FieldA" is empty (or "Null")
In the expression field I have tried: "FieldA" = 0; or "FieldA" = NULL
Neither of these work obviously and I am at a bit of an impasse.


texts: Multipoint routing logistics software solution?

texts: I work for a state health program and am looking for the most efficient way to route inspectors from restaurant to restaurant. This is something like the travelling salesman problem, but with additional complexities. We have 75 inspectors, making less than 8 hour round trips from 17 different offices, inspecting ~18,000 businesses at 1.25 hours an inspection (not counting travel time). These businesses either need inspected annually, 3 times a year, or 5 times a year depending on their assigned risk level. Does anyone have any ideas/suggestions for this? I am talking to Telogis and Viamente as far as getting quotes for using their products but am interested in looking into all avenues. It's possible this is not on the right Stackexchange, if you know of a more appropriate venue I'd appreciate knowing that as well.


texts: Striping on linear fit scale image in Google Earth Engine

texts:
I have put together a linear fit image for two arbitrary geometries in North America and Europe. I am seeing striping in these images when I look at the Grayscale single band "scale". Does anyone know why this is occurring? It doesn't make sense to me over this time series.

Here is a sample of code for yearly data and then the ndvi collection that was made for the linear fit.

var terra_2011 = terra.filterDate('2011-01-01', '2011-12-31').map(addTime);
var terra_bands_2011 = terra_2011.select('tmmx','pr','system:time_start');
var precip_accumulation_2011 = terra_bands_2011.select('system:time_start','pr').sum();
//Create min, max and Average tempearture bands and rename bands
var temp_minmax_2011 = terra_bands_2011.select('tmmx','system:time_start').reduce(ee.Reducer.minMax());
//print(temp_minmax_2011, 'max 2011');
var temp_avg_2011 = terra_bands_2011.select('tmmx','system:time_start').reduce(ee.Reducer.mean());
//print(temp_avg_2011, 'temp avg');
//Combine precip and temp bands into one "image"
var merge_climate_2011 = precip_accumulation_2011.addBands(temp_minmax_2011).addBands(temp_avg_2011);
//print(merge_climate_2011);
//Create NDVI composite image
var collection_2011= ee.ImageCollection('LANDSAT/LT05/C01/T1_TOA')
  .filterDate('2011-01-01','2011-12-31')
  .map(addQualityBand);
var only_NDVI_2011 = collection_2011.select('nd','system:time_start');
var ninety_five_percent_2011 = only_NDVI_2011.reduce(ee.Reducer.percentile([95]));
var eighty_five_percent_2011 = only_NDVI_2011.reduce(ee.Reducer.percentile([85]));
var median_2011 = only_NDVI_2011.reduce(ee.Reducer.median());
var mask95_2011 = ninety_five_percent_2011.select('nd_p95').gte(0.1);
var lte95_2011 = ninety_five_percent_2011.updateMask(mask95_2011);
var mask85_2011 = eighty_five_percent_2011.select('nd_p85').gte(0.1);
var lte85_2011 = eighty_five_percent_2011.updateMask(mask85_2011);
var mask50_2011 = median_2011.select('nd_median').gte(0.1);
var lte50_2011 = median_2011.updateMask(mask50_2011);
//Map.addLayer(lte95_2011, imageVisParam3, "95 mask_2011");

//Create Collection and Apply Linear Fit
var ndvi_collection = ee.ImageCollection([lte95_1990, lte95_1991,lte95_1992,lte95_1993,
  lte95_1994,lte95_1995, lte95_1996,lte95_1997,lte95_1998, lte95_1999,lte95_2000, lte95_2001,lte95_2002,
  lte95_2003,lte95_2004, lte95_2005,lte95_2006, lte95_2007,lte95_2008, lte95_2009,lte95_2010, lte95_2011]);
var linear_fit_ndvi = ndvi_collection.select(['system:time_start_p95', 'nd_p95'])
  .reduce(ee.Reducer.linearFit());
var fit_year = linear_fit_ndvi.expression(
  '31536000000 * SCALE', {
  'SCALE': linear_fit_ndvi.select('scale')
  });
////print(linear_fit_ndvi, 'linear fit');
//Map.addLayer(linear_fit_ndvi.clip(geometry), imageVisParam4, 'ndvi fit');
Map.addLayer(fit_year.clip(geometry), imageVisParam5, 'scaled fit');

Here is the link to the full script: https://code.earthengine.google.com/9d0aba5699739485b44d9dc2c2a2f79c


texts: How to make lines invisible for displaying solely labels in GeoServer?

texts:
I am creating a road map for the city. I would like to display only road names on labels. For this purpose, i have created a style. In this style, i have changed all of stroke-width's to 0 (zero) but it didn't work.

          <CssParameter name="stroke-width">
            <ogc:Literal>0</ogc:Literal>
          </CssParameter>

How can i do it?


texts:
According to Rasterizing shapefiles with GDAL and Python? it seems that gdal allows rasterization of a layer according to a certain attribute

ds = gdal.RasterizeLayer(target_ds, [1], source_layer, options = ["ATTRIBUTE=ID"])

I would like to rasterize a layer with small square polygon geometries using rasterio

rasterio has a features.rasterize method (see docs), but it seems to only be able to output numpy arrays with 1's where the polygons occur. Since rasterio is a higher level library than gdal I am a bit surprised that selecting an attribute to burn the raster with is not supported functionality.

Is there any built-in method in rasterio to burn a raster by an attribute from the original vector layer?


texts: I'm script challenged and trying to repair a shapefile with a missing .shx. This tool sounds ideal: https://code.google.com/p/pyshp/, but the four step instructions are not enough. Can someone give me a quick walk through how you would run this tool.


texts: How do I import QGIS C++ libraries into my Qt/C++ project?

texts:
I'm working on Ubuntu 14.04 and I downloaded QGIS Desktop version 2.18.14...

I also have Qt Creator version 3.0.1 (based on Qt 5.2.1)...

When I downloaded QGIS Desktop via the Ubuntu software center, shared libraries such as libqgis-core2.0.1 were also downloaded....

I can't seem to find the paths to these libraries....

How can I import/include these libraries for a C++ project in Qt Creator or is there some other way to use the QGIS C++ API?


texts:
I'm trying to setup QGIS server as a local host for now so that I can test it out and see if it's what after. I'm completely new to server setup, but I've managed to get the server to respond to a Get Capabilities request such as:

http://localhost/cgi-bin/qgis_mapserv.fcgi.exe?SERVICE=WMS&VERSION=1.3.0&REQUEST=GetCapabilities

But when I try to launch a Get Capabilities request for a QGIS project stored on the server such as this one:

http://localhost/cgi-bin/?MAP=/cgi-bin/home/qgis/projects/world/world.qgs&SERVICE=WMS&VERSION=1.3.0&REQUEST=GetCapabilities/qgis_mapserv.fcgi.exe

Then I get an "Access Forbidden 403" response and the server produces the error log that I've attached below.

Where am I going wrong here?


texts:
For a web GIS application, we need to get all of our data and basemaps in the same projection.

We would like to use OpenTopoMap basemap which is in Web Mercator projection (EPSG 3857) but for our web GIS app, it has to be in Lambert93 (EPSG 2154).

How is it possible to reproject OpenTopoMap easily, if possible, with open source solution to create my own webservice of OpenTopoMap in Lambert93?


texts:
I am attempting to determine the best way to calculate the Minimum and Maximum elevation for a defined area.

The area is defined by an IPolygon

The elevation is defined as an ISurface

Possible Solutions:

The only way I can think of is to check every cell of the surface.  Is there a method that performs this?
In the comments, @whuber, suggested a Zonal Maximum.  Expanding on this idea, I've begun looking into the IZonalOp.

Are there other options?


texts: Have Legends been implemented for ESRI's REST API?

texts:
ESRI's 10.0 documentation for legend indicates that this url should return a legend:

http://map.floridadisaster.org/GIS/rest/services/Facilities/Critical_Facilities/MapServer/legend

But it doesn't work.  Is there something special I need to do when publishing a mapservice to enable legends?  If not, does anyone know of a url that returns a legend?


texts: ArcGIS not connecting to PostgreSQL db with all upper-case name

texts:
We've been trialling ArcGIS at work and are just about to buy. However, we have a PostgreSQL database which contains a lot of stuff we work with, which is named in all upper-case ("GISSTUFF", for instance). ArcGIS won't connect to it.

We've contacted the ESRI support desk and they say it needs renaming to lower-case. Which is easy enough to accomplish, but unfortunately throws out our other connections. GISSTUFF suddenly becomes gisstuff and their connections are lost.

And that's a lot of connections across a whole range of things. It can be managed, but we'd have to think of every point first where the connection would be, re-establish it and test it. I have no doubt that, since none of the staff here were the original system admins and documentation/training was sparse, things would fall through the gaps.

So, does anyone know of a way to make ArcGIS (just the standard Desktop stuff - I think that would be Pro nowadays) handle PostgreSQL databases in all upper-case?


texts: How can I switch to layout view in ArcScene and export a map in jpg format?


texts: Wildlife surveys using tablet and GPS, import/export to/from ArcMap

texts:
I've found some info about this on a previous post (Advanced android app for mapping and GPS logging?), but it didn't provide enough info, I'm hoping some other tech experts can help!

I have a tablet (Android) and it has integrated GPS, I would like to run a (preferably free) program into which I can input polygon layers (boundary units) and then track our flight path as we move through those polygon units.  We are talking fairly large scale (units are at least 30 km²).  There's no wifi, and we require something that keeps up to 100kph.  Our current GPSs do this fine, but we can't import polygons into the GPSs.  We have used ArcPad in the past, but it is too slow.

To summarize - we want a tablet-based program to which we can import GIS polygon layers and then watch our track on the tablet as we fly transects in the polygon units.  We would like to export these track to ArcMap when we're done, so pretty much any format will work.  I've looked into a couple options (Ozi, Google Earth), but I'd like to know if anyone has any experience with tracking GPS on a tablet at higher speeds and across large areas.

Thanks!


texts: Layers with the same CRS in QGIS do not line up

texts:
I am busy with a project where I have to layer a .dxf and a .kml file but the layers are not lining up.  The .dxf file's coordinates are completely wrong and the scale also does not appear to math up with the .kml file's. I have tried separately saving each layer to a new CRS, enabled "On the Fly" and setting the project CRS.  The layer's properties show that they are both in the same CRS.
The original CRS was undefined so it was automatically assigned WGS84. I later changed it to Hartebeeshoek 94.
It might be possible that the company who created the .dxf used a custom projection, is there any way to tell?  Can anyone please help?


texts:
I have a raster layer which contains mean erosion values in tonnes/hectare/year. The layer is in EPSG:4283-GDA94 with degrees as units (the green layer below). When I save it (i.e. reproject it) to EPSG:28356 - GDA/MGA zone 56 - projected, the values in the raster change (the blue raster below). Please see the image below.

My question is: What calculation has been done on the values in the raster during the reprojection from degrees to metres. I need to know how the values in the reprojected raster are related to the source raster (which is in degrees) so that I can use the data in some raster calculations.

Obviously during the reprojection the cell size will change but how does QGIS re-calculate each pixel values ?


texts:
i am trying to run GDAL in Ipython Notebook but i don't get it.

I did conda install GDAL and it was aparently ok, but when i try to do the import i get this:

ImportError                               Traceback (most recent call last)
<ipython-input-6-27bf4694dd2b> in <module>()
----> 1 import gdal

C:\Users\Usuario\Anaconda\lib\site-packages\gdal.py in <module>()
      1 # import osgeo.gdal as a convenience
----> 2 from osgeo.gdal import deprecation_warn
      3 deprecation_warn('gdal')
      4
      5 from osgeo.gdal import *

C:\Users\Usuario\Anaconda\lib\site-packages\osgeo\__init__.py in <module>()
     19                 fp.close()
     20             return _mod
---> 21     _gdal = swig_import_helper()
     22     del swig_import_helper
     23 else:

C:\Users\Usuario\Anaconda\lib\site-packages\osgeo\__init__.py in swig_import_helper()
     15         if fp is not None:
     16             try:
---> 17                 _mod = imp.load_module('_gdal', fp, pathname, description)
     18             finally:
     19                 fp.close()

ImportError: DLL load failed: No se encontró el proceso especificado.

Any idea please?


texts: Select SDE.ST_GEOMETRY polylines that have arcs?

texts:
Oracle 18c/10.7.1 EGDB:
I have an SDE.ST_GEOMETRY polyline feature class.
Some of the polylines in the feature class have arcs.

Question:
Is there a way to use Oracle SDE.ST_GEOEMTRY functions to select polylines with arcs?
I've looked through the docs, but haven't been able to find any functions that pertain to segment types.


texts: Keep raster data in original projection and just display it using a new transformation (at the poles) in R vs ArcMap

texts:
Desired output: Raster with original characteristics but in transformed projection using R (I.e. See second paragraph here.). NB, I need to do this for 100s of different rasters.

You can achieve this in ArcMap for a single raster doing the following. 1. Load original raster. 2. Right-click background of data frame and change coordinate system to new desired coordinate system (see figure below). Panel on the left shows original raster with cell dimensions (1 deg long by 2 deg lat). Panel on the right shows desired output, where the overall map projection has been updated and cell dimensions of raster have been preserved.

To achieve this in R, all roads point me toward using the function: projectRaster.

Where I have specified the following (r_toconvert = the raster I wish to convert; Shp_with_desired_projection = the shapefile with my desired projection):

r.new <-  projectRaster(r_toconvert, crs=CRS(projection(Shp_with_desired_projection)), method = "bilinear")

Note: my original raster is in projection:

+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0

And the desired projection for display purposes I'd like to be:

+proj=stere +lat_0=-90 +lat_ts=-71 +lon_0=0 +k=1 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0

However, when using projectRaster I produce rasters that do not preserve the same characteristics as my original raster (plotted in ArcMap for ease of exploration).

What I would like to achieve:

But what I actually end up achieving using projectRaster and trying both methods:

This seems to be somewhat non-trivial based on the ways rasters require re-projection (as outlined here section 6.6)

However, it clearly seems to be possible to preserve the original raster's characteristics (as evident above and in the ESRI webpage linked above) and plot with a new transformation.

WHAT HAVE I TRIED:

tmap package (which can also project on the fly):

tm_shape(antarctica, is.master = T)+ # Set is.master to specify this is the master projection
  tm_borders(col="black")+
  tm_polygons(col="grey")+
  tm_shape(r_toconvert)+ # My original raster
  tm_raster("Var", palette = terrain.colors(10))

Where the output is far from ideal:

The output is improved upon when using the new raster following projecRaster, and having the method set to ngb (method set to ngb does not seem to be the correct thing to do for a continuous variable though...):

r.new <-  projectRaster(r_toconvert, crs=CRS(projection(antarctica)), method = "ngb")

tm_shape(antarctica, is.master = T)+ # Set is.master to specify this is the master projection
  tm_borders(col="black")+
  tm_polygons(col="grey")+
  tm_shape(r.new)+ # My original raster
  tm_raster("Var", palette = terrain.colors(10))

But there still appears to be the same issue with new raster now losing characteristics of original raster (hence the jagged cell shapes above).

I've also played with some R base plot options, and looked in the SOmap package currently in development

I seek an R solution.


texts: Is there a R package that does spatial regularization?

texts:
By spatial regularization I mean the method described in

Aitkinson and Curran (1995) "Defining an Optimal Size of Support
for Remote Sensing Investigations", or
Section 4.8 of Webster and Oliver (2007) Geostatistics for Environmental Scientists

which is a way of calculating the variogram of a spatial variable that has a different support to the observed one.


texts: How to change the style of a layer in a mapbox?

texts: I would like to change the style of an arbitrary layer in a mapbox, but I can't. I can get the desired layer name, but how to set up this layer style?


texts:
I have been working with GeoWebCache on a out of the box build GeoServer 2.17.0.

I have set "Enable direct integration with GeoServer WMS" in the Caching Defaults.

When I call GeoWebCache directly (http://localhost:8080/geoserver/nurc/gwc/service/wms?service=WMS&version=1.1.0&request=GetMap&layers=nurc:Arc_Sample&bbox=-180.0,-90.0,180.0,90.0&width=768&height=384&srs=EPSG:4326&format=image/jpeg)

I get the following error: 400: Requested horizontal resolution: 1.40625 , best match: 0.703125 exceeds 10% threshold. Perhaps the client is configured with an incorrect set of scales (resolutions),or the DPI setting is off compared to the one in GWC ?

The 0.703125 number is the resolution for the EPSG:4326 srs which is pre-defined and set in stone AFAICT. I am not sure where the 1.40625 number comes from, I assume it is calculated based on the request.

When I make a normal WMS request (http://localhost:8080/geoserver/nurc/wms?service=WMS&version=1.1.0&TILED=true&request=GetMap&width=360&height=180&layers=nurc:Arc_Sample&bbox=-180.0,-90.0,180.0,90.0&srs=EPSG:4326&format=image/jpeg) I get the following in the Return headers:

geowebcache-cache-result: MISS
geowebcache-miss-reason: request does not align to grid(s) 'EPSG:4326'

I have tried a few different requests but the result is always the same.

I have seeded the first 4 zoom levels for the layer

Set 'Compute from SRS bounds' for the layer's bounding box settings and dropped and recreated the layer and the cache layer

Alas none of it has changed anything.

Do you have any insight or ideas?


texts: How to get size of file geodatabase feature class on disk?

texts:
Is there a simple way to determine the amount of hard drive space a feature class takes up?

I feel like I'm missing something simple, but I don't recall a method that does so.

You would think that right-clicking a dataset in ArcCatalog would do it for you.

The best I've ever been able to do was get the size of the entire gdb by looking in Windows Explorer.

If there isn't a way to do this in the ArcCatalog UI, I'd still be interested in a method to do so in code.

Any ideas?


texts: Spatially join geojson point and polygon


texts: Merge layers that have different coordinate systems

texts:
I need to merge two shapefiles and the issue is that one shapefile has WGS 84 coordinate system and the second OSGB 1936.

When I try to do the trick with MMQGIS -> Combine -> Merge layers I get an output shapefile but the two layers do not match at all.


texts: Reclassifying raster layer with certain colors for value ranges in QGIS?

texts:
I'm using QGIS. I want a raster layer with 6 different colors each color representing a certain value range (e.g. -1 to 0 red, 0 to 0.2 orange).

In the layer style menu I can only approximate colors.


texts:
As OpenLayers supports raster reprojection I would like to make it reproject xyz-tilelayers e.g. from OSM, Mercator projection to Mollweide projection (ESRI:53009).

From this example I took the necessary code snippets for using Mollweide projection, although in the example vector data is used.

I do not understand why the extent is set to [-9009954.605703328, -9009954.605703328, 9009954.605703328, 9009954.605703328], but it seems to work as the whole world is displayed.

In my case, reprojection itself seems to work, too, but all I get is a limited view of the world:

Furthermore, tiles disappear completely on some panning and zooming (out) interactions.

What am I doing wrong?

var projCode = 'ESRI:53009';
proj4.defs(projCode, '+proj=moll +lon_0=0 +x_0=0 +y_0=0 +a=6371000 +b=6371000 +units=m no_defs');
var mollweide = new ol.proj.Projection({
        code: projCode,
        extent: [-9009954.605703328, -9009954.605703328, 9009954.605703328, 9009954.605703328],
        worldExtent: [-179, -89.99, 179, 89.99]
});

var map = new ol.Map({
        target: 'map',
        layers: [
            new ol.layer.Tile({
                source: new ol.source.OSM()
            })
        ],
        view: new ol.View({
            center: [10,10],
            projection: mollweide,
            zoom: 0
        })
});


texts: Using similar color on points and lines in ArcMap or FME?

texts:
I have two feature classes/shapefiles. One containing points and one containing lines.

They have both a field called Gate_number. I can use Symbology to give every point with the same gate_number similar color, and i can do the same with the shp containing lines.

My problem is that I want to have the same color on the lines and ponits which have similar gate_number.

Every point with gate_number=2155 should be green and ever line with gate_number=2155 should be the same color green.

Can this be done in ArcMap or FME?


texts: Getting Latitude and Longitude from name of intersection point?

texts:
For example, my intersection is Cross Avenue @ Sacramento Street, in the City of Tulare, CA.

I need a tool that gives me a GPS location for like 1500 intersections.


texts: Error when calculating X/Y Coordinate for point feature class

texts:
I am receiving an error when I try to calculate X/Y coordinates, and I cannot understand why the process is not working. I am working with a point feature class in an SDE database. The point feature class has about 3,255 records. Within the point feature class I have two fields to represent the X coordinate and the Y coordinate. Both of these fields are a double field type with precision 9 and scale 5. Both fields are set to allow NULL values, and no domains are assigned.

When I go to calculate geometry I choose the coordinate, use the current data projection, and select US feet. Everything up to this point seems fine, however, as soon as I initiate the calculation, ArcMap kicks back the following error:

In this case, the feature with Object ID of 3256 is the first record in the feature class. If I select yes then I get this error message:

Any ideas as to why this is happening? Am I using the wrong field type or do I have the wrong precision/scale for the double fielt type?


texts: Is it also possible to graduate labels like it is possible for markers?

texts:
QGIS Dufour provides an enhanced vector-styling interface but is it also possible to graduate the labels like it is possible for the markers?

I have cities with population data: I can have the marker size according to the population value but is it also possible to adjust the lettersize of the label according the population?


texts: Striping in ArcMap picture symbol PDF export?

texts:
I made two PNGs to use as point symbols for state highways. When I export them, the symbols have stripes.

I've tried embedding all document fonts and converting marker symbols to polygons with similar results.

How I can export these symbols properly?


texts: Exporting layout in .pdf A3 size using PyQGIS?

texts:
I am using PyQGIS with GGIS 3.

I'd like to export my map in a A3 size .pdf (so that my printer prints it directely in A3).

How do I do that?

I did not find the right settings for exporting the pdf.

layout = QgsProject.instance().layoutManager().layoutByName('Example')
pc = layout.pageCollection()
pc.page(0).setPageSize('A3', QgsLayoutItemPage.Orientation.Landscape)

exporter = QgsLayoutExporter(layout)

pdf_settings = QgsLayoutExporter.PdfExportSettings()
pdf_settings.dpi = 300
#pdf_settings.papersize = 'A3'  --> is there something alike?

exporter.exportToPdf(result_folder +"\\myA3.pdf", pdf_settings)


texts: Is there a way to join specific features of 2 seperate line layers into 1 new layer?

texts: Specifically, I am looking to put speed limits (which currently reside in my road network layer), into my bus route layer (which follows some of the same roads, but not all within the road network). The merging vector layer option could be useful (once I manually filter the results in excel), but I am wondering if there is a better way to do this.


texts: OpenLayers v3.10.1 problem getting vector layer features unless layer visibility initially set to true

texts:
This is driving me crazy because I had it working with an earlier version: I'm trying to use the feature extents of a hidden (off) layer to set my zoom-to code for a given feature, but the only way I can get it to work now is to initially load the layer with visibility 'true':

areas_extents_src = new ol.source.Vector({
            crossOrigin : 'anonymous',
            projection : 'EPSG:3857',
            format : new ol.format.GeoJSON(),
            url : wfsURL + 'jbwq:sample-envelope-200m&outputFormat=json'
});

sampling_areas_extents = new ol.layer.Vector({
            title : 'Sampling Area Boundaries',
            name : 'sampling_areas_extents',
            source : areas_extents_src,
            visible : true,
            style : extentsOn
});// end SAE wfs layer

then do this:

sampling_areas_extents.once('postcompose', function(event){
    map.render();
    sampling_areas_extents.setVisible(false);
});

which just seems wrong! But no amount of hacking - turning the layer on before the rest of the zoom function, trying to clear/reload source, loading the layer with a hidden style, etc., will accomplish what I need. Is this a bug, or am I missing something?


texts: QGIS: CSV point grid to polygon - not raster

texts:
I have a CSV file with various data in 1 degree wide by 0.5 degree high rectangles, for which I have the N,S,E,W, & centroid lat/longs, i.e. every bit of spatial data required. I'm trying to convert this into a polygon grid, so that I can display the associated data in the variety of ways that raster doesn't do nearly as well (for ref: I have successfully converted to raster using this method (point > shape > raster)).

Other Q&As recommend points2one plugin, however this results in a big mess of triangles; in addition, no matter how I import the csv, qgis will only display one point per cell, either the centroid (manually selecting centroid lat/long) or the NE corner (default, qgis picks up the NORTH & EAST field names).

I've just successfully done this using WKT, creating a formula in excel to generate the polygons from the NESW coordinates, which seems like the best way so far.

My question is: is there a better way of doing this natively in QGis? It feels like the software has all the info available to be able to do this with one or two clicks - am I just missing it?

Thanks!

p.s. if anyone reads this & is in the same boat, my formula was:
="POLYGON (("&E2&" "&D2&", "&E2&" "&F2&", "&G2&" "&F2&", "&G2&" "&D2&", "&E2&" "&D2&"))"
with S/W/N/E in columns D/E/F/G respectively, and the formula in H. Remember to copy * paste back over the output as values before you save.


texts: Set same dimensions multipolygon and point data

texts:
I'm trying to plot point data (coördinates) on a multipolygoon map. But I can't get them on top of each other and in the right shape.

Underlying projection data of the Netherlands multipolygon data:

Simple feature collection with 1 feature and 3 fields
geometry type:  MULTIPOLYGON
dimension:      XY
bbox:           xmin: 13565.4 ymin: 306846.2 xmax: 278026.1 ymax: 619232.6
epsg (SRID):    28992
proj4string:    +proj=sterea +lat_0=52.15616055555555 +lon_0=5.38763888888889 +k=0.9999079 +x_0=155000 +y_0=463000 +ellps=bessel +towgs84=565.4171,50.3319,465.5524,-0.398957,0.343988,-1.87740,4.0725 +units=m +no_defs

Downloadlink to a .zip with 4 Geopackage files from the Netherlands

Underlying point projection data of coordinates:

Simple feature collection with 1466 features and 4 fields
    geometry type:  POINT
    dimension:      XY
    bbox:           xmin: 3.389737 ymin: 50.77106 xmax: 7.187327 ymax: 53.44402
    epsg (SRID):    28992
    proj4string:    +proj=sterea +lat_0=52.15616055555555 +lon_0=5.38763888888889 +k=0.9999079 +x_0=155000 +y_0=463000 +ellps=bessel +towgs84=565.4171,50.3319,465.5524,-0.398957,0.343988,-1.87740,4.0725 +units=m +no_defs

some rows of data:

points <- structure(list(brouwerijnaam = c("Almelosche", "HertogJan", "HertogJan",
"Bavaria", "Bavaria", "Brand"), begin_datum_jaar = c(1987, 1981,
1981, 1719, 1719, 1340), eind_datum_jaar = c(1990, 0, 0, 0, 0,
0), plaats = c("Almelo", "Arcen", "Arcen", "Lieshout", "Lieshout",
"Wijlre"), geometry = structure(list(structure(c(52.35192, 6.660562
), class = c("XY", "POINT", "sfg")), structure(c(51.473427, 6.181219
), class = c("XY", "POINT", "sfg")), structure(c(51.473427, 6.181219
), class = c("XY", "POINT", "sfg")), structure(c(51.518607, 5.597416
), class = c("XY", "POINT", "sfg")), structure(c(51.518607, 5.597416
), class = c("XY", "POINT", "sfg")), structure(c(50.832695, 5.895593
), class = c("XY", "POINT", "sfg"))), class = c("sfc_POINT",
"sfc"), precision = 0, bbox = structure(c(xmin = 50.832695, ymin = 5.597416,
xmax = 52.35192, ymax = 6.660562), class = "bbox"), crs = structure(list(
    epsg = NA_integer_, proj4string = NA_character_), class = "crs"), n_empty = 0L)), row.names = c(NA,
-6L), sf_column = "geometry", agr = structure(c(brouwerijnaam = NA_integer_,
begin_datum_jaar = NA_integer_, eind_datum_jaar = NA_integer_,
plaats = NA_integer_), .Label = c("constant", "aggregate", "identity"
), class = "factor"), class = c("sf", "tbl_df", "tbl", "data.frame"
))

Code I used to shape the x and y coördinates in a sf object with the same crs of multipolygon data.

library(sf)
library(tidyverse)

points <- df %>%
  filter(!is.na(lon)) %>%
  st_as_sf(coords = c("lat", "lon"))

points <- st_set_crs(brouwerijen_sf, 28992)

ggplot() +
  geom_sf(data = multipolygon) +
  geom_sf(data = points) +
  coord_sf(datum = NA)


texts: Pathfinder Office - how to tell if an external antenna was used with GPS

texts: Does pathfinder office store the antenna name in any of the log files?  I'm trying to go back and compare old files to determine if data collected with an external antenna did any better than those using only the internal antenna.


texts:
Possible Duplicate:
  How to install QGIS on Ubuntu 12.04?

Title says it all, and if so, how do I go about getting it? I currently have 1.7 Wroclaw installed via the Ubuntu Software Center, but I need the more advanced Field Calculator of 1.8.

I am using Ubuntu out of necessity because the Windows version has a bug preventing large raster images being exported. I'm unfamiliar with all the ridiculous hoop-jumping and command typing seemingly required to simply install software in Ubuntu, so answer as if to a complete idiot and you won't be far off the truth!


texts: Problems with shared projects that include files in OneDrive

texts:
I used to share QGIS projects with colleagues by storing our files in a server, which everyone would have access to. Recently we started using OneDrive and this has become a problem with the projects that we all share, because if we store the files in a OneDrive and then the rest synchronise their folders, we will still have different paths. e.g.
We store the folder in my OneDrive
C:\Users\XXX My username XXX\OneDrive - XX Company Name XX\Projects\Working\XXX.qgz

Even if my colleagues synchronise the file in their local drives, they will have something like
C:\Users\XXX Their username XXX\XX Company Name XX\XX My name XX - XX Project Folder XX\XXX.qgz

This means QGIS will not be able to find the files, and they will not be able to use the project.
My question is, does anyone know how to use shared projects with others with files stored in OneDrive (or similar storage platforms)? Is there anything I can do in QGIS to address this issue? Or is the only solution to go back and store files in the server?


texts: Is it possible to edit draw panel in Leaflet Draw?

texts:
Can i edit leaflet draw panel by adding some new buttons?

For example, i want two buttons for drawing polylines, three for markers and etc.


texts: Leaflet: Is there a way to use click event just like mouseout?

texts:
I made a map with leaflet by following the official Leaflet tutorial on interactive choropleth maps.

I would like to make the map compatible with smartphones, so I added a click event, which acts like hovering on a feature.

function onEachFeature(feature, layer) {
layer.on({
    mouseover: highlightFeature,
    mouseout: resetHighlight,
    click: highlightFeature
});
}

The consequence: When a user is clicking, the hightlightFeatures remain highlighted.

Is there a way to prevent this? In other words: Is there an event on the lines of mouseout, in order to activate resetHighlight by clicking another feature?


texts: Synchronisation problems using QFieldSync

texts:
When resynchronizing data from QField GPKG project to QGis SHP project, using QFieldSync, attribute data are shifted to the right column. There is a lag, I have never seen something like that before.

Still, visualizing GPKG data using QGis Data Manager, without resynchronize, it can be seen that GPKG data are in compliance, it is found in the correct column.

QField GPKG Project screenshot

QGis SHP Project screenshot (therefore after resynchronising)

So, the issue comes from the resynchronising process using QFieldSync.

I don't really care but on-site coworkers that are not primarily directed at GIS will not understand it. Has anyone else come across a similar problem ?


texts: GeoNode WFS - restrict some layers to certain users

texts:
I´m using GeoNode 2.4, which uses GeoServer 2.7 and allows GeoNode staff and admin users to access GeoServer admin panel. My aim is to make some layers in GeoNode only downloadable by a few GeoNode users, and only accessible via WFS to those few users, while keeping all other layers free to everyone else. I can manage the GeoNode download permission, but the layers are still accessible in WFS service (any user can load them and "save as"). Here´s what I´v done:

I separated those layers in a new created workspace and created a new role. Then made a rule that allows reading the new workspace only to that role and admin. Now I should assign that role to those specific GeoNode users,  but I don´t know how to. In GeoServer admin panel I can´t even see the list of GeoNode users (even when I can log in to the GeoServer panel using GeoNode user credentials). When accessing WFS with a regular geonode user, I can´t see the restricted layers (which is fine) and when accessing with admin geonode user I can see all layers (which is fine too). I don´t want to give admin privileges to those users so they can access the restricted layers via WFS.

Any help?


texts: Changing chart properties, ArcMap 10.3 Desktop

texts:
I am attempting to change the properties of the attached chart to create something a little more aesthetic. After searching ESRI Help, message boards, google... I have not found any material on how to make these charts look better (I can't even find a way to change the font!).

Usually, I would just make the chart in MS Excel, however this chart will be part of a time series and needs to be linked to the data to update as time moves, so it cannot be static.

Do you know of any way to change these properties without unlinking it from my dataset?


texts: How can I simulate an event in OpenLayers 3?

texts:
I have the following OpenLayers 3 code for this:

map.on('dblclick', function(evt) {
    // Do Something.
});

The idea is catch the event principally. I've tried to find something but without success.

Any idea?


texts:
I have a shapefile with about ~5k polygons. I have a database of about 3.8M points.

I want to know, for every point, which polygon it's in.

For example, Point A may be in Polygon A, Polygon B, and Polygon C.

The problem I have is two fold:
1) How do I tie a specific point match to a specific polygon? Right now, I can only get the sum total of the number of points within  a polygon using the following code:

colSums(gContains(all_polygons, spdf, byid = TRUE)) -> res

2) How do I make this computationally easier? I am less familiar with GIS analysis in R, but in traditional SQL language, we may do a cross join of all polygons with all points, and then throw a flag 1,0 if the point is in the polygon. This is incredibly computationally intense, as it generates a table that is 5k* 3.8M. Is there another way around this?


texts: ArcGIS 10.2 objectid returning non sequential values

texts:
I am making a SQL insert in a geotable. I'm getting the objectid via:

sde.gdb_util.next_rowid('my_user', 'my_table')

However, I've noticed that numbers given are not always in sequence. For instance, if I run the script 6 times, I get:

run|oid
---|------
  1| 1601
  2| 1602
  3| 1603
  4| 1917
  5| 1918
  6| 1604

Is this normal?

I know that for ids the only important thing is uniqueness, and sequentiality is optional, but I'm intrigued.


texts:
I have an ESRI Geodatabase (Version unknown) which is not displayed in ArcCatalog as a gdb but as a normal directory with a ".gdb" ending. It could be, that whilst it was being unzipped, something got corrupted.

I have tried restoring it with the ArcGIS Pro "repair GDB" Tool but to no avail.

Is there a way of determining what is wrong or missing so that I could repair it?


texts:
I am creating an app with leaflet that shows a route on a map based on the button clicked.  For example when you click on "button1", the associated route is successfully displayed on the map (i.e. map.addLayer(route1Data); works). When you click on either button2 or button3, the route for button1 should be removed. For some reason map.removeLayer(route1Data) is not working.

I am testing button1 at this time.

Why is the layer not being removed?

Steps, 1. whichRoute() --> 2. Tracks() --> vlocation()

See scripts below:

HTML

<title>Route Chooser</title>
    <link rel="stylesheet" href="http://ajax.googleapis.com/ajax/libs/jqueryui/1.8/themes/black-tie/jquery-ui.css" type="text/css" />

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

    <button class='rtbttn' id='rt1' name='rt1' data-name="route_1">Route 1</button>
    <button class='rtbttn' id='rt2' name='rt2' data-name="route_2">Route 2</button>
    <button class='rtbttn' id='rt3' name='rt3' data-name="route_3">Route 3</button>

Javascript

function tracks(xx){
$.ajax({
    url:"../php/vTrackerv2.php",
    type: "POST",
    dataType: 'json',
    //cache: false,
    success:function(rdata)
    {
        var locData = [];
        for (x=0; x< rdata.length-27; x++) {
            var loc = [];
            loc.push(parseFloat(rdata[x].latitude));
            loc.push(parseFloat(rdata[x].longitude));
            locData.push(loc);
        }
    vlocation(locData,xx);
    }, //end success
  });
};

function whichRoute(){
$('#bttnSelect').on('click', function(event) {
    // Retrieve the 'name' data attribute of the  tag that the user clicked.
    var bttnClicked = $(event.target).closest('button').data('name');
    if (bttnClicked == 'route_1'){
        console.log("You clicked on:", bttnClicked);
        tracks(1);
    }

    if (bttnClicked == 'route_2'){
        tracks(2);
        console.log("You clicked on:", bttnClicked);
    }

    if (bttnClicked == 'route_3'){
        tracks(3);
        console.log("You clicked on:", bttnClicked);
    }

  });
}


texts: Hot Spot Analysis for an NDVI Time Series

texts:
I am working with a collection of Landsat 8 images with pixel-wise NDVI computations over a period of time and I'm hoping to pin point areas within the region of interest where NDVI change is greatest. I have looked into different Hot Spot Analysis techniques and ArcGIS has a Hot Spot Analysis tool that seems like it may be a viable route to go. However, I am not sure how to move forward with my Google Earth Engine derived data set to compute hot spots of change between two image periods and further down the road, hot spots of change over a full grazing period for the region of interest.

Does anyone have experience with this type of change detection analysis and would you be able to offer guidance? Thanks!


texts:
I imported shp to postgis 2.0.1 with shp2pgsql and chose to use geography instead of geometry. However, now when I try to drop the table, the pgadmin/sql console/phpPgAdmin freezes. I can't neither truncate it or change the columns.

Any help would be appreciated :-).


texts: Place label along line near start and end of line in QGIS

texts:
I have four fields managing the numbering of the start and end of road sections.

I created two label rules, one for the start and one for the end of the road sections.

I set the position of the labels to "parallel / on the line".

I tried to configure the position of the labels according to these expressions
x(end_point ($ geometry)) and y(end_point ($ geometry)) for the end of the road sections, and a similar expression for the start of the road sections.

The problem is that the labels no longer take into account the position "on the line" and certain labels no longer appear due to their overlapping.

This tool exists in arcmap which allows you to manage the labeling from the start or the end of the road section. Is there a similar tool in QGIS?

If not, do you have an idea for a script allowing you to configure these labels?


texts: Select by location is not working in the processing modeler

texts: When I run the modeler in qgis 2.14 I get no results from the select by location algorithm.  I tried outputting my buffer polygon and using the select by location plugin and it worked fine.  I have the same setting in both but for some reason the modeler doesn't do anything.


texts: How to get the TOC to update after adding a new Layer in the ESRI JSAPI?

texts:
I am developing an Web map with the ESRI's JavaScript API v 3.9.

I am using the TOC widget available at http://www.arcgis.com/home/item.html?id=9b6280a6bfb0430f8d1ebc969276b109

The documentation shows how you can initialize the TOC when creating the map. In my case, the user can dynamically add any ArcGIS Dynamic Mapservice by providing the URL.

I have tried the following:

Add the Layer to the toc's layerInfos and then calling toc.startup(); This does not do anything. Neither is there an error, nor do I see the new layer.
Delete the HTML inside the TOC Div, and try to recreate the widget with the new layerInfos. This gives me an error saying that the div has already been registered.

How do I show this new Layer in the TOC when it is already initialized?


texts:
I am trying to run this code:

# create virtual multi-band raster
gdalbuildvrt -separate rgb.vrt red.tif green.tif blue.tif
# write multi-band raster into geotiff
gdal_translate -of GTiff -co PHOTOMETRIC=RGB rgb.vrt rgb.tif

I want to know how to be able to do this in a terminal outside of QGIS.


texts:
My script should multiplicates the pixel value of a raster, but produce a raster with only 0.

import os
import os.path
import processing
from qgis.core import *
from qgis.analysis import QgsRasterCalculator, QgsRasterCalculatorEntry

rasterVH_filepath = "/home/giacomo/Desktop/divisione/VH1/VH1.tif"
output = "/home/giacomo/Desktop/divisione/out/prova.tif"

entries = []

myLayer1 = QgsRasterLayer(rasterVH_filepath, 'raster')
my1 = QgsRasterCalculatorEntry()
my1.ref = 'my@1'
my1.raster = myLayer1
my1.bandNumber = 1
entries.append( my1 )

calc = QgsRasterCalculator('my1@1 * 10',
output,
'GTiff',
myLayer1.extent(),
myLayer1.width(),
myLayer1.height(),
entries)

calc.processCalculation()


texts: I'm trying to populate a table with parcel information and assign sequential numbers while dealing with some parcels that have anywhere between 2-50+ entries.  One owner has many parcels.  Is there a way to query out just duplicates to identify them.


texts: Geographic constraints on tables in MapInfo

texts: In MapInfo, how can I remove polygons from a table according to whether or not points from another table fall within their boundaries?


texts: How can I set up a parameter as an attribute drop down menu on a script tool?

texts:
I have built a script tool in ArcGIS Pro and I would like to set up one parameter as a drop down menu (the options in the drop down menu should be all the attribute values of one of the map layers).

I can't seem to find the right 'Data type' to pick up. If I choose 'Layer' as 'Data Type' and I run the tool I get a drop down menu with all map layers. I'm looking for something similar but specific to a layer field. If I choose 'field' in 'Data Type' and I run the tool absolutely nothing happens. What is this 'field' referring to? How could I achieve what I'm looking for?


texts: Creating Raster to Point data using ArcMap ModelBuilder Iterate raster tool

texts:
Creating Raster to Point data using ArcGIS ModelBuilder Iterate raster

Raster material is in folder (ffaa)
Using Iterate rasters in model builder
Select folder 1 (ffaa) in step 2
Raster to point function input and input raster, output point feature setting
Set field to blank

Status
ArcToolbox's raster to point function creates 1 raster to point data.
However, it fails in ModelBuilder.
As a result of confirmation, there is no field selection in the raster to point of ModelBuilder.
I tried entering a value in the field name, but it failed.
Error message

ERROR 010151: No features found in raster. Possible empty feature
class. Failed to execute (Raster to Point).
ERROR 999999: Error executing function. Failed to execute (Raster to
Point).

desired goal
When I want to use the raster to point function in the ModelBuilder, I want to select the field I want in the field selection.
So I directly used Value, "", ==, %%, %Value%, !Value!, $Value$, all, && etc in the field field, but it didn't work.


texts: Can you make a travel region polygon with with Google Maps API?

texts:
I am working on a local design project and need to determine if someone is local to my area to allow them to participate. Local, as defined by the group, means that the person is within 2 hr driving time based on Google maps to 32 Vassar St., Cambridge, MA.

While I can figure out how to get travel times from point to point, I cannot seem to figure out if there is a way to create a polygon that displays the region which we consider local.

Any ideas on how to do this using Google tools?


texts:
SpatialJoin_analysis is not working with delimiter semicolon (;).  It only worked with a comma delimiter (,)

arcpy.SpatialJoin_analysis(Unique_Intersections,
                       Streets_Unsplit,
                       Signed_Intersections,
                       "JOIN_ONE_TO_ONE",
                       "KEEP_ALL",
                       "IntersectionName \"IntersectionName\" true true false 150 Text 0 0 ,Join,\" , \",G:\\xStreetNew\\Test.gdb\\RoadCenterline_Unsplit,STREETNM,-1,-1",
                       "INTERSECT",
                       "1 Meters", "")

If I change the delimiter (,) (highlighted in yellow above) to (;), it won't work.  Any ideas?


texts: Geocoding - get lat/long from 11000 address

texts:
What I have:

an Excel file with 11 thousand address (street + zip code + city)
all address are from Portugal
the address are not accurate (some address have the name street in a short way)

What I am looking for:

get "geocodes" (latitude/longitude) of all address.
using a opensource API, or a service of a website

What I will do with the data collected:

with the "geocodes" I will create a file like a KML file with all address and some others information about each address
the "geocodes" will be use in a commercial product by the company I am working.

In my research I have found some API and websites. But I am not sure if they suit my requirements

Google Maps API

I think it was a limit if 2500  requests per IP per day.
I think I can not collect the latitude/longitude, only show the location in a Google Map webpage

Yahoo PlaceFinder API

I think it is not free

website Gisgraphy

The results on my test are not good for Portugal address

website GPS Visualizer's Address Locator

uses Google and API Yahoo API, So have the same restrictions above

Some links I found:

http://blog.programmableweb.com/2012/06/21/7-free-geocoding-apis-google-bing-yahoo-and-mapquest/#comment-189818
https://developers.google.com/maps/articles/geocodestrat?hl=pt-PT#quota-limits


texts: What is a purpose of parameter 'separater' in mapnik csv plugin datasource?

texts:
What is the purpose of the parameter 'separator' in the Mapnik CSV plugin DataSource

Example of the code:

         parameters p;
        p["type"]="csv";
        p["inline"]=quotedStr;
        p["separator"]="|";

What is the purpose of last line?


texts: Turning first vertex of a line into a point

texts:
import arcpy
import math
import pythonaddins

class Line(object):
"""Implementation for Projects_addin.Project1_2 (Tool)"""
def __init__(self):
    self.enabled = True
    self.cursor = 3
    self.shape = "Line" # Can set to "Line", "Circle" or "Rectangle" for interactive shape drawing and to activate the onLine/Polygon/Circle event sinks.
def onLine(self, line_geometry):
    cur = arcpy.da.InsertCursor("NewL", ["SHAPE@"])
    cur.insertRow([line_geometry])
    arcpy.RefreshActiveView()

    fields = arcpy.ListFields('NewL')
    length = False
    heightM = False
    heightF = False
    for field in fields:
        if field.name == 'Length':
            length = True
        elif field.name == 'HeightM':
            heightM = True
        elif field.name == 'HeightF':
            heightF = True
    if length == False:
        arcpy.AddField_management('NewL', 'Length', 'DOUBLE')
        arcpy.RefreshActiveView()
    if heightM == False:
        arcpy.AddField_management('NewL', 'HeightM', 'DOUBLE')
        arcpy.RefreshActiveView()
    if heightF == False:
        arcpy.AddField_management('NewL', 'HeightF', 'DOUBLE')
        arcpy.RefreshActiveView()

    arcpy.SelectLayerByAttribute_management('NewL', 'CLEAR_SELECTION')
    arcpy.RefreshActiveView()

    fid = int(arcpy.GetCount_management('NewL').getOutput(0))-1
    arcpy.SelectLayerByAttribute_management('NewL', 'NEW_SELECTION', ''' FID = {} '''.format(fid))
    arcpy.RefreshActiveView()

    arcpy.CalculateField_management('NewL', 'Length', '!shape.length@meters!', 'PYTHON_9.3')
    arcpy.RefreshActiveView()
    arcpy.SelectLayerByLocation_management('Source_Package', 'INTERSECT', 'NewL')

    with arcpy.da.SearchCursor('Source_Package', ['SunElevati']) as cursor:
        for feat in cursor:
            new_heightM = math.tan(math.radians(feat[0]))
    arcpy.SelectLayerByAttribute_management('Source_Package', 'CLEAR_SELECTION')
    arcpy.RefreshActiveView()
    arcpy.CalculateField_management('NewL', 'HeightM', '!Length!*{}'.format(new_heightM), 'PYTHON_9.3')
    arcpy.RefreshActiveView()
    arcpy.CalculateField_management('NewL', 'HeightF', '!HeightM!*3.28089895', 'PYTHON_9.3')
    arcpy.RefreshActiveView()

This code generates a line with attribution depending on the polygon it is overlapping. Is it possible to take the first vertex of the line and create a point in another feature class, while simultaneously transferring the attribution?

I was thinking (line_geometry.firstPoint.X, line_geometry.firstPoint.Y) would create the coordinates, but I am still considering other possibilities.

Edit 1:

point = arcpy.Point(line_geometry.firstPoint.X, line_geometry.firstPoint.Y)
pt_cur = arcpy.da.InsertCursor(Point Feature Here, ['SHAPE@'])
pt_cur.insertRow([point])
arcpy.RefreshActiveView()

I added this to the end and I have created the point I desired, now I am working on transferring the attributes.

Edit 2:

I have also successfully added the attribution aspect.


texts: Leaflet with custom map (CRS) - how to convert JPG to TileLayer tiles

texts:
I have a JPG file with a map which isn't in a real scale, and can't be represented by real geographical latitude and longitude, so I so the Leaflet CRS solution and this seems like the way to go.

The problem is that the map file (JPG) size is ~30MB and can't be loaded in a single request, so I read a bit more and found the TileLayer solution.
I even found a demo .

Now, my only problem is how can I take the JPG file and convert it to layers? Is there an online / offline tool for that?
As far as I understand, all I need is to split the large image to 256px x 256px tiles, place it in directories in my server and TileLayer will do the rest?


texts:
In a leaflet application, I have a base map, and several overlay layers, each of which are showing GeoJSON. The GeoJSON is going update (client receives over a websocket) at 4Hz.

The GeoJSON consists of a good number of polylines. The area covered is pretty small, so I don't think that tiling the GeoJSON will help.

I don't think that their will be so much data that it will overwhelm the browser/client. But I want to know if anyone has tips on what could decrease the load on the browser.

I.e. I could should to delete the layer content each update and load fresh content.

Or I could as the server to give me vertex list for the new layout of the polylines and then do an update in the browser.

There are probably other ways to minimize work. (e.g. are there plugins for this problem?) nothing in the leaflet plugin list leaped out at me.

What tips do you have?


texts:
Oracle 18c/10.7.1 geodatabase:
My geodatabase has an owner/schema/user called OPERATIONS.
Operations staff have a direct db connection in ArcMap/Catalog to the OPERATIONS owner. This connection is needed so that staff can create and alter tables in OPERATIONS.
I want to prevent Operations staff from deleting tables in the OPERATIONS owner.
What are my options for doing this?


texts: Earth engine - ee.Terrain.slope() output at wrong resolution after DEM mosaicking

texts:
The big picture: I'm using Google Earth Engine (Javascript) to do a terrain correction on Sentinel-1 imagery, as in Andreas Vollrath's work. However, I want to use a higher-resolution Canadian digital elevation model, CDEM (see documentation here) instead of the SRTM DEM.
CDEM comes as an Image Collection rather than an image, so needs to be mosaicked into an image before being put into a larger terrain correction function, where I use function terrain.slope() to get the slope from the DEM. However, when I mosaic this dataset, terrain.slope() does not produce the slope output at the same high resolution of the CDEM dataset (though the values and data types look correct). However, if I do not use mosaic but call terrain.slope() on the mapped CDEM ImageCollection, it turns out correctly (but with no-data borders around individual images because of edge effects).
My questions: can I make a mosaic of the CDEM that I can use in the terrain.slope() function, and have an appropriately scaled slope layer as output that can be used for terrain correction? What part of the mosaicking process results in the lower resolution slope layer? Note that I've tried mosaicking by using mosaic() and a median reducer, and both lead to this problem.
You can view code that illustrates my problem here: https://code.earthengine.google.com/331ddb0c3324256ebd4ec1ea0db4402c
The section that is problematic is:
var dataset2 = ee.ImageCollection('NRCan/CDEM');
var elevation2 = dataset2.mosaic()   //somehow this line causes problems for the slope layer***
Map.addLayer(elevation2, elevationVis, "CDEM")  //as a DEM it shows up fine on the map

//calculate slope in radians
var alpha_sRad3 = ee.Terrain.slope(elevation2).select('slope')
            .multiply(Math.PI/180).setDefaultProjection(proj);
Map.addLayer(alpha_sRad3, {}, 'slope CDEM') //lower resolution than the elevation2 layer - why??***

Whereas, if you do the slope function on each individual image in the CDEM dataset, it works:
var alpha_sRad4 = dataset2.map(function(image){
  return ee.Terrain.slope(image).select('slope')
            .multiply(Math.PI/180).setDefaultProjection(proj);
  })
Map.addLayer(alpha_sRad4, {}, 'slope CDEM collection')

Note that 1) this problem also exists for terrain.aspect(), 2) I moved this question here from StackOverflow, and 3) the terrain correction is also discussed here: Terrain correction (flattening) of Sentinel 1 images
Thanks!!


texts: ERROR ArcPY parsing error SyntaxError: invalid syntax (line 2)

texts:
I'm trying to convert multiple lines to points so I typed the code bellow but it keeps on showing me this error:

  Parsing error SyntaxError: invalid syntax (line 2)

liste_ligne =["E1_H20A_H35_D21E","E1_I_H16A_H18B","E1_K_D16A_D18B","E1_L_H14_H16B","E1_H2_D70_D40C","G1_D50_D55_H40C","E1_H1_D60_D65_H70","E1_M_D12_H12","E1_Bbis_H45","E1_C_H18A","E1_C_H40","E1_D20A_D21A_D35","E1_E_D18A","E1_E_H50","E1_Bbis_D20E1","E1_G2_D21C_H60_H65","E1_LM_LoisirMoyen","E1_Bbis_D20E","E1_F2_D40_D45","E1_LL_LoisirLong","E1_ABIS_H20E_H21A","E1_F1_H55_H21C","E1_Lbis_D16B_D14"]
# loop the elements on the list
for ligne in liste_ligne:
# link to shapefile to join
inFeatures = "E:\shapes\"+ligne+".shp"
# link to output
outFeatureClass = " E:\E1TracesPoints\+ligne
# Execute FeatureVerticesToPoints
arcpy.FeatureVerticesToPoints_management(inFeatures, outFeatureClass, "MID")

What have I done wrong?


texts: How can I construct a query (for example for Overpass) for getting all buildings where inside the building is a node (for example a shop)?


texts: Function st_clip(raster, integer, geometry) does not exist

texts:
I am trying to replicate this example by alphabetasoup for summarizing pixel values according to census polygons, but am getting this error.

ERROR:  function st_clip(raster, integer, geometry) does not exist
LINE 8:   FROM (SELECT sa_id, ST_SummaryStats(ST_Clip(rast,1,geom)) ...
HINT:  No function matches the given name and argument types. You might need to add explicit type casts.
SQL state: 42883
Character: 176

The code is only changed slightly to reflect my particular schemas, the results I'm querying, and to make it more legible.
WITH
   polys AS (
       SELECT id As sa_id, geom
       FROM prac4.small_areas
    ),
   poly_stats AS (
       SELECT  sa_id, (stats).*
        FROM (SELECT sa_id, ST_SummaryStats(ST_Clip(rast,1,geom)) As stats
        FROM prac4.no2_raster
        INNER JOIN polys
        ON ST_Intersects(geom,rast)
    ) As foo
    )

SELECT sa_id, SUM(mean*count)/SUM(count) As avg_pval
FROM poly_stats
GROUP BY sa_id
ORDER BY sa_id;

A similar question has been posted here but I think the answers there had to do with incorrect SQL. I'm probably missing something obvious.
I'm on postgis 3.0.2., postgresql 12.4. The postgis_raster extension has been added to the database.


texts: How can I view all features in the attribute table?

texts:
I downloaded a large map area from OpenStreetMap. There are features I can select, but the I cannot view in the attribute table or filter. Per the link below, I tried to change the attribute table settings to "Show All Objects" and other options, but whatever I try, there are some features that are invisible to the table.
QGIS Attribute Table Not Displaying Data


texts: Spatial analysis on raster in order to get watersheds of approximate 2 ha size

texts:
I'm trying to analyse a raster in order to get watersheds at approximate 2 ha size, using flow direction/accumulation.

I've been using the steps outlined below, with initial conditions for cells starting at 2 ha. However whenever I do a raster calculator, I keep getting a raster with only one big black 0.

Fill

Flow Direction

Flow Accumulation

Raster Calculator @ 12500 cells

Focal Statistics

Watershed

Raster to polygon

Am I doing something fundamentally wrong in this build up? Basically I need watersheds as polygons with a size of 2 ha, using a DEM Raster.

Additional information about raster DEM:
Raster is about 40 mb, in a projected CRS and a cell size of 0,4 x 0,4 m covering about 420 ha.

X: 5992 Y: 4280

Float32 - Thirty two bit floating point


texts: Categorized stroke color of a circle in QGIS

texts: Is there a way in QGIS to make circles that have strokes with different colours and the colour is based on the categories in one of my data fields?


texts: Clip geometries in table based upon value

texts:
I have defined a view in postgres that returns a set of polygons and can be queried as follows:-

SELECT wkb_geometry, grade_count FROM ng_schema.grade_count_report [sql.where] order by grade_count asc

This returns me a set of polygons with varying values of grade_count

What I want to do next is find geometries that overlap and cut out the overlap from the geometry that has the lower grade_count. Where they were equal I would cut the geometry with the largest area.

How best to do this?

UDATE: To be clear when I say "has the lower grade_count" this is the lower grade count of the two intersecting geometries NOT the lower grade_count of all of them as the answer below has calculated.

So I need to find ALL intersections across ALL geometries and clip them all with the LOWEST grade_count of the two losing part of its geometry.


texts:
I'm trying to upload a 37 GB raster dataset to PostGIS. The raster is out-of-db (-R option), so the actual pixel values are not stored in PostGIS. I'm using this command:

raster2pgsql -s 25833 -C -I -R -t auto /data/norway/terrain/10m/*.dem dem.10m | psql -h localhost -U postgres -p 5432 -d norway

Everything works fine, but it's veeery slow: There are 254 raster files in my folder, and in 24 hours only 8 are added. With this speed it will take me a month to upload the dataset!

Is there any way I could speed up the process?


texts: Extract mean value of a raster into a grid shapefile

texts:
I have a soil map of an island and a Grid polygon. I want to extract the mean value of the soil map to polygons squares of grid.

I read this and this but Ι can't find the results.  I use Geospatial Modeling Environment

isectpolyrst(in="\~\Polygon_grid.shp", raster="C:\~\TAXGWRB_1km_T453_ProjectRast1.tif", prefix="SOIL", thematic=TRUE, proportion=TRUE);

269 polygons processed.

But can’t find the export or the results.


texts: Info from the getResultImageLayer missing

texts:
I was finally able to get the result layer overlaid on my application but now I am facing another problem.

The output of the layer url does not show any field or attribute details. I am not sure why this is happening? The output looks like this,

As you can see that it is only describing the Extent and other Map properties but the attribute details, which are visible in the jobs folder of the newly created final_shp, is missing on the url and so I am not able to fetch anything onto the map.

Any Idea how I can resolve this?


texts: Calculating number of points within polygon which are filtered according to a given rule

texts:
Description of my problem (somewhat simplified)

I have two datasets:

One polygon layer "villages" which represents the boundary of villages
One point layer "buildings" which represent buildings. This layer has the attributes "year of construction" and "type".

I would like to count the number of buildings (points) inside each village (polygone) which are of a given type and age, namely:

Number of individual houses with year of construction < 2000
Number of individual houses with year of construction > 2000
Number of appartment buildings with year of construction < 2000
Number of appartment buildings with year of construction > 2000

The following figure shows an example, with the blue polygon representing the village and each point a building, with styling applied to distinguish the 4 categories of buildings mentionned above (purple: buildings <2000; green: buildings >2000; triangles: individual houses; squares: appartment buildings):

I know how to count points within a polygon using "Analysis Tools" in the "Processing framework" plugin, but I do not how to count only the points which respect a given rule.

Any suggestions?


texts: Changing Buffer table geometry dynamically in PostGIS using plain SQL

texts:
I have two tables in PostGIS, a point and a point_buffer table. The point table has a buffer_distance field, with a default value say 200. Now I want to change the buffer table geometry whenever I change the buffer_distance value in my point table. I can do this for a single row in my point_buffer table using the following:

UPDATE point_buffer
SET the_geom = (SELECT ST_Buffer(the_geom,500) FROM point WHERE gid = 1)
FROM point
WHERE point.gid = point_buffer.gid

But whenever I am trying to change the whole of point_buffer table (droping WHERE clause in my sub-query) I am getting error message:

  'ERROR:  more than one row returned by a subquery used as an
  expression'.

My question is, can I change the entire point_buffer table in one go? I know one option is using a for lopp, with the upper bound being the count value of my point table and incrementing the point.gid value . But I want to do this in plain SQL.


texts: Calculate field using variables as fields gives ERROR 000539?

texts:
When I calculate a field like this:

field_0="field_name0"
arcpy.CalculateField_management("Layer", field_0, '!field_name1! / !field_name2!', "PYTHON_9.3")

It works! The field "field_name0" in the attribute table is populated with the value calculated from fields "field_name1" and "field_name2".

But, if I use variables instead of field names:

field_0="field_name0"
field_1="!field_name1!"
field_2="!field_name2!"
arcpy.CalculateField_management("Layer", field_0, 'field_1 / field_2', "PYTHON_9.3")

I get this error:

  Runtime error  Traceback (most recent call last):   File "",
  line 3, in    File "c:\program files
  (x86)\arcgis\desktop10.2\arcpy\arcpy\management.py", line 3354, in
  CalculateField     raise e ExecuteError: ERROR 000539: Error running
  expression: field_1/ field_2  Traceback (most recent call last):
  File "", line 1, in  NameError: name 'field_1' is
  not defined  Failed to execute (CalculateField).

I want to use loops later on, that's why I want it to work with variables, and not simply with direct field names.


texts: Determining the Geoid model used by GPS receiver

texts:
I have the following NMEA output by Holux M-241 GPS receiver:

$GPGGA,144014.000,2227.0987,N,11410.1382,E,2,6,1.67,32.6,M,-1.3,M,0000,0000*7B

I know that WGS84 height is higher than Geoid height 1.3 meter. But how can I know that which Geoid model it is using?

Also, is it a database stored in the GPS chip for Geoid height and WGS84 height conversion?


texts: Text size is too small in Mapbox style in R Leaflet map

texts:
Does anyone have any experience with mapbox basemaps in R leaflet?
The font is too small on my map. I have tried to change the zoomOffset and tileSize as detailed in this question, and in the tileOptions here, but no success.
library(leaflet)

leaflet_url <- "https://api.mapbox.com/styles/v1/williamlai/ckiy242l16cli19mvze8sqlpp/tiles/{z}/{x}/{y}?access_token={access_token}"

# mapbox basemap - font is too small
leaflet() %>%
  addTiles(urlTemplate = leaflet_url,
           tileOptions(tileSize = 512, zoomOffset = -1)) %>%  # the options don't seem to work
  addMarkers(lng=174.768, lat=-36.852)


texts: PyCharm doesn't fully recognise PyQGIS 3.18 API

texts:
I use PyCharm CE 2020.3.3 (also tested with 2020.2.5) and just upgraded QGIS to the new version 3.18 (first version with Python 3.9). When I updated my start script for PyCharm and updated the Python interpreter in PyCharm, the program recognises some parts of the API, but not all:

It seems that it's only finding those classes that are directly defined in the init.py of the qgis.core module.
Is that a QGIS bug or a PyCharm bug or related to how I loaded PyQGIS into PyCharm?
This is my start script:
@echo off
call C:\OSGeo4W64\bin\o4w_env.bat
call "%OSGEO4W_ROOT%\apps\grass\grass78\etc\env.bat"
@echo off
path %PATH%;%OSGEO4W_ROOT%\apps\qgis\bin
path %PATH%;%OSGEO4W_ROOT%\apps\grass\grass78\lib
path %PATH%;%OSGEO4W_ROOT%\apps\Qt5\bin
path %PATH%;%OSGEO4W_ROOT%\apps\Python39\Scripts

set PYTHONPATH=%PYTHONPATH%;%OSGEO4W_ROOT%\apps\qgis\python
set PYTHONPATH=%PYTHONPATH%;C:\Users\jludwig\PycharmProjects\QGIS-LRS-Rail;%OSGEO4W_ROOT%\apps\Python39\Scripts
set QGIS_PLUGINPATH=%QGIS_PLUGINPATH%;C:\Users\jludwig\PycharmProjects\QGIS-LRS-Rail;%OSGEO4W_ROOT%\apps\qgis\python\plugins
set PYTHONPATH=%PYTHONPATH%;%QGIS_PLUGINPATH%
set PYTHONHOME=%OSGEO4W_ROOT%\apps\Python39

path %PATH%;C:\Program Files\Git\bin

start "PyCharm aware of QGIS 3" /B "C:\Program Files\JetBrains\PyCharm Community Edition 2020.2.5\bin\pycharm64.exe" %*

And this is the list of paths known to the Python interpreter:


texts: How to set field attribute of selected features to NULL with PyQGIS 3

texts:
In my standalone QGIS script I'm selecting a couple of polygons and want to set its value in the field 'distance' to NULL.

I checked the PyQGIS Documentation, past SE Questions and came up with the following code:

for feature in layer.selectedFeatures():
    feature.setAttribute('distance', '')
    layer.updateFeature(feature)

return layer

The selection works perfectly fine but it doesn't change the attributes. According to the docs the function setAttribute() takes the field name as first variable and than the input you like it to have.

I also tried to use None instead of ' ', but it didn't work either.

What am I doing wrong?

I'm using QGIS 3.10


texts: Better ways to check if GRASS r.water.outlet works out fine output

texts:
I use my python script to call GRASS module r.water.outlet from outside, and here's my script(I perform this script in pyscripter):

import os
import sys

#set up GRASS environment variables
sys.path.append(os.path.join(os.environ['GISBASE'], 'etc', 'python'))
import grass.script as g
import grass.script.setup as gsetup
gisbase = os.environ['GISBASE']
gisdb = 'C:\Users\Heinz\Documents\grassdata'
location = 'newLocation'
mapset = 'TC'
gsetup.init(gisbase, gisdb, location, mapset)

#csv file reading and importation
import csv
rows = list(open('tc_sta.csv'))
totalrows = len(rows) - 1

#loop through the csv(coordinates) file in r.water.outlet module
f = open('tc_sta.csv', 'r')
element = list(csv.reader(f))
i = 0
j = 0
while True:
    if i <= totalrows:
        g.run_command('r.water.outlet', drainage = 'dra', basin = 'b' + str(i + 1), east = element[i][j], north = element[i][j + 1])
        i = i + 1
    else:
        break

#print all files in grass database to check results of the module
print g.read_command('g.list', _type='rast')

and here's the output in console,

The output shows that the GRASS module successfully runs and generates 8 new basin area raster files(b1 to b8), but I still want to check if there are something wrong. So I launch GRASS GUI (GRASS command line can't display map under Windows platform), using d.rast command to get visual of those output maps.

But I can only see b4 in display window. Others' are blank in display window. I 've checked the attribute of every output map using r.category, and they all have the same value as b4's.

I display b4 and input coordinates of module r.water.outlet(thus, batch point) together to check if other coordinates except b4's are out of range of basin area, but it shows not,

I am wondering that are there better ways to check this？(Seems that I can't use python from outside to display maps.) Because somehow maybe the outputs are all right.

I work with GRASS 6.4.3 and python 2.7.6 under Windows 8.1 64-bit.


texts: How to Change Vector Datasource in QGIS 3.2.0?

texts:
QGIS 2.18 has a very good function Change Vector Datasource:

In QGIS 3.2.0 there is Update SQL layer but it doesn't work like Change Vector Datasource. When I open it, select and do update, it acts in the way shown below on screenshots.

Please help how to change vector datasource as easy as in 2.18.


texts: Displaying direction of water system with arrows in QGIS?

texts:
I would like to put direction arrows next to my Water system just like in OpenTopoMap (see the image under).

I already watched the source code but they made it with Mapnik. However, I am working on my maps in QGIS.

In QGIS there are arrows line. I can put it next to my flow, but they are stick one after the another. I would like to have separated arrows by a certain distance.


texts: ModuleNotFoundError: No module named 'PyQt5.Qsci'

texts:
When I run my script I get this error:

from qgis.core import *
from qgis.gui import *
from qgis.utils import *
from PyQt5.QtCore import *
from PyQt5.QtGui import *

app = QgsApplication([], True)
path = "C:/OSGeo4W64/apps/qgis-dev"
app.setPrefixPath(path, True)
app.initQgis()

canvas = QgsMapCanvas()
title = "PyQGIS Standalone Application Example"
canvas.setWindowTitle(title)
canvas.setCanvasColor(Qt.white)

layer_info = 'LineString?crs=epsg:4326'
layer = QgsVectorLayer(layer_info, 'MyLine', "memory")
pr = layer.dataProvider()
linstr = QgsFeature()
wkt = "LINESTRING (1 1, 10 15, 40 35)"
geom = QgsGeometry.fromWkt(wkt)
linstr.setGeometry(geom)
pr.addFeatures([linstr])
layer.updateExtents()

QgsProject.instance().addMapLayer(layer)
canvas.setExtent(layer.extent())
canvas.setLayers([layer])
canvas.zoomToFullExtent()
canvas.show()
exitcode = app.exec_()
# QgsApplication.exitQgis()
sys.exit(exitcode)

Python 3.6 -
PyQt5 system versions.


texts: How to close/remove QGIS welcome page (recent projects) with Python?

texts:
Is it possible to close the QGIS welcome page with Python?

I have a permanently docked Python plugin, that uses showRubberBand to directly write into map canvas. Therefor I have to close/remove the welcome page first.


texts: Georeferencing a jpeg in QGIS but the scale is way off

texts:
I've gone through the really helpful walk-through(s) of how to georeference a jpeg in order to make it a layer, but after inputting the coordinates (either by typing in actual coordinates or clicking "From map canvas" and aligning points) and loading the new layer, the scale is way off (and I think the jpeg isn't even remotely in the right place--hard to tell because of how off the scale is).

Any ideas?


texts: Issue while using "Zoom to Rectangle" in Basemap Matplotlib

texts:
I am using the following code to plot the Basemap map on entire tkinter canvas. The code is as under:

from tkinter import *
from tkinter.ttk import *

import numpy as np
from matplotlib.figure import Figure

import matplotlib.pyplot as plt
from matplotlib.backends.backend_tkagg import (NavigationToolbar2Tk)
from matplotlib.backends.backend_tkagg import FigureCanvasTkAgg as FigureCanvas

import os

os.environ['PROJ_LIB'] = r'C:\Users\AppData\Local\Continuum\anaconda3\pkgs\proj4-4.9.3-vc14_5\Library\share'

from mpl_toolkits.basemap import Basemap

w=Tk()

fig_plot = plt.Figure(figsize=(7.5,4))
ax = fig_plot.add_subplot(111)

map = Basemap(projection='cyl',
    llcrnrlon = -180, llcrnrlat = -90,
    urcrnrlon = 180, urcrnrlat = 90,
            resolution='l', ax=ax)

map.drawmapboundary(fill_color='aqua')
map.fillcontinents(color='#ddaa66',lake_color='aqua')
map.drawcoastlines()

map.drawparallels(range(-90, 90, 30), labels = [False, True, True, False])
map.drawmeridians(range(-180, 180, 30), labels = [True, False, False, True])

Canvas = FigureCanvas(fig_plot, w)

ax.set_aspect('auto')

fig_plot.subplots_adjust(left=0.0, bottom=0.0, right=1.0, top=1.0)

Canvas.draw()

Canvas.get_tk_widget().pack(side='bottom', fill='both', padx=100, pady=10, expand = 1)

toolbar = NavigationToolbar2Tk(Canvas, w)
toolbar.update()
Canvas._tkcanvas.pack(side='bottom', fill='both', padx=10, pady=10, expand = 1)

w.mainloop()

With the above code I am successfully able to plot the basemap on entire tkinter canvas area as shown in the figure below.

Issue which I am facing is that when I am zooming to a particular section/region of the map by using Zoom to rectangle feature of matplotlib as shown in the figure below, then the basemap plot is getting stretched which I don't want.

Can someone help me out as to why this image stretching is happening and how to avoid it. Any help would be appreciated.
EDIT:
I just tried the following code in which if I am using Zoom to Rectangle in matplotlib, then the image is not getting stretched.
from tkinter import *
from tkinter.ttk import *

import numpy as np

import matplotlib.pyplot as plt
from matplotlib.backends.backend_tkagg import (NavigationToolbar2Tk)
from matplotlib.backends.backend_tkagg import FigureCanvasTkAgg as FigureCanvas

from matplotlib.figure import Figure

import os

os.environ['PROJ_LIB'] = r'C:\Users\AppData\Local\Continuum\anaconda3\pkgs\proj4-4.9.3-vc14_5\Library\share'

from mpl_toolkits.basemap import Basemap

w=Tk()

fig_plot = plt.Figure(figsize=(7.5,4))
ax = fig_plot.add_subplot(111)

map= Basemap(projection='cyl',
    llcrnrlon = -180, llcrnrlat = -90,
    urcrnrlon = 180, urcrnrlat = 90,
            resolution='l', ax=ax)

map.drawmapboundary(fill_color='aqua')
map.fillcontinents(color='#ddaa66',lake_color='aqua')
map.drawcoastlines()

map.drawparallels(range(-90, 90, 20), labels = [False, True, True, False])
map.drawmeridians(range(-180, 180, 20), labels = [True, False, False, True])

C22 = FigureCanvas(fig_plot, w)

ax.axis('equal')

fig_plot.subplots_adjust(left=0.0, bottom=0.0, right=1.0, top=1.0)

C22.draw()

C22.get_tk_widget().pack(side='bottom', fill='both', padx=100, pady=10, expand = 1)

toolbar = NavigationToolbar2Tk(C22, w)
toolbar.update()
C22._tkcanvas.pack(side='bottom', fill='both', padx=10, pady=10, expand = 1)

w.mainloop()

After running the above code I am getting the following tkinter screen:

But now the issue is that the basemap map is not occupying the entire plot area and there are blue blank spaces between the Meridian & Parallel Lines and the plot area as shown in the figure above.


texts: Using ST_AsMVT with ST_TileEnvelope clipping in mapbox

texts:
I have been searching for a while for a way to correctly project the ST_AsMVT to ensure it fits to the map. The below query was the only suggestion I found that will correctly transform to the mapbox EPSG:3857 coord system.

However, I have noticed at the global extent zoom levels the bounding box is clipping it seems. Just wondered if there is anything obvious that I have missed here given that other than this it projects perfectly. The screenshot is using an example with natural earth country polygons imported into postgis.

I also noticed that sometimes there is a square cut out the middle too which when researching looked like a smilier issue with the way different systems cut the grid of the earth?

I have also seen a few Postgis BBox functions but none of them work like the sql below.

Any advice would be greatly appreciated.

    WITH webmercator(envelope) AS (
      SELECT ST_TileEnvelope(${z}, ${x}, ${y})
    ),
    wgs84(envelope) AS (
      SELECT ST_Transform((SELECT envelope FROM webmercator), 4326)
    ),
    geometries(wkb_geometry) AS (
      SELECT ST_Transform(wkb_geometry, 3857)
      FROM ${table}
      WHERE wkb_geometry && (SELECT envelope FROM wgs84)
    )
    SELECT ST_AsMVT(tile) as mvt FROM (
      SELECT  ST_AsMVTGeom(wkb_geometry, (SELECT envelope FROM webmercator))
      FROM geometries
    ) AS tile

Zooming and moving viewport, reveals the polygon (Just to show it does have all the data)


texts: sentinelsat Python API error in download

texts:
I am having a go at using the sentinelsat python API to download satellite imagery. However, I am receiving error messages when I try to convert to a pandas dataframe. This code works and downloads my requested sentinel satellite images:

from sentinelsat import SentinelAPI, read_geojson, geojson_to_wkt
from datetime import date

api = SentinelAPI('*****', '*****', 'https://scihub.copernicus.eu/dhus')

footprint = geojson_to_wkt(read_geojson('testAPIpoly.geojson'))

products = api.query(footprint, cloudcoverpercentage = (0,10))

#this works
api.download_all(products)

However if I instead attempt to convert to a pandas dataframe

#api.download_all(products)

#this does not work
products_df = api.to_dataframe(products)

api.download_all(products_df)

I receive an extensive error message that includes

  "sentinelsat.sentinel.SentinelAPIError: HTTP status 500 Internal Server Error: InvalidKeyException : Invalid key (title) to access Products
  "

(where (title) is also replaced with processed, platformname, processingbaseline, etc.). I've tried a few different ways to convert to a dataframe and filter/sort results and have received an error message every time (note: I have pandas/geopandas installed). How can I convert to a dataframe and filter/sort with the sentinelsat API?


texts: ArcObjects: How to Speedup Spatial Query to Feature Layer?

texts:
I have ROAD (POLYLINE) Shape File with more than 400,000 Records. I am working on a code which goes like this:

Loop through all Feature of the Shape File
Take First Feature
Find the LAST VERTEX of that feature. I get a POINT
Query to FEATURE Layer all the POLYLINES Connected to that POINT and get back result in feature cusros
Do some Process on returned feature of RESULT Cursor
Go to Step1 until Loop is not completed.

When I do STEP-4, it obviously take lot of time because, I create a BUFFER of Say 1 meter on that POINT and Make a spatial FILTER and Run LAYER.SEARCH QUERY on entire featurelayer using that SPATIAL FILTER. Everytime querying 400,000 records making it very slow.

Is there a better way to do this thing ?

I have to use ArcObjects VBA.
I do not want to Create NETWORK from
the layer.


texts:
I'm attempting to use gdal_merge to stitch together four JP2 aerial orthophotos. the merge is successful, but the resulting image has washed out color, as you can see here:

All color in the four original tiles looks fine before the merge.

The line I'm using to merge them is:

gdal_merge -o merged.tif -of GTiff 16SEJ790950_201103_0x1500m_4B_1.jp2 16SEJ790965_201103_0x1500m_4B_1.jp2 16SEJ805950_201103_0x1500m_4B_1.jp2 16SEJ805965_201103_0x1500m_4B_1.jp2

GDALINFO color information shows:

COLORSPACE = RGB

COMPRESSION_RATE_TARGET=2

Band 1 Block=256x256 Type=Byte, ColorInterp=Red
  Description = Red
  Overviews: 5000x5000, 2500x2500, 1250x1250, 625x625, 312x312, 156x156

Band 2 Block=256x256 Type=Byte, ColorInterp=Green
Description = Green
  Overviews: 5000x5000, 2500x2500, 1250x1250, 625x625, 312x312, 156x156

Band 3 Block=256x256 Type=Byte, ColorInterp=Blue
  Description = Blue
  Overviews: 5000x5000, 2500x2500, 1250x1250, 625x625, 312x312, 156x156

Band 4 Block=256x256 Type=Byte, ColorInterp=Undefined
  Description = AllOpacityPremultiplied
  Overviews: 5000x5000, 2500x2500, 1250x1250, 625x625, 312x312, 156x156


texts: Is there a way in ArcObjects to restrict attribute value when user modifies it. Let's say for integer value between 0-10 or for string : no/yes?


texts: Getting error while building Virtual Raster

texts: The following layers were not correctly generated.C:/Users/Shreelakshmi C R/AppData/Local/Temp/processing_45117873c29b4e15a154e8405acddbfb/99d397046dcf4d578d7fc4ade40e9ba8/OUTPUT.vrtYou can check the 'Log Messages Panel' in QGIS main window to find more information about the execution of the algorithm.


texts: Reduce raster resolution based on cell ratio

texts:
I have 2m impervious/non-impervious data. I would like to reduce the resolution to 30m and have the cell value for the 30m data be determined by the ratio or proportion of impervious / non-impervious 2m cells that will fall within the 30m block.

I am using ArcGIS 10.2.1
The diagram in red is what I am trying to achieve .
I= impervious
N=Non-impervious


texts: Automating ArcGIS Online Shortlist Story Map description customization?

texts:
Is it possible to automate the descriptions under each image/event in the Shortlist Story Map?

I want the word "Year:" to appear before every yyyy.

I want this type of label and colon to be bold and appear to the left of every item in my description box.

Is there a way to customize this without having to edit by hand every single one of my 45 events?

I'm currently using a template, but I'm open to imputing some type of code if that's what's needed. I have never made a custom-code for a Story Map before.

Link to the Shortlist https://azdema.maps.arcgis.com/apps/Shortlist/index.html?appid=df7dc1d3908545b58188f8939dd7b115