texts:
I calculated the EVI for a Landsat image from -1 to 1. However, when I use the inspector tool in GEE and click around on my image, I see I have values less than -1. Is there something wrong with my code?
var collectionID = 'LANDSAT/LC08/C01/T1_TOA'
var startDate = '2018-06-01'
var endDate = '2019-07-01'

var data = ee.ImageCollection(collectionID)
          .filterBounds(area)
          .filterDate(startDate, endDate)
          .filterMetadata('CLOUD_COVER', 'less_than', 10)
          .median();
print(data)

var FCC ={ bands: ['B5','B4','B3'],min:0, max: 0.3};
Map.addLayer(data,FCC,'Landsat8');

var evi = data.expression(
  '2.5 * ((NIR - RED) / (NIR + 6 * RED - 7.5 * BLUE + 1))', {
    'NIR' : data.select('B5'),
    'RED' : data.select('B4'),
    'BLUE' : data.select('B2')
});
print(evi)

Map.addLayer(evi, {min: -1, max: 1, palette: ['red', 'brown', 'yellow', 'green']}, 'EVI');


texts: Computing NDVI in ArcMap from Landsat 8 "Analysis-Ready" Surface and TOA Reflectance Products

texts:
The internet is full of various contradictory advice on corrections one can make to get NDVI computed from Landsat 8 to fall between -1 and 1.

This is for instructional purposes, so I hope to use the Analysis-Ready Data products (surface and TOA reflectance), which have already been atmosphere- and sun-angle-corrected.

I assume the surface reflectance product would be the best for NDVI, but the values I get (raster calculator: 1.0*("NIR"-"red")/("NIR"+"red")) range from -1.14 to -0.6, which seems unrealistic.  The TOA product gives slightly better numbers, -1.1 to 1.38, but I still don't know why we're exceeding the -1 to 1 range.  Any ideas?

The values on the SR product range from -522 to 10133 (red) and -47 to 10558 (NIR) and on the TOA product they are more positive: 216 to 9775 (red)and 66 to 10600 (NIR).


texts: Issues with creating an address locator using multi-point features in Arcmap/catalog 10.4?

texts:
The error I get when trying to create this locator is "000042: Failed to create the address locator which states that I must -

have permission to write to the disk/database
have enough disk space
have a valid connection
have valid geocoding rule files
have valid address locator style

I have write permissions, I have enough disk space (1+ tb), and I do not require a connection which leaves #4 and #5 as possible reasons why it fails.

I am trying to use the gazetteer locator style as each of the features represent a geographic place, and a multi-point feature (which I suspect is the reason for the error).

So in short, does Arcmap/catalog 10.4 have issues with creating an address locator using multi-point features? Is there a work around for this? Is the gazetteer locator style the wrong style to use in this case?


texts:
I'm reading into how to customise projections with Proj4 - specifically to achieve tilted projections such as the example below (from vis4.net)

The Proj4 literature is pretty heavy going and having explored it I'm not clear which parameters I'd need to change, or whether a custom datum would need to be defined.


texts: Importing OpenStreetMap polygons to QGis with relation attributes

texts:
I'am trying to import large OSM data to QGis, namely all buildings polygons in Czech Republic.
Finally, I imported polygons using PostGIS layer, because everything else had strong issues. See my previous question: Importing large OSM data to QGis failed (using QuickOSM, SpatialLite layer or add vector layer).

Now, I have another problem: some attributes (tags) are missing - probably the attributes, that are placed in relations.

Is there method, how to apply OSM relation tags to it's sub polygons and import them to QGis?

Example: this building (relation): openstreetmap.org/relation/3679169 has interesting tag "building:levels" but PostGIS layer in QGis contains only  polygon without it.


texts: European Rail Freight Map (shapefile or other ArcGIS format)?

texts:
I need European Rail Freight Map (shapefile or other ArcGIS format) for a GIS project.

Does anyone have this data?


texts: Remove additional lines from vectorized stream network

texts:
I have a large raster stream network (excerpt below) in QGIS that has been thinned using r.thin (light blue squares) and then vectorised using r.to.vect (pink lines). However, because some of the streams are close together, there are extra unwanted lines (in yellow) that are created as part of the vectorisation process.

Is there a process for removing these yellow lines without affecting the rest of the network?
Ideally there might be a process that only involves using the vector layer itself. However, if that's not possible, it could involve the use of one of the underlying raster layers that was used to derive the stream network (eg hydrologically controlled DEM or D8 flow accumulation raster).
EDIT: the thinning process removes some cells from the raster that are part of the drainage channels (think of a T-intersection for example - the head of the T is removed, making the intersection a Y). So using the DEM to determine line direction is not accurate, because the "downhill" end of the line may be higher than the "uphill" end.


texts: How to sort a list with dictionaries into ascending order

texts:
I used reduceRegion to get statistics of certain provinces in Germany. I get a list with dictionaries (province ID, histogram and median) and would like to sort the list by the median values. I have issues to get the list sorted.  This is what I have so far:

var geometry = ee.Geometry.Polygon(
        [[[2.165739166604279, 55.7242449869763],
          [2.165739166604279, 45.702626323936244],
          [19.56808291660428, 45.702626323936244],
          [19.56808291660428, 55.7242449869763]]], null, false),
    DEU = ee.FeatureCollection("users/gaertnerp/DEU_gadm1");

var NO2 = ee.ImageCollection('COPERNICUS/S5P/NRTI/L3_NO2')
.filterBounds(geometry).select('NO2_column_number_density')
.map(function(img) {return img.multiply(1e6).copyProperties(img,['system:time_start']); })
.filterDate('2020-01-01', '2020-02-15').median().rename('NO2');

var empty = ee.Image().byte();

var outline = empty.paint({
  featureCollection: DEU,
  color: 'ID1'
});

var reducers = ee.Reducer.histogram().combine({
  reducer2: ee.Reducer.median().unweighted(),
  sharedInputs: true
});

var hist = NO2.addBands(outline).reduceRegion({
  reducer: reducers.repeat(1).group({
    groupField: 1,
    groupName: 'NAME_1'
  }),
  geometry: geometry,
  scale: 5000,
  maxPixels: 1e15
})

var groups = ee.List(hist.get('groups'))
groups = groups.slice(0,5)
print('groups',groups)

Here is the list I would like to sort by the median value.


texts: FeatureClassToFeatureClass_conversion in ArcGis 10.1

texts:
I have many > 60 .dwg files I want to convert just the polyline parts of these and have put together a list of python commands using excel to get something like this:

However when I try running this in the python window or from a saved script I keep getting input features and output location do not exist.

What am I doing wrong?


texts: Fixing Segmentation fault: 11 from PyQGIS

texts:
I am trying to run the default script that PyQGIS Cookbook provides along with adding few lines to read my project. It tries to read the read project but at the end gives a Segmentation Fault. How can I fix this?

import sys
sys.path.append('/Applications/QGIS3.app/Contents/Resources/python')

from qgis.core import *

# supply path to qgis install location
QgsApplication.setPrefixPath("/Applications/QGIS3.app/Contents/MacOS", True)

# create a reference to the QgsApplication, setting the
# second argument to False disables the GUI
qgs = QgsApplication([], False)

# load providers
qgs.initQgis()

project = QgsProject.instance()
project.read('polygondata.qgs')
print(project.fileName())
print(project.mapLayers())

# Write your code here to load some layers, use processing
# algorithms, etc.

# When your script is complete, call exitQgis() to remove the
# provider and layer registries from memory

qgs.exitQgis()


texts: Building spatialindex in Postgres with JNDI for GeoNetwork?

texts:
Which steps should be taken to build the spatialindex table inside PostgreSQL when using a JNDI connection?

Is the table built by GeoNetwork or does it need to be created manually?

What I have tried so far:

For building the spatialindex table in GeoNetwork (3.0.6) in the Postgres database as opposed to using a shapefile internally, I switched database connection to JNDI, as recommended in the advanced configuration docs.

Creating JNDI connection first

To do so, I had to change the database type used in srv.xml to <import resource="../config-db/jndi-postgres-postgis.xml"/> and configure my datasource in context.xml like so:

<Resource
        name="jdbc/geonetwork"
        type="javax.sql.DataSource"
        auth="Container"
        driverClassName="org.postgresql.Driver"
        url="jdbc:postgresql://localhost:5432/geonetwork"
        username="geonetwork_user"
        password=""
        maxActive="20"
        maxIdle="20"
        maxWait="200"
        initialSize="20"
        removeAbandoned="true"
        removeAbandonedTimeout="120"
        logAbandoned="true"
        testOnBorrow="true"
        defaultReadOnly="false"
        defaultAutoCommit="false"
        validationQuery="SELECT 1"
        timeBetweenEvictionRunsMillis="10000"
        minEvictableIdleTimeMillis="1800000"
        testWhileIdle="true"
        numTestsPerEvictionRun="3"
        poolPreparedStatements="true"
        maxOpenPreparedStatements="1200"
        provideDataStore="true"  --> tried this to build spatialindex table
/>

I added the provideDataStore parameter according to this section. Actually it looks like it must go into the /config-db/jndi-postgres-postgis.xml file, but this is rejected by the application.

I adapted some configuration in the /config-db/jndi-postgres-postgis.xml file as mentioned here:

 <bean id="jdbcDataSource" class="org.springframework.jndi.JndiObjectFactoryBean">
        <property name="jndiName" value="java:/comp/env/jdbc/geonetwork"/>
        <property name="cache" value="true"/>
        <property name="exposeAccessContext" value="false"/>
 </bean>
 <bean id="jpaVendorAdapterDatabaseParam" class="java.lang.String">
        <constructor-arg value="POSTGRESQL"/>
 </bean>

With this configuration, I was able to switch the database connection to JNDI.

Build spatialindex into Postgres

But nevertheless, it does not build the spatialindex into Postgres.

I tried:

Restarting tomcat
Harvesting new datasets
Playing around with the provideDataStore parameter
Creating the table spatialindex myself to make GeoNetwork use it on restart

The table is defined according to the definition from the 2.10 branch's create-db-postgis.sql file:

CREATE TABLE spatialindex
 (
    fid int not null,
    id  varchar(250),
    primary key(fid)
  );
CREATE INDEX spatialindexNDX1 ON spatialindex(id);
SELECT AddGeometryColumn('spatialindex', 'the_geom', 4326, 'MULTIPOLYGON', 2 );
CREATE INDEX spatialindexNDX2 on spatialindex USING GIST(the_geom);

But my spatialindex.shp file in the data directory still gets altered instead of the database table.

Is there anything else neccessary to make GeoNetwork build the spatialindex table in Postgres?


texts: QLR loaded via PyQGIS not loading correctly but works when using layer--add qlr

texts:
Based on Load .qlr with loadNamedStyle I was able to add a QLR to a python tool via PyQGIS -however, it appears broken...

            elif ftype=='qlr':
                try:
                    print ('Loading QLR layer')
                    #ref:https://gis.stackexchange.com/questions/240691/load-qlr-with-loadnamedstyle/240694?noredirect=1#comment511181_240694
                    print("QLR %s called and loading into group %s" % (filePath, fname))
                    root=QgsProject.instance().layerTreeRoot()
                    print('adding group')
                    group=root.addGroup(fname)
                    QgsLayerDefinition().loadLayerDefinition(filePath,QgsProject.instance(),group)
                except:
                    print ('Error loading QLR')

fname is the name of the group and file, filepath is the path/file.qlr

The QLR is loading a shp that's in the same folder.

It looks like

 <!DOCTYPE qgis-layer-definition>
<qlr>
  <layer-tree-group checked="Qt::Checked" name="" expanded="1">
    <customproperties/>
    <layer-tree-layer id="events_6d4a5418_1dce_4ba0_201903071250415994831" checked="Qt::Checked" name="events" source="./events.shp" providerKey="ogr" expanded="1">
      <customproperties/>
    </layer-tree-layer>
  </layer-tree-group>
  <maplayers>
    <maplayer autoRefreshEnabled="0" labelsEnabled="1" refreshOnNotifyMessage="" maxScale="0" minScale="1e+08" hasScaleBasedVisibilityFlag="0" readOnly="0" styleCategories="AllStyleCategories" autoRefreshTime="0" type="vector" simplifyDrawingTol="1" refreshOnNotifyEnabled="0" simplifyDrawingHints="0" simplifyAlgorithm="0" simplifyLocal="1" simplifyMaxScale="1" geometry="Point">
      <extent>
        <xmin>142.5193149000000119</xmin>
        <ymin>-37.29052469999999886</ymin>
        <xmax>143.2661899999999946</xmax>
        <ymax>-36.60563100000000247</ymax>
      </extent>
      <id>events_6d4a5418_1dce_4ba0_201903071250415994831</id>
      <datasource>./events.shp</datasource>
      <keywordList>
        <value></value>
      </keywordList>
      <layername>events</layername>

and it comes out like

However...

It works fine if I load it in the QGIS python console using the same commands...

Note screen grabs are of a complex and simple QLR I tried but both end up the same...


texts: Where are the Expression-based Labeling options in QGIS?

texts:
I have imported a layer into my map to show the location of towns/cities. However I had to set graduated styling so that only cities over a certain population show.

I now want to label the ones that show but if I add labels all of the labels show rather than just the ones which I want to appear in my map.

I think I can solve the problem using Expression-based labeling. Following this tutorial: https://underdark.wordpress.com/2011/10/27/expression-based-labeling-for-qgis/

But I can't find the settings. When I go into the labelling options I see this:

What do I do next? Any help would be amazing.


texts: Why does Select by Location select points outside search radius?

texts:
I'm using Select by Location to find wind turbines within 5-km of a center point feature -- why does this command also select nearby turbines?

The image below shows my selected points with a 5-km buffer circle drawn for comparison. As you can see, turbines that are not even close to touching the circle are included.


texts: How do I plot points as graduated/proportional circles in R?

texts:
I'm back again with more R mapping questions! This one is concerning symbolizing points by a certain attribute with graduated/proportional circles.

DATA: My full code can be found here if you're interested. My CSV of interest is here. Ethiopia shapefile can be downloaded here

I've plotted points from the CSV over a shapefile. For the sake of brevity, here's the heavily abridged version of my code:

library(raster)
library(rgdal)

#set your working directory

#read in eth shapefile
eth <- readOGR(dsn = "D:/Mapping-R/Ethiopia", layer = "ETH_adm0")

#read in total branch location csv
branches <- read.csv("Branches_Africa.csv", header = TRUE)

#transform branches data frame to spdf for mapping
coordinates(branches) <- ~Lon + Lat

#plot shapefile with point overlay
plot(eth)
points(branches$Lon, branches$Lat, col = "blue", pch = 16, cex = .5)

As you can see in the branches CSV, there are several different attributes. In my case, I would like to symbolize the "share" attribute with graduated/proportional circles, as seen here:

Ideally I would also like to change the colors of each point as well to match the QGIS symbology as closely as possible.

If you look at my full code, you'll notice that I've not used ggplot to map, so I would prefer the solution not come from that package (or would ggplot make things easier once I've transformed all the data to fit the package?).  I tried using spplot as per this post with:

spplot(branches, "share", col.regions = brewer.pal(9, "Reds"), cuts = 8, scales = list(draw = T))

but the code gave me the following error: Error in fill.call.groups(args.xyplot, z = z, edge.col = edge.col, colorkey = colorkey,  : number of colors smaller than number of factor levels

I also tried to work with mapCircles from the rCarto package, but I had no luck there either.

SO, how can I change the symbology of branches$share so that graduated/proportional (and preferably differently colored) circles are displayed instead of just one size/color?


texts: Tiling Rasters based on an existing Polygon Shapefile (Vector Grid)

texts:
I have a tile index representing the tile scheme used to tile Lidar data.  I want to tile imagery rasters using the same tiling scheme.  I have not found anyway using open source tools to tile a raster based on a shapefile template.  All the tools I have found require you to know the tile size and origins.

Are there any tools (especially GDAL) that will do this out of the box?


texts:
I am trying to control how many road symbols and how do I do that ? I have geodatabase that has Transportation and in it has a featureclass, "Major_Rds" and it also has three symbology for Interstate, NM state roads, and US roads.

What I would like to do is to limit a few symbols to be display at a specific scale.

You can see there are too many road symbols in Southeast New Mexico but you see a few in Northwest Texas.


texts: Dropping a table from PostGIS database using pgAdmin III

texts:
A table in my PostGIS database (PostgreSQL 9.3.5) had the wrong geometry-type, therefore I wanted to change the geometry type.
I tried to do this by dropping the geometry column like this:

ALTER TABLE xxx.my_table
 DROP COLUMN geometry_column;

After which I would have added a new geometry column of the correct type.
In hindsight, I should probably have used ALTER COLUMN, however this didn't occur to me at the time.

The query started to run, however it never finished, so I had to stop the execution.
Following this, every time I click on the table from within pgAdmin, the program crashes. I have tried using DROP TABLE to delete the table, but the query just runs without finishing.
I also tried deleting the entire schema from within pgAdmin, however this also causes a crash.

The program gets stuck when "retrieving details" on the table, so I suspect the table is corrupt in some way. We have a backup of the schema in a working state, however using the restore functionality in pgAdmin doest not seem to work either.

Is there anyway to "fix" the table or force delete it in any way?


texts: Display XY data not displaying for some rows in a CSV, but displaying for the rest

texts:
I am using ArcMap 10.7 to map data points (decimal degrees) across eight countries. I have an Excel .csv file that contains the data points (548 rows), but when I import the data (as a .csv) and try to "Display XY data", some of the points do not display. Note: all of these points are in one country. Within Excel, I double-checked that the latitude and longitude values for this country are correct (they are). I have plotted these same data coordinates (decimal degrees) in a smaller excel .csv before, and they displayed in ArcMap.
Here are the steps I followed to import my data:
Option 1: Catalog -> Navigate to folder containing .csv file -> Drag .csv file into Data frame -> Use "Make XY Event Layer" from the ArcToolbox -> Select Latitude (Y) and Longitude (X) -> Add to map
Option 2: File -> Add Data -> Add XY Data -> Select file -> Select Latitude (Y) and Longitude (X) -> Add to map
Both of these scenarios did not display all data.
Problem:
Some locations not appearing in Display XY data
Solutions I've tried:

Import dataset as a .csv file into an empty ArcMap document (i.e. no basemap or other shapefiles present).
Re-save Excel .csv file as a new .csv, re-import into ArcMap
Clear formatting of .csv, re-save as a new .csv, re-import into ArcMap
Delete all other rows in .csv besides Latitude, Longitude and "Country", then re-save .csv and re-import into ArcMap
Click .csv file in Table of Contents, click Data -> Export -> Save as dBase table -> Add to map as layer
Open up old .csv where the data for this country appeared, copy these data into a new .csv with the rest of the coordinates that work, re-import into ArcMap.

None of these solutions have worked. Is there anything else I can do to fix this issue?


texts:
I am working on a one-band raster with pixel values from 1 to 7, each representing a land-use class. I also have a grid of polygons and what I need is to obtain the percentage of each class in each cel of the grid.
I already used the Zonal Statistics tool and with rasters representing only one variable each it adds three columns to the attribute table of the grid: count of values, sum and average. If I do the same with the land-use raster, I only obtain the same count, sum and average values but only for the overall pixels with no distinction for the classes.


texts:
I set up Open Source Routing Machine (OSRM).  The documentation does a pretty good job of laying everything out and I was able to get a server set up and running on my machine.  My question is about handling the JSON output from my server.  As explained in the documentation the routing information is encoded.  My application was built accepting JSON from Mapquest's Open API which is unencoded.

My question is do I need to perform a special function to unencode the data?  Or better yet, is their some feature baked into OSRM I am not using properly to get the unencoded results?

edit

Solved: Adding ?compression=false to the request query will output unencoded data.


texts: I need to create a subset of data based off their distance from particular points in ArcMap. Is it possible to create Feature layers that have a defined buffer, that only show data within that buffer? Or must I create a physical subset of the data based my buffered selection?


texts: Exporting several (6 tiles) Landsat NDVI Google Earth Engine output from median images to my drive

texts:
I am new to Google Earth Engine. I managed to produce NDVI on the map window from 6 Landsat tiles from median values .median(). With the .first() I managed to export it but with the median() I get several tif files without data. What is the problem? Here is the code I used.

// var surfaceReflectanceL5 = ee.ImageCollection('LANDSAT/LT05/C01/T1_SR');
var img = ee.Image(ee.ImageCollection('LANDSAT/LC08/C01/T1_SR')
    .filterBounds(roi)
    .filterDate('2018-01-01', '2018-12-31')
    .sort('CLOUD_COVER')
    .median());
Map.addLayer(img, {bands: ['B5', 'B4','B3'],min:0, max: 3000}, 'True colour image');

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

// Make a palette: a list of hex strings.
var palette = ['FFFFFF', 'CE7E45', 'DF923D', 'F1B555', 'FCD163', '99B718',
               '74A901', '66A000', '529400', '3E8601', '207401', '056201',
               '004C00', '023B01', '012E01', '011D01', '011301'];

// Display the input image and the NDVI derived from it.
Map.addLayer(img.select(['B5', 'B4', 'B3']),
         {gain: [0.1, 0.1, 0.1]});
Map.addLayer(ndvi, {min: 0, max: 1, palette: palette}, 'NDVI');

// the problem is here: why it is exporting several layers without data (sizes less than 200kb)
// Export the image, specifying scale and region.
Export.image.toDrive({
  image: ndvi,
 description: 'NDVI_2018_A1',
 scale: 30,
 region: roi
 });


texts:
I'm building a web application, which involves (currently) a REST backend and a frontend SPA. The backend is hosting massive geo-enabled data, and the frontend is displaying it on a mapbox map. I'm hitting performance issues, because of the stack I set up. Here's the stack, when it comes to displaying a map.

    MapboxGl.js player
        ^          ^
        | (HTTP)    \
 Node.js backend     Mapbox APIs vectors
        ^
        |
     MongoDB

The HTTP endpoints serve both GeoJSON and REST data. The GeoJSON part is dynamically generated from a set of MongoDB collections (say, gpspoints, images, lines etc.). The Mapbox API is my Mapbox Studio style.

I'm looking for the best alternative to this design. And reading about maps, vectors, tiles, rendering and stuff, I came up to the idea that I should trade my GeoJSON server for a Vector tiles server (only for the map part). I'd like to end up with something like.

         MapboxGl.js player
       ^         ^           ^
      /          | (HTTP)     \
Vector tiles  Node.js backend   Mapbox APIs vectors
   Server            ^
     ^               |
     |               |
Other backend ?   MongoDB

I'm struggling at the point of making design choices.
The questions I have in mind are :

About generating :

I came to the idea that I must generate my vector tiles from a GeoJSON. Is that correct ? Would it be possible to generate vector tiles from scratch, I mean, say, from a node.js routine script, parsing datas from various sources (.csv, other databases, etc.) and make operations on it (using turf.js), then generating vector stuff, like pushing, maybe updating, deleting, etc ?
I found geojson-vt which seems to be able to translate GeoJSON to vector tiles, but only in json format. I found vt-pbf, am I on the right direction ?

About serving

Once I would have been able to generate whatever-format-vector-tile, I am going to be willing to serve it. I found Tilelive.js which seems very good, and has plenty of modules, including backend ones. Is that a good option ?
Speaking generaly, would it be possible not to generate any static file stuff, and maybe dynamically serve any-vector-format from a database, through something like tilelive ?

Additional info :

My datas are "mutable", but I have routine scripts that run at night and parse new datas. Thus, live-map-displayed-data mutability is not a requirement : I can update / re-parse a whole dataset at night to serve it on the morning
Both generating and hosting on Mapbox studio is not an option, since my data is critical and must be hosted in specific coutries. A self-hosted option is mandatory, any cloud-based stuff is not ok :/

This question is redundated in softwareengineering.stackexchange


texts: OpenStreetMap: need a tile service with English labels in non-English countries

texts:
I want to use OpenStreetMap tile server in its QtLocation plugin.

Alas, if I use the host provided by default, there are non-English labels in non-English countries (Arabic, for example).

I don't know Arabian, so I need a tile service where there are English labels. Can I get one?


texts: Move Python Add ins to ArcGIS Pro without having to learn a new programming language?

texts: I have created a Python Add in for ArcMap and it has been tested and functions very properly. However, I am looking online, and it seems like I have to go through a lot of processes just for it to function in ArcGIS pro. I need to be able to do this quickly but I just haven't seen anywhere that would allow this to be a smooth transition.


texts: Missing 'Z' component from points extracted from polyline in QGIS

texts:
I have polylines which were created with XYZ coordinates and after extracting the points (nodes) I want to export these to a .csv file with their XYZ coordinates. When I try the 'save as' function and specify the layer options 'AS_XYZ' I get the file with all z coordinates as '0'.
Can anyone help me get the missing component?
I know the height data is associated with the shapefiles as they display in 3D correctly, so why can i not retrieve it?

I am using QGIS las palmas on a mac. Also ArcGIS on PC (using XTools), which gives the same problem!


texts: How can I merge two bands with different sizes using Rioxarray?

texts:
I have two images that I want to combine as follows
Expected result:

I tried to use the function rioxarray.merge.merge_arrays() in the following way but I got this error that I can't understand:
>>> import rioxarray
>>> from rioxarray import merge
>>> b4_1 = rioxarray.open_rasterio("LC08_L1TP_204030_20200607_20200607_01_RT/LC08_L1TP_204030_20200607_20200607_01_RT_B4.TIF")
>>> b4_2 = rioxarray.open_rasterio("LC08_L1TP_204031_20200607_20200607_01_RT/LC08_L1TP_204031_20200607_20200607_01_RT_B4.TIF")
>>> merged = merge.merge_arrays([b4_1, b4_2], nodata=0.0)
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/home/ruben/anaconda3/envs/rs-aCoisa/lib/python3.7/site-packages/rioxarray/merge.py", line 113, in merge_arrays
    **{key: val for key, val in input_kwargs.items() if val is not None},
  File "/home/ruben/anaconda3/envs/rs-aCoisa/lib/python3.7/site-packages/rasterio/merge.py", line 250, in merge
    boundless=False, masked=True, indexes=indexes)
  File "/home/ruben/anaconda3/envs/rs-aCoisa/lib/python3.7/site-packages/rioxarray/merge.py", line 43, in read
    Image.fromarray(data_window[iii]).resize((width, height))
  File "/home/ruben/anaconda3/envs/rs-aCoisa/lib/python3.7/site-packages/PIL/Image.py", line 1905, in resize
    return self._new(self.im.resize(size, resample, box))
ValueError: image has wrong mode
>>> b4_1
<xarray.DataArray (band: 1, y: 7811, x: 7691)>
[60074401 values with dtype=uint16]
Coordinates:
  * band         (band) int64 1
  * y            (y) float64 4.899e+06 4.899e+06 ... 4.665e+06 4.665e+06
  * x            (x) float64 5.082e+05 5.082e+05 ... 7.389e+05 7.389e+05
    spatial_ref  int64 0
Attributes:
    STATISTICS_MAXIMUM:        63430
    STATISTICS_MEAN:           13305.927976843
    STATISTICS_MINIMUM:        0
    STATISTICS_STDDEV:         13401.114806381
    STATISTICS_VALID_PERCENT:  100
    transform:                 (30.0, 0.0, 508185.0, 0.0, -30.0, 4899015.0)
    scale_factor:              1.0
    add_offset:                0.0
    grid_mapping:              spatial_ref
>>> b4_2
<xarray.DataArray (band: 1, y: 7831, x: 7701)>
[60306531 values with dtype=uint16]
Coordinates:
  * band         (band) int64 1
  * y            (y) float64 4.741e+06 4.741e+06 ... 4.506e+06 4.506e+06
  * x            (x) float64 4.71e+05 4.71e+05 4.711e+05 ... 7.02e+05 7.02e+05
    spatial_ref  int64 0
Attributes:
    transform:     (30.0, 0.0, 470985.0, 0.0, -30.0, 4740615.0)
    scale_factor:  1.0
    add_offset:    0.0
    grid_mapping:  spatial_ref

What can I do to fix the problem? Or what alternative to I have to do this but still ending up with a xarray.DataArray?


texts:
Is there a default solution for choosing certain profile before GUI of QGIS app will load ? I want to be able to choose one profile from set of profiles always before GUI loads.
Is there a way to disable GUI in QGIS app before receiving master password for certain profile ?


texts:
We just moved from 10.0 to 10.2 and woke up to the wonderful news that all our geolocators are broke and need to be recreated. Most were made in 9.3. This is fine, but the locator building process has changed significantly and to be honest I've inherited the GIS coordinator/admin role with very little knowledge of ESRI products.

I'm trying to build a locator or composite locator that will check for the address, street name, and then if it can't find anything show the block range where this address should exist. I've got the first two created and working, but I cannot get a locator to do block ranges. The fields line up properly and we have block ranges in our streets layer. What am I missing? I tried the dual range type locator to no avail.


texts: Cannot split vector across the 180degree line in QGIS

texts:
I am trying to split the vector I imported from geojson across the 180degree line and it gives me an error:
No features were split : If there are selected features, the split tool only applies to those. If you would like to split all features under the split line, clear the selection.
I noticed that I can split the polygon before the 180degree line but once I split it beyond it, I receive an error and cannot see a feature.
If I cut before the 180degree line (see the vertical split I did, which does not pass through 180 degrees), it splits the polygon as expected:

But, if I cut beyond the 180degree line, it returns an error:

I am using EPSG:3832 (PDC Mercator, centering the Pacific), QGIS version 3.16.1-Hannover
Has anyone encountered the same problem and know a solution?


texts: Overwrite existing raster layer using PyQGIS

texts:
QgsVectorFileWriter allows you to overwrite an existing layer with your output, but I have not found a similar option with QgsRasterFileWriter. Here's my code that is writing a raster layer out to a geopackage:
ds = ogr.Open(gpkg, True)
pipe = QgsRasterPipe()
provider=layer_in.dataProvider()

file_writer = QgsRasterFileWriter(gpkg)
file_writer.setOutputFormat('gpkg')
file_writer.setCreateOptions([f"RASTER_TABLE={layer_name}", 'APPEND_SUBDATASET=YES'])

if pipe.set(provider.clone()):
    projector = QgsRasterProjector()
    projector.setCrs(provider.crs(), provider.crs())
    if pipe.insert(2, projector) is True:
        if file_writer.writeRaster(pipe, provider.xSize(),provider.ySize(),provider.extent(),provider.crs()) == 0:
                feedback.pushInfo(f"ok")
        else:
                feedback.pushInfo(f"error: {file_writer.writeRaster(pipe, provider.xSize(),provider.ySize(),provider.extent(),provider.crs())}   {help(file_writer.setCreateOptions)}")
ds = None


texts: Issue using Shapefiles generated with QGIS in CadnaA

texts:
I frequently use the shapefiles from QGIS and ArcMap in other no GIS software, predominantly a noise software CadnaA.

CadnaA has the ability to import shapefiles, look for matching attribute names and associated the corresponding attributes to the CadnaA object type. Attributes such as 'Name', 'ID', 'LW_LI' etc... can be picked up and automatically inserted into the softwares object.

ArcMap does this perfectly. QGIS does not.

I am assuming that there is some sort of difference in the shapefile, but cannot determine what. I have seen elsewhere that it might be something to do with the shapefile type (3D perhaps?).

Edit:

Sorry this took so long. Here are examples of the shapefiles. The Arcmap ones import perfectly but the QGIS ones do not.

ArcMap Shapefiles

QGIS Shapefile

Anyone have any insight into this?


texts:
I have two types of image, one of them is a TRMM image with 0.25 degree cells and the other one is a MODIS image with 1 km spatial resolution.
How can I change my TRMM image resolution to 1 km without GCP points?


texts: Iterating over pyshp writer object records

texts:
The current branch of pyshp has removed the Editor class that combines reading and writing records. The justification is that removing the Editor has been planned for a while now, and might as well remove along with the other major new changes. Is there any way I can call a shapefile.Writer object and invoke  a records() type functionality that lists all present records in the writer object?
The only records introspection in pyshp currently is the .shpNum property that is just an incrementor:

def __dbfRecord(self, record):
    """Writes the dbf records."""
    f = self.__getFileObj(self._dbf)
    self.recNum += 1


texts: Map size very small in layout view in ArcMap 10.3

texts:
I am very new to GIS and ArcMap.
 I am currently working on small part of a DEM using ArcHydro toolkit. I have a map which is a clipped part from DEM. When I open it in data view, it looks fine. But for adding legends and scale when I use layout view, it is just visible at the corner of the page. Can't extend it for the print. What can I do?


texts: Interpolating z Values (height) from TIN (Raster) to line (vector) using QGIS

texts:
I have a raster (TIN) and I want to interpolate the heights (z-Values) to a line Shape. Everything that I find is interpolation between some points and making a raster. For example, Raster creation tools of SAGA or IDW or TIN interpolation, they all get a vector layer as an input file and make a Raster out of the values from attribute table.
Any idea?


texts: Offset a label inside a polygon with mapbox-gl style?

texts:
Im trying to offset labels on the border of a polygon to the inside of the polygon in a mapbox-gl style. My data has a series of polygons with shared borders, and the labels on the borders need to be offset in to indicate which polygon is which.

The carto-css equivalent of the effect im looking for is text-placement:line; text-dy:-5;.

I would expect that "symbol-placement": "line", "text-offset": [0,1] in the layout properties would offset the labels in to the polygon, but it offsets the labels down, not inside.

Adding "text-keep-upright": false offsets the labels in the correct direction, but it makes the text upside down on the bottom edge of the polygons.

Is there any way to offset the labels in, and keep them upright?

For the sake of completeness, here's a link to my style as it is now https://gist.github.com/JesseCrocker/17e11ccbe12c2b5393032f4585c13cdf


texts: In Google Earth Engine, can I export a layer with visual parameters added to a GeoTIFF?

texts:
I am trying to get a LandSAT-8 image of the certain region with cloud cover 2. I am using some visual parameters and adding my layer. I want to export this layer to GeoTIFF but at the moment I am exporting the image which is without the visualization parameters I have chosen. So is it possible to export the visual layer instead of the image? Following is my code:

var collection = ee.ImageCollection('LANDSAT/LC08/C01/T1_SR')
.filter(ee.Filter.eq('WRS_PATH', 25))
.filter(ee.Filter.eq('WRS_ROW', 31))
.filter(ee.Filter.lte('CLOUD_COVER', 2));
var image = ee.Image(collection.sort('CLOUD_COVER', false).first()).select(['B2']);
var visParams = {
  bands: ['B2'],
  min: 400,
  max: 1500,
  gamma: [1.5]
};
Map.addLayer(image, visParams, 'MyImage');
Export.image.toDrive({
    image: image,
    description: '25_31_B2',
    scale: 30,
    maxPixels: 1e13
});


texts: I have a polygon feature layer I am updating within a personal geodatabase. I have ordered my attributes in the table in a certain way for data entry. Every time I re-open the feature for editing, I have to re-order the fields to my preferred order. Could someone help me figure out how to make my preference the default?


texts: How to debug QGIS Server WMS "Download of capabilities failed"?

texts:
Currently I am using a QGIS server with a Linux CentOS system. I am attempting the create a WMS in QGIS and upload it to my website. I have followed the instructions in the "Quantum GIS User Guide: version 1.7.0 Wroclaw". I am trying to call the newly created WMS in QGIS by "Add Layer(s) from a server" button. I am very certain I am entering the correct address but I always get this error below

Could not understand the response. The wms provider said:
Download of capabilities failed: Error downloading

http://***************/cgi-bin/projects/qgismapserv.fcgi?SERVICE=WMS&REQUEST=GetCapabilities
Server replied Internal Server Error


texts: I have often used the Zonal Statistics plugin to calculate the mean pixel value of a raster. However, this time I have an additional Boolean (pixels are equal either to 1 or to 0) raster, which is overlapping with the other. I would like to calculate some kind of "conditional mean" of the former raster, i.e. the mean for those pixels which in the latter raster are equal to 0 and to 1. Can anyone suggest an approach? Is maybe R more appropriate that the the GRASS-GIS environment for doing so?


texts:
When I launch GeoNetwork, all changes in files of style (geonetwork.css, gn_search_less, gn_view.less,...) are lost, and the system generates a new file /static/gn_search_default.css.

How can I customize the GeoNetwork style (colors, style, ..) and have them persist every time I launch GeoNetwork?


texts: How to efficiently determine which of thousands of polygons intersect with a linestring

texts:
I have a shapefile that contains several thousand Polygons/MultiPolygons and another shapefile that has a variable number of LineStrings.  In the example below the blue Polygons are those that intersect with the red LineString and the brown Polygons are those that do not.

I am using the Shapely intersects binary predicate to check whether each LineString intersects with any one of the Polygons:

import fiona
from shapely.geometry import LineString, Polygon, MultiPolygon

# Open each layer
poly_layer = fiona.open('polygon_layer.shp', 'r')
line_layer = fiona.open('line_layer.shp', 'r')

# Create each Polygon/MultiPolygon
the_polygons = []
for poly in poly_layer:
    if poly['geometry']['type'] == 'Polygon':
        the_polygons.append((poly['properties']['poly_id'],
                             Polygon(poly['geometry']['coordinates'][0])))
    elif poly['geometry']['type'] == 'MultiPolygon':
        all_polygons = []
        for poly_coords in poly['geometry']['coordinates']:
            all_polygons.append(Polygon(poly_coords[0]))

        the_polygons.append((poly['properties']['poly_id'],
                             MultiPolygon(all_polygons)))
    else:
        print poly['geometry']['type']
        raw_input()

# Create each LineString
the_lines = []
for line in line_layer:
    the_lines.append(LineString(line['geometry']['coordinates']))

# Check for Polygons/MultiPolygons that the LineString intersects with
covered_polygons = {}
for poly_id, poly in the_polygons:
    for line in the_lines:
        if poly.intersects(line):
            covered_polygons[poly_id] = covered_polygons.get(poly_id, 0) + 1

The above code takes a long time to run, and being new to Fiona and Shapely I am not sure that this is the most efficient way to accomplish things.  So is there a better way to do what I want?

Edit: The key to speeding up the computations was to use a spatial index.  The linked question does not have what I believe is the fastest solution - using GeoPandas.


texts: How to filter Legend to show only visible features?

texts:
Is it possible to get the legend in a pagelayout to only show patches/swatches for features fall within the visible extent of the Map?

The solution needs to work on FeatureLayers and ArcGIS Server layer. It's possible to set a property on the AutoVisibility to filter the Legend to only show the layers that are visible, but if the layer is visible and only a few of multiple symbols in the renderers is visible the legend gets too many items to be useful.

Update: This seems to be a question that has been asked several times before on the ESRI-forum, found one solution that seems to works on a FeatureLayer if it has a IUniqueValueRenderer Compressing a legend to show displayed features only but we are using both FeatureLayer and ArcGIS Server layers, so that only solves part of our problem.

We are trying to find a way to hide
and show the IElement from ILegendItem.Graphics
, but it feel like a deadend.
On the server ILegendItem.Layer it
might be possible to cast to a
MapServerIdentifySublayer and use
Identify to somehow find the visible
objects.

If that's possible its only a matter of connecting the two.


texts: How to find direction points using Google Maps API 3?

texts:
I am drawing a route in Google Maps. I have passed the start point, way points and end points. My route is drawing correctly. I am using Google Maps API 3.

request = {
  origin: originAddress,
  destination: destinationAddress,
  waypoints: waypoints,
  optimizeWaypoints: false,
  travelMode: google.maps.DirectionsTravelMode.DRIVING
};

directionsService.route(request, function(response, status) {
  if (status == google.maps.DirectionsStatus.OK) {
    directionsDisplay.setDirections(response);
  }
}

If I change the route direction, then how to find the changed direction points' latitude and longitude?

Any help will be appreciated.


texts:
I have download latest version (0.46) of osmosis from
https://wiki.openstreetmap.org/wiki/Osmosis#Downloading
and mapsforge-map-writer-0.9.1.jar from http://search.maven.org/#search%7Cga%7C1%7Cmapsforge-map
and map ("AM.osm") from
http://be.gis-lab.info/data/osm_dump/dump/latest/

I'm on Windows-7 64-bit Professional and JDK 1.8u45
I unpacked osmosis to: c:\Progs\osmosis
I created a folder called plugins in c:\Progs\osmosis (>> so: c:\Progs\osmosis\plugins)
In the plugins folder I put the file mapsforge-map-writer-0.9.1.jar
I put the file AM.osm.pbf (and AM.osm) into c:\Progs\osmosis

I then opened a command-prompt window and executed the following commands:

PS C:\Progs\osmosis> osmosis --rb file=AM.osm.pbf --mw file=AM.map

PS C:\Progs\osmosis> osmosis --rx file=AM.osm --mw file=AM.map type=hd

I got this error:

ьр  23, 2018 9:46:43 AM org.openstreetmap.osmosis.core.Osmosis run
INFO: Osmosis Version 0.46
ьр  23, 2018 9:46:43 AM org.openstreetmap.osmosis.core.Osmosis run
INFO: Preparing pipeline.
ьр  23, 2018 9:46:43 AM org.openstreetmap.osmosis.core.Osmosis main
SEVERE: Execution aborted.
java.lang.NoClassDefFoundError: gnu/trove/procedure/TShortIntProcedure
        at org.mapsforge.map.writer.model.MapWriterConfiguration.loadTagMappingFile(MapWriterConfiguration.java:351)
        at org.mapsforge.map.writer.osmosis.MapFileWriterFactory.createTaskManagerImpl(MapFileWriterFactory.java:58)
        at org.openstreetmap.osmosis.core.pipeline.common.TaskManagerFactory.createTaskManager(TaskManagerFactory.java:60)
        at org.openstreetmap.osmosis.core.pipeline.common.Pipeline.buildTasks(Pipeline.java:51)
        at org.openstreetmap.osmosis.core.pipeline.common.Pipeline.prepare(Pipeline.java:112)
        at org.openstreetmap.osmosis.core.Osmosis.run(Osmosis.java:86)
        at org.openstreetmap.osmosis.core.Osmosis.main(Osmosis.java:37)
        at sun.reflect.NativeMethodAccessorImpl.invoke0(Native Method)
        at sun.reflect.NativeMethodAccessorImpl.invoke(Unknown Source)
        at sun.reflect.DelegatingMethodAccessorImpl.invoke(Unknown Source)
        at java.lang.reflect.Method.invoke(Unknown Source)
        at org.codehaus.plexus.classworlds.launcher.Launcher.launchStandard(Launcher.java:330)
        at org.codehaus.plexus.classworlds.launcher.Launcher.launch(Launcher.java:238)
        at org.codehaus.plexus.classworlds.launcher.Launcher.mainWithExitCode(Launcher.java:415)
        at org.codehaus.plexus.classworlds.launcher.Launcher.main(Launcher.java:356)
        at org.codehaus.classworlds.Launcher.main(Launcher.java:47)
Caused by: java.lang.ClassNotFoundException: gnu.trove.procedure.TShortIntProcedure
        at org.java.plugin.standard.StandardPluginClassLoader.loadClass(StandardPluginClassLoader.java:330)
        at java.lang.ClassLoader.loadClass(Unknown Source)
        ... 16 more


texts: Python select by attribute using a variable counter

texts:
I am brand new to python and trying to learn how to do some simple loops but I have hit a snag.  Below is what I am working with, I am trying to select the polygon row by row, but I keep getting EOL errors and various other syntax issues as I try to troubleshoot.  If I remove """ findMe """ and replace with """ "pgNum" = 1 """ it works fine.  How can I change 1 to be x?

I am using the python window within ArcGIS

Thanks.

import arcpy
workSpc = (r'C:\Users\adams\Downloads\tmp')
smooth = 'MNR_89_Smooth'
grid = 'test_grid'
x = 1

if x == 1:
    findMe = '"pgNum" = ' + str(x)
    print findMe
    arcpy.SelectLayerByAttribute_management("test_grid","NEW_SELECTION",""" findMe """)
    x = x + 1
    print x


texts:
I have a project with just an OSM topo map and want to export it via print composer with a scale and a frame. Very basic stuff. However, at some point when trying to get some ticks for lon and lat at the borders of the map, QGIS didn't respond anymore. When opening the project, it always crashes or doesn't respond.
That's not the first time I encounter this problem and not the only version of QGIS. Atm the issue occurrs in QGIS 2.18.14.


texts: How to hide system tables of oracle in sde connection?

texts: I connected to oracle 11g r2 with ArcGIS 10.2 using SDE connection. The connection is ok but in the connection i see system tables of  APEX_030200 and other  schemas. Same question is in the Geonet but nobody answer the question in the forum. How can i hide system tables in sde user ?


texts: R: how to export my Kriging result to raster?

texts:
I want to export results of my Kriging calculation to raster or some other format. My aim is to take this converted file and intersect it with another SpatialPointDataFrame (location of the point doesn't have to be the same, so that is the reason of conversion). The goal is equivalent of "Extract values To Points" in ArcGIS - thus, merge information from kriging and another SPDF for one locations.

For second part, I imagine just to use raster::intersect(krig, SPDF).
But how to get raster data for my Kriging? I tried to raster::rasterize

r <- raster(ncols = length(ux), nrows = length(uy))  # length(ux), length(uy)
krig.rst<-rasterize(krig.est, r, "var1.pred")

and

krig.rst<-rasterize(krig.est, r, "var1.pred", fun = function(x, ...) c(min(x))

The raster seems to be calculated, but I can't display it.

Please, is the problem in my raster dimension (nrow, ncol, ncell)? How can I assigned it? My first thought was that 1 cell of matrix (from Kriging) will correspond to 1 pixel of raster with assigned value, but I am obviously wrong.

Long reproductible exemple adapted from here: https://stackoverflow.com/questions/37177017/set-automatic-colors-levels-in-filled-contour-kriging

library(sp)
library(gstat)
library(RColorBrewer)

# https://rpubs.com/nabilabd/118172
# packages for manipulation & visualization
suppressPackageStartupMessages({
  library(dplyr) # for "glimpse"
  library(ggplot2)
  library(scales) # for "comma"
  library(magrittr)
})

data(meuse)

# create spdf
meuse.spdf<-meuse

# convert to spdf
coordinates(meuse.spdf) <- ~ x + y

# calculate voronoi tesselation - will be needed to create underlying point data
voronoipolygons = function(layer) {
  require(deldir)
  crds = layer@coords
  z = deldir(crds[,1], crds[,2])
  w = tile.list(z)
  polys = vector(mode='list', length=length(w))
  require(sp)
  for (i in seq(along=polys)) {
    pcrds = cbind(w[[i]]$x, w[[i]]$y)
    pcrds = rbind(pcrds, pcrds[1,])
    polys[[i]] = Polygons(list(Polygon(pcrds)), ID=as.character(i))
  }
  SP = SpatialPolygons(polys)
  voronoi = SpatialPolygonsDataFrame(SP, data=data.frame(dummy = seq(length(SP)), row.names=sapply(slot(SP, 'polygons'),
                                                                                function(x) slot(x, 'ID'))))
}
meuse.voro <- voronoipolygons(meuse.spdf)

# create underlying grid
s.grid <- spsample(meuse.voro, type = "regular", n = 6000)

# calculate kriging
 # create variogram
lzn.vgm <- variogram(log(zinc)~1, meuse.spdf) # calculates sample variogram values
lzn.fit <- fit.variogram(lzn.vgm, model=vgm(1, "Sph", 900, 1)) # fit model

plot(lzn.vgm, lzn.fit) # plot the sample values, along with the fit model

# calculate kriging
lzn.kriged <- krige(log(zinc) ~ 1, meuse.spdf, s.grid, model=lzn.fit)

# extract the unique x and y locations in the grid
ux<-unique(coordinates(lzn.kriged)[,1])
uy<-unique(coordinates(lzn.kriged)[,2])

# extract the predicted values and format var1.pred into a matrix of gridded values
predmat3 <- matrix(lzn.kriged$var1.pred, length(ux), length(uy))

# display the data

filled.contour(x = seq(0,1, length.out = nrow(predmat3)),
               y = seq(0,1, length.out = ncol(predmat3)),
               z = predmat3,
               col = brewer.pal(6,"Purples"), nlevels = 5)


texts:
It might be a stupid question but I need to ask it. On the latest version of OpenLayers which is using OSM tiles, when you are quickly zooming map the tiles which are not loaded yet gets blurry:

On very, very old OpenLayers 2 it was looking like this:

I like the new behaviour a lot but one of the customers is asking me if there is a way to make it as it was before? He doesn't like those blurred tiles but prefers white squares. I also remember that in the past there was some kind of "x" icon in one of the tiles squares when they have not been loaded.

Bing tiles load fast and there is no problems but OSM is working slowly those days.


texts:
I want to design Editor Forms in QGIS 3.4 programmatically via QGIS' Python API.

Step 1: Setting Editors and human readable names for fields works fine via calls to

QgsVectorLayer.setEditorWidgetSetup()
QgsVectorLayer.setFieldAlias()

Step 2: Introduce Tab boxes in my Form works as well
I tried to switch to using QgsEditFormConfig.EditorLayout.TabLayout, but this does not work:

editFormConfig = layer.editFormConfig()
editFormConfig.setLayout(1)
editFormConfig.addTab(....)
layer.setEditFormConfig(editFormConfig)

My question is: how can I connect editorWidgetSetup settings from step 1 with those Tab boxes, setting the order of elements and which item is on which tab box?


texts:
Querying through the EPSG database, I can see the tables epsg_coordoperation where transforms and conversions between different CRS are stored, and epsg_coordoperationmethod which stores the mathematical methods used by the transformations. Correct me whether I'm wrong.

To make an example, looking through the details of the conversion form ED50 to WGS84 over Sicily (http://www.epsg-registry.org/, --> retrieve by code --> 1143), I can see that the transformation involves simple shifts on the 3 cartesian axis. Operation is ED50 to WGS84, operation method is Geocentric Translation.

My statements which I would like to be either confirmed or corrected are:

coordinate operations are always meant to reach a target CRS;

coordinate transformations involve a change of datum, e.g. [geocentric<-->geocentric] or [projected<-->projected], etc.
..as coordinate conversions do, but regarding a [projected<-->geographic] error-free translation;

one could apply arbitrary operation methods to a georeferenced dataset, say a simple arbitrary shift on one or more dimensions in the most simple case, given that the requirementes are met.

My questions:

Considering only the operation methods, does applying method X to CRS A always take to a different CRS? But which one? E.g. a Lambert Conic Conformal (2SP) - EPSG:9802 - projects ellipsoidal coordinated to a projected XY map, but to which CRS? Or again, a shift on an ED50-based image may take to a WGS84-based image, but only if the shifts are the ones defined in the correct operation (--> Parameter values).
Bob wants to reproject a coverage from a source to a target CRS. Are there cases in which he would prefer to specify its own path of CRS transform concatenation instead of the default method offered by foo library?

Thanks in advance for any effort.


texts:
I am trying to work out the relationship between the 'processing cell size' and the result from Tabulate Areas. If one uses the same as the raster cell size, and this is the same for both rasters, the area output is correct.

However if one uses another value, for example use a processing size of 50 to a cell size of 25x25. The result is roughly four times higher. This suggests it is counting the overlapping cells and calculating the area by multiplying it with the 'processing cell size' value, over riding the rasters own cell size.

Can this be right? I would have expected the processing cell to size to be the unit of measurement, e.g. to allow adjustment of the resolution of the area calc - if for example one uses different resolution rasters or a vector set.

Any clarification would be welcome.


texts: Create toggle buttons to activate individual query and identify tasks

texts:
Can you please help with toggling query and identify tasks?
For example,
Parcel: button
Footpring: button
CIP: button

whenever one of the buttons is checked, all others are  inactive.

D.E.Wright provided a test link that is no longer active here: Toggle Toolbar Button and Activate Identify

Thanks for any help!

EDIT: I am attempting to build a javascript api application that involves several functionalities: draw tools, query, find, etc. Now, whenever I click on the map it evokes the identify task, whether I want to or not. I would like to control when this task in called.

I am trying to figure out the best way to display all of my code for your viewing. I have an html file and a series of .js files along with the css.


texts: How to style custom elements to match the ArcGIS widgets?

texts:
Using the ArcGIS Server v4 JavaScript API, the instructions for creating a custom widget are very complex.

In my case I simply need to add some custom buttons, and style them to match the existing widgets, so these instructions are overkill.

Is there a simple way to style buttons to match the out-of-the-box widgets?


texts: Gdal_translate command for SDS having spaces in its name

texts:
Can anyone help me with using gdal_translate on an HDF4 dataset that a space?

One SDS in the the HDF is:

HDF4_EOS:EOS_GRID:"MOD13Q1.A2004001.h08v05.005.2007234070731.hdf":MODIS_Grid_16DAY_250m_500m_VI:250m 16 days NDVI

I want to use gdal_translate to extract this as a tif or envi file and I use the following command:

HDF4_EOS:EOS_GRID:"MOD13Q1.A2004001.h08v05.005.2007234070731.hdf":MODIS_Grid_16DAY_250m_500m_VI:250m 16 days NDVI ndvi.envi

This displays a message "too many command options" and displays the general gdal_translate template and "The following format drivers are configured and support output:" list.

Where am I going wrong? I have used the same command another dataset that does not have spaces like that of between "16 days NDVI" and that works perfectly.


texts: Add Specific Map Layout Name to Document Path in ArcGIS Pro

texts: In ArcGIS Pro (using 2.6.2), the <dyn type="project" property="path"/> dynamic text only provides the path of the entire project. It does not include the actual map or layout within that project. I would like to include this so when I export to PDF, another user can know exactly which map the PDF is sourced from. Is there a way to append this information to the file path, using dynamic text or another method?


texts: Creating multiple input features in ModelBuilder?

texts:
I am a novice GIS user.

I am trying to clip multiple input features using the ModelBuilder and expecting multiple outputs.

How do I create multiple input (I've seen the several overlapped oval as shown below)?


texts:
I am new to GIS, in school for it actually, but I come from a photography / graphic design background.  So that said, I always thought that vector graphics were not lossy (scale up or down, doesn't matter) and raster graphics were lossy (or have the potential for being lossy if attempted to scale up).  It's certainly true in the graphic design world, do the same rules apply to GIS?
Earlier when I was explaining to someone on another website that one of the differences between "raster vs vector" was that they were "lossy vs not-lossy" respectively, another person felt the need to jump in with this:

"You are not correct.
First, raster is usually performed at a fixed
resolution and is only lossy with some compression techniques. It's
the compression technique, not the raster, which is lossy. Because of
compression and difficulty of displaying a variable dataset in vector
without losing information, some rasters are more efficient than
vectors in displaying info.
Second, vectors can be lossy based on line and area simplification techniques
where unnecessary points are removed because they aren't relevant for
display scale. It can also be lossy in the precision of a coordinate"

Not really sure what some of that language even means.
Can someone here clarify which one is lossy and which is not lossy?


texts: GeoServer not wrapping mutliline text on \n

texts: I am trying to display labels from attribute which contents multiline texts. The lines in texts are divided with \n. Can GeoServer wrap these texts on \n?


texts: Calculate river distance to headwater for every pixel?

texts:
I already posted a similar question a few weeks before. But since my dataset changed and I still havent solved it, I thought its good to post a new question instead of editing the old.
Using ArcGIS 10.4

I am now working with a DEM. With the ArcHydro Toolbox I calculated subcatchments which I need for running a habitat model later. Now I need the maximum distance of every subcatchment to the source of the river (headwater).

I was thinking a tool like "path distance" or "flow length" should do the trick. And then maybe zonal statistics to calculate it for my subcatchments. But I cant get it right.

Any ideas?

UPDATE:
With Lovette's tip with the weightraster I finally got it done. But since my subcatchment are not 100% correct (I guess that is normal) and my pixels are a little big (I need 50x50), the results of the zonal statistics are then a little weird sometimes. In picture 1 you can see that the blue (larger distance) pixel pf the main stream also touches the subcatchment of the side stream. Because the zonal statistics always takes the maximum, the whole catchment of the small side river gets a distance value that is way too high. The side streams should always have shorter distances.

Any ideas what to do about that? Can I maybe give the subcatchments values manually. This would take me a while, but from the numbers it would be doable.


texts: Displaying KML attributes (Images) on Google My Map

texts: I don't know if i am the right place to ask this question, i want to display images attached with each attribute in the kmz file on GoogleMyMaps. I can view the images on google earth but when i upload the same kmz file on GoogleMyMap other attributes are displayed but it is not displaying the image.


texts: What is a Format of csv plugin datasource string in mapnik

texts:
I have tried that Format Is that is a correct format to pass to CSV plugin,the sample C++ code is given Below:

std::string csv_string("wkt,Name POLYGON((5.123456 51.6, 6.8 52.6, 7.8 51.6, 5.123456 51.6)),Roads" );

  {     parameters p;
        p["type"]="csv";
        p["inline"]=csv_string;
        layer lyr("Roads");
        lyr.set_datasource(datasource_cache::instance().create(p));
   }


texts: Usage of interpolation methods with DEMs: Linear vs. natural neighbors

texts:
Using the ArcGIS tool "TIN to Raster" (3D Analyst Tools --> Conversion  --> From TIN) to work on my Digital Elevation Model I was wondering about the interpolation methods.

I can choose either or:

linear interpolation ("Default. Calculates cell values by applying linear interpolation to the TIN triangles")
natural neighbors interpolation ("Calculates cell values by using natural neighbors interpolation of TIN triangles...The algorithm used by the Natural Neighbor interpolation tool finds the closest subset of input samples to a query point and applies weights to them based on proportionate areas to interpolate a value")

[From ESRI Dokumentation]

So far, I understand the technical difference between both and how the algorithms work. But I am struggling how to interpret that knowledge. If I choose one of the interpolation methods, what  consequences does the choice of one method mean to my result? Are there use cases, when to take either the one or the other method? What are advantages or disadvantages of each method?


texts: Google earth: adding multiple images to a placemark

texts:
I have a couple of images that I want to add to a placemark. The problem is when I add the multiple images, google earth displays them one below the other.
Is it possible that the image are displayed side by side whenever I click the placemark icon.

And If there a dozen of images for the same placemark can I create a "slideshow" of the imagery that displays only one image at a time and swaps when user clicks the placemark button again?

Any help is appreciated. Thank you.


texts: Exporting XY data to .shp and I lose my points

texts: I have coordinates (sampling points) as part of a larger excel file that I am trying to convert to a .shp file to overlay on a map.  When I add them to Arc as XY data they display just fine, but when I export to .shp file I lose my points.  When I select zoom to layer it zooms way way way out and still displays nothing.  Any ideas about what might be going wrong?  I am using ArcMap 10.1


texts: Distance from a point to a poly line and attributes from the poly line GIS

texts:
I'm trying to figure out what tool would be best to get the distance from multi-points from multi polylines and see which attributes are closest to what point.

I'm doing a project where I have historical tornado data and would like to know what intensity and how close these historic tornados were close to rural schools.


texts: How to locate file path to original layers from feature layers in AGOL web map?

texts: I was given an AGOL web map that was created by someone else within my organization that is no longer here, and I have no idea where the original editable layers in the web map exist on our network that were published to this web map. Is there a way in AGOL to see the original file paths for published feature layers? I want to be able to bring in the live data from this web map into future mxds and not have to just export the data as shapefiles from the web map.


texts: How to Convert Just a Sulawesi Map to Leaflet

texts:
I'm still develop a Web GIS that show population distribution in form of heat map. I'm use Quantum GIS with Open Street Map & Q2Leaflet Plugin as a tool to accomplish it.
The problem is, I just want to show the population distribution of Sulawesi, but when I download data via Open Street Map, I ended up download the World Map. I just want to use Sulawesi Map, convert it to leaflet & embed it to my Web GIS.

How should I do to convert ONLY Sulawesi Map into leaflet?


texts: MXD Layer Management When Publishing Services with ArcGIS

texts:
Currently have an MXD with 130 layers in the table of contents. All of these layers need to exist in the service when it gets published each time an update is performed. I'm in the understanding that when a service is published to ArcGIS for Server/Enterprise, that it does a complete replacement of the service with the new service. Is there a way to just add to the existing service?

What I'm trying to accomplish is to remove all 130 layers in the MXD that I don't need since they already exist in the service. I just want to push the new layers up and add them to the current service. This is a rolling update process, so eventually, I will have 500+ layers, and I don't want to manage all that in the MXD each time I update.

How can I do this? What do I need to do with the MXD and/or GDB to accomplish this?

EDIT
The 130 layers (stored in a GDB as FCs) exist in an MXD and represent county level health data. Every layer has 3140 or so records with various fields. No application of symbology at the MXD level as that is being done at the service endpoint level. The service exists to feed an application tool where users search for data, select it and display a thematic map (only one at a time). There will be no intention to show two or more at the same time.


texts: Does contrast stretching help image segmentation and classification?

texts:
I am trying to identify and classify marine litter from drone RGB images at beaches. Due to the small size of litter and only having RGB images available, my supervisor recommended to carry out contrast enhancement. Does this actually help to get better image segmentation and classification?
I have read that it only changes the display view, so I am not sure how this could somehow help.


texts: Defining buffer size when labeling text in QGIS?

texts:
I need to set label buffer size to near zero but I dont want to overlap labels. How to define buffer size when labeling text?

Currently I am using labeling settings which is allowing all available text labeling showing in map. Text labeling are now fine for me, because there is both text classes, address numbers and street names. Unfortunately there is now also overlap labels.

When I remove "Show all labels for this layer" -settings, overcrossing labels disappeared but also street name text although it should not touch address numbers.

How to set label buffering settings better if I want text classes very near together but no overlapping?


texts: Loading Tiger Data, and populating missing table data (such as tabblock)

texts:
I'm working with postgis and followed many tutorials out there to try to use tiger data.

this one seems to be the best:

https://bluehatrecord.wordpress.com/2015/07/06/installing-the-tiger-geocoder-with-rhel-6-6-postgresql-9-4-postgis-2-1-7-and-the-tiger-2013-census-data/

To save time and make sure things are working, I reduced it to a smaller set of states to populate my Tiger data.

Looking in the output, I'm seeing things like this:

Unable to open tl_2016_44007_addr.shp or tl_2016_44007_addr.SHP.
tl_2016_44007_addr.dbf: shape (.shp) or index files (.shx) can not be opened, will just import attribute data.

and

Unable to open tl_2016_36019_featnames.shp or tl_2016_36019_featnames.SHP.
tl_2016_36019_featnames.dbf: shape (.shp) or index files (.shx) can not be opened, will just import attribute data.

and sometimes:

HTTP request sent, awaiting response... 304 Not Modified
File ‘www2.census.gov/geo/tiger/TIGER2016/PLACE/tl_2016_11_place.zip’ not modified on server. Omitting download.

I assume the second error means that it already downloaded the file, so it skipped it?

I ran the script again with just DC as a state, and even went with the author's example for a sanity check:

SELECT g.rating, ST_X(g.geomout) As lon, ST_Y(g.geomout) As lat, (addy).address As stno, (addy).streetname As street, (addy).streettypeabbrev As styp, (addy).location As city, (addy).stateabbrev As st, (addy).zip FROM geocode('1600 Pennsylvania Avenue Northwest, Washington, DC 20500') As g;

Success! I got back the matching data, which is good, but when I ran his second example to confirm:

SELECT g.*, ST_X(g.geomout) As lon,ST_Y(g.geomout) As lat,(addy).*, bg.tabblock_id AS FIPS FROM geocode('1600 Pennsylvania Avenue Northwest, Washington, DC 20500') As g
JOIN tabblock As bg
ON (g.geomout && bg.the_geom)
WHERE ST_Within (g.geomout, bg.the_geom);

it returned back 0 records. I did a select * from tabblock and it returned 0 records.

What table is that, and should it normally download based on the script?

If not, how is this table populated?


texts: Preserving attributes of KML file when loaded into QGIS?

texts:
How do I upload a KML file into QGIS with the attribute data intact? It is present in Google Earth but not when imported into QGIS?

http://data.gov.uk/dataset/sites-of-biological-importance-sbi-in-greater-manchester-as-at-2015-04-28

This is the website which contains the data I need. I need to show the site names on my map, but I am currently unable to do this.
My coordinate reference system is: OSGB 1936/BNG (OTFR Enabled).


texts: QGIS: add a widget to a toolbar other than the Plugins toolbar

texts:
I have a QToolButton that I want to add to an existing toolbar. Using the Python Console in QGIS 2.14.12 to debug, here's the code fragment that I'm having a problems with:

iface.toolButton=QToolButton()
iface.toolButton.setIcon(QIcon("C:/Users/Rudy/.qgis2/python/plugins/sandbox/icons/action1.png"))
iface.addToolBarWidget(iface.toolButton)

The toolButton displays in the Plugins toolbar, just like the documentation for addToolBarWidget says :-). However, I want it to display in another existing toolbar, not the Plugins toolbar. There is the addUserInputWidget call, but all that happens when I use that call is that a user input box flashes momentarily on the screen, then disappears. The result is that the toolButton doesn't show up anywhere. Isn't there an API call that would place widgets in any toolbar that I specify?


texts: How to convert a column of numeric IDs to letters?

texts:
I'm working with qgis 1.8.0. I divided a map into 500 features (small rectangles), and now I have the doubt if there is any way to change the ID of the attributes from numbers to letters? When I did the division of the original map, each of the 500 divisions got one numeric ID (please refer to column 3 of image) and now I'm trying to load that map to STATPLANET PLUS and when working with that program, the ID may not be numerical. Any ideas?

Regards!!
actual attribute table


texts: How to read the app.config settings in ESRI ArcObjects ArcMap Add-In

texts:
I am developing custom toolbar and buttons for ArcMap 10.1 utilizing VS2010 and ArcObjects SDE 10.1 extensions.
I put some application settings in the app.config file and I want to make them available in the runtime.
When I execute

private static readonly string appSetting1 = System.Configuration.ConfigurationManager.AppSettings["appSetting1 "];

the value is null.

When I execute

AppDomain.CurrentDomain.SetupInformation.ConfigurationFile

the value is

  "C:\Program Files (x86)\ArcGIS\Desktop10.1\bin\ArcMap.exe.Config"

Should I put my application settings is the ArcMap.exe.Config in order make them available to the ArcMap 10.1 custom toolbar or there is another way to do this?
Thanks


texts: CartoDB: cluster circles aren't proportional

texts:
I'm trying to visualize thousands of records on a map on CartoDB using the cluster feature but circles with 100 or 1000+ records have the same dimensions.
I'd like to have more discrimination in circle dimensions, is there a way to set limits for cluster classification? Or, even better, draw them exactly proportionally?

CSS code:

/** cluster visualization */

#fumo_globale{
  marker-width: 12;
  marker-fill: #FD8D3C;
  marker-line-width: 1.5;
  marker-fill-opacity: 1;
  marker-line-opacity: 1;
  marker-line-color: #fff;
  marker-allow-overlap: true;

  [src = 'bucketE'] {
    marker-line-width: 1;
    marker-width: 24;
  }

  [src = 'bucketD'] {
    marker-line-width: 1;
    marker-width: 34;
  }

  [src = 'bucketC'] {
    marker-line-width: 1;
    marker-width: 44;
  }

  [src = 'bucketB'] {
    marker-line-width: 1;
    marker-width: 54;
  }

  [src = 'bucketA'] {
    marker-line-width: 1;
    marker-width: 64;
  }
}

#fumo_globale::labels {
  text-size: 0;
  text-fill: #fff;
  text-opacity: 0.8;
  text-name: [points_count];
  text-face-name: 'DejaVu Sans Book';
  text-halo-fill: #FFF;
  text-halo-radius: 0;

  [src = 'bucketE'] {
    text-size: 12;
    text-halo-radius: 0.5;
  }

  [src = 'bucketD'] {
    text-size: 17;
    text-halo-radius: 0.5;
  }

  [src = 'bucketC'] {
    text-size: 22;
    text-halo-radius: 0.5;
  }

  [src = 'bucketB'] {
    text-size: 27;
    text-halo-radius: 0.5;
  }

  [src = 'bucketA'] {
    text-size: 32;
    text-halo-radius: 0.5;
  }

  text-allow-overlap: true;

  [zoom>11]{ text-size: 16; }
  [points_count = 1]{ text-size: 0; }
}


texts: Update Point using WFS Transaction with GeoServer

texts:
I am trying to update a Point inside my PostGIS database that is connected with GeoServer.

I specify the new coordinates in the WFS-Transaction post message like this:

<gml:Point srsDimension="2" srsName="urn:x-ogc:def:crs:EPSG:4326">
<gml:coordinates decimal="." cs="," ts=" ">0,0</gml:coordinates>
</gml:Point>

This works just fine.

However I am confused about the next step. The Points inside my PostGIS database are not in EPSG:4326 but EPSG:28992. So I tried changing the srsName to EPSG:28992 but this just makes GeoServer throw an exception.

How am I supposed to update Points using the EPSG:28992 projection?

If I update the points using srsName="urn:x-ogc:def:crs:EPSG:4326" GeoServer throws an PointOutsideEnvelopeException which is understandable.

Am I missing something here? Can I only update points using EPSG:4326?


texts:
In QGIS 2.13.2 I have a google map and a shapefile. Now I want to print it with all the necessary items, but I have a big problem with the coordinates.
First of all, the coordinates just appear when I use the WGS84 / Pseudo Mercator (EPSG: 3857). If I try to use the normal WGS 84 nothing appears. Why is that?

The other problem is, that my coordinates show something completely different than they should be. I looked them up and its somewhere in Africa, although I need the coordinates for a place in Germany. Did I make a mistake with georeferencing the shapefile?


texts: How to change 0,360 WKT string to -180,180?

texts: I'm reading POLYGON WKT strings using OpenLayers.Format.WKT().read(). The polygons use 0,360 longitudes but my map is -180,180. When I display the polygons on the map they are displayed outside of the map. How can I display them on their correct location? I tried OpenLayers transforms and Proj4js transforms (plus  defining the 0,360 CRS using +lon_wrap=180) but nothing seems to work! Help would be appreciated ;-)


texts: How to specify data defined rounded background corners for labels in QGIS 3?

texts:
In QGIS 2 I could specify data defined rounded corners for labels background, using a string, for example '5, 5' :

In QGIS 3, this same syntax is correctly recognised by the expression generator, but labels corners remain square :

Does anyone know the right syntax or is this a bug ?


texts: GRASS "New Mapset" and "Open Mapset" etc is greyed out

texts:
All tutorials I have looked at say to just click "New Mapset" to begin a new one but obviously I can't do this to even begin.

I'm not sure if there's something you need to do before creating a new Mapset: I've tried it with a (OS OpenData) map already in QGIS and with a blank page and nothing seems to allow me to click on it.

What am I doing wrong?


texts: Using Join Field in ArcPy to merge multiple tables

texts:
I'm trying to merge multiple tables using the Join Field function in ArcPy. This script currently runs perfectly to merge the multiple tables but gives the columns new names. However, I will like the file name appended to the new columns as they are added to the table.
How can I modify my code to add columns with the file name appended?
import arcpy
arcpy.env.workspace = "d:/INDICES5/Pixelnum/SAVI"
tables = arcpy.ListTables()
in_data = "d:/INDICES5/Pixelnum/SAVI/PixelnumSAVI2020.dbf"
for table in tables:
     arcpy.JoinField_management(in_data,"Pixel_num", table,"Pixel_num")


texts: ArcGIS integration with OpenText eDocs DM?

texts:
I work for a municipality and all of our internal documents are stored in a Document Management system, OpenText eDocs DM.
I would like that in my internal Mapping Application (ArcGIS Server / Silverlight API), I could link my Parcels with its legal documentation stored in eDocs.

Has anyone been able to implement such a thing?


texts: Using FeatureCollection to mask an image in Google Earth Engine?

texts:
I want to classify crop types using Sentinel-2 image in GEE. Before starting the classification, however, I want to update the mask of the image, so that water, urban areas, forest and other non-arable land is masked. For that purpose I uploaded a shape file with the arable land and imported it as a FeatureCollection.

My question is is it possible to use this FeatureCollection to update the mask of the image dataset. I know of the updateMask() function but it takes an image as its argument.


texts:
Urban rate (the % of population living in urban areas) is a rather common demographic indicator, at nation-level resolution.
I'm trying to locate sub-national urban rate data (NUTS level 2 to be precise)

The raw data is accessible more or less in country census data, yet appears in a plethora of languages and sources.

So if anyone knows of a source for urban rates, I would appreciate it.


texts: Align OpenLayers raster with OS Mapping Vector

texts:
I have an OS vector map of a town in Ireland. When I load the Google Satellite raster from the OpenLayers plugin, the results are slightly off, by about 2m, probably due to aerial perspective.

I'm trying to digitise information from the Satellite raster onto the OS Mapping and it would be easier if I could force the Satellite raster to move temporarily so they align in that perspective.

Is this possible?


texts: Checking status of layer calculation using API

texts:
I wrote a python script that updates the values in a dataset and then returns a static image of the map it's connected to.

The problem is if I download the image too soon after updating the dataset, one of the visualization layers on the map that requires a lot of calculations won't show up.

I don't want to have to set a manual wait time in the code because the size of the data I'm uploading could range from 100 to 10,000 data points.

Is there an API url I could call that would know when a layer is being updated and when it has finished?

I have tried using vizualization url but could not find anything in the JSON data about an update status.

img_url1 = 'https://%s.carto.com/api/v2/viz/%s/viz.json/?api_key=%s&' %(self.options['u'], map_id, self.options['k'])
r1 = requests.get(img_url1)
print(json.dumps(r4.json(), indent=2))

On a similar thread someone sugested using this to refresh the the_geom data but I had no luck there either:

refresh_url = "https://{username}.carto.com/api/v2/sql?q=UPDATE table SET the_geom = cdb_latlng(latitude, longitude)&api_key={api_key}".format(username=self.options['u'], api_key=self.options['k'])
r4 = requests.get(refresh_url)
print(json.dumps(r4.json(), indent=2))

I have tried using the named map API but the 'Last Updated' field shows the time when the update started not finished, and the 'Date' field just shows the current time.

img_url2 = 'https://{user}.carto.com/api/v1/map/named/{template}/jsonp?auth_token={auth_token}&callback=callback'.format(user=self.options['u'], template=tpl_name, auth_token=self.options['k'])
r2 = requests.get(img_url2)
print(r2.headers)

'Date': 'Mon, 10 Feb 2020 20:27:14 GMT',
'Content-Type': 'text/javascript; charset=utf-8',
'Transfer-Encoding': 'chunked',
'Connection': 'keep-alive',
'Vary': 'Accept-Encoding, Authorization',
'Access-Control-Allow-Origin': '*',
'Access-Control-Allow-Headers': 'X-Requested-With, X-Prototype-Version, X-CSRF-Token, Authorization',
'Carto-Rate-Limit-Limit': '21',
'Carto-Rate-Limit-Remaining': '20',
'Carto-Rate-Limit-Reset': '0',
'Cache-Control': 'public,max-age=5,must-revalidate',
'Surrogate-Key': '{removed for privacy}',
'Last-Modified': 'Mon, 10 Feb 2020 20:14:40 GMT',
'X-Layergroup-Id': '{removed for privacy}',
'X-Content-Type-Options': 'nosniff',
'X-Varnish': '47747576',
'Age': '0',
'X-Cache': 'MISS',
'Content-Encoding': 'gzip'


texts: Showing only raster attribute table using GDAL?

texts:
I have a binary ArcGRID raster (AIG/Arc/Info Binary Grid) with an attribute table. I would like to get only this attribute table but gdalinfo outputs also projection information and bounds and then wanted GDALRasterAttributeTable XML part:

$ gdalinfo us_120fbfm13 -nogcp -nomd -noct filename
Driver: AIG/Arc/Info Binary Grid
...
Coordinate System is:
PROJCS["unnamed",
...
Center      ( 1540680.000, 1553400.000) ( 78d42'42.88"W, 35d46' 8.50"N)
Band 1 Block=256x4 Type=Byte, ColorInterp=Undefined
  Min=1.000 Max=99.000
  NoData Value=255
<GDALRasterAttributeTable>
...

How to show only raster attribut table using GDAL command line tools gdalinfo or perhaps gdal_translate?

Additional question which might be separate but might be also answered together with this is whether there is a possibility to get the attribute table in format such as CSV and not XML.

I'm on Linux, so I can actually do something like:

gdalinfo us_120fbfm13 -nogcp -nomd -noct -nofl | grep "<.*>"

but I would rather use some more general and possibly cross-platform solution.


texts: Grass 7 on Ubuntu 14.04. Can't start Grass

texts:
I attempted to install GRASS GIS 7 on ubuntu 14.04 using these instructions.

When I call grass70 from the command line I get a welcome message and a prompt to continue. I then get two error messages:

Traceback (most recent call last):
  File "/usr/lib/grass70/gui/wxpython/gis_set.py", line 31, in <module>
    from core import globalvar
  File "/usr/lib/grass70/gui/wxpython/core/globalvar.py", line 29, in <module>
    from core.debug import Debug
  File "/usr/lib/grass70/gui/wxpython/core/debug.py", line 25, in <module>
    import grass.script as grass
  File "/usr/lib/grass70/etc/python/grass/script/__init__.py", line 5, in <module>
    from db     import *
ImportError: No module named 'db'

which says there is an error in the GUI startup, switching to text-mode and then

File "/usr/lib/grass70/bin/grass70", line 251
    os.mkdir(tmpdir, 0700)
                    ^
SyntaxError: invalid token

which leaves me back in the bash shell.

Does anybody have some advice on what could be going wrong and how it could be fixed?

UPDATE: I installed the Grass7 stable version that has just been released using the following:

sudo add-apt-repository ppa:ubuntugis/ubuntugis-unstable
sudo add-apt-repository ppa:grass/grass-stable
sudo apt-get update
sudo apt-get install grass7

I still get the ImportError: No module named 'db'.

There is a module called db.py in the path /usr/lib/grass70/etc/python/grass/script/.

/usr/bin/python --version outputs Python 2.7.8 but python --version outputs Python 3.4.2 :: Anaconda 2.1.0 (64-bit).


texts: Slow IE8 responsiveness with OpenLayers and many features

texts:
Background

Adding about 700 sites to an OpenLayers Map. Firefox 11 and Google Chrome display the sites just fine. IE8, which is the organization's standard and will not be upgraded for the foreseeable future, bogs down when all 700 sites are displayed.

Each site has a representative icon, from a set of 15 icons. Clicking on each site displays an informative popup. The popup window is loaded on a per-click basis. The HTML for the popup window is loaded into memory once (and only once) from a single template. Clicking the site reveals the popup, which dynamically loads JSON data and fills in the template on-the-fly (to reduce the memory footprint).

This application uses OpenLayers v2.12 (723 kb).

Code

The sites are being added to the "sites" layer as follows:

/**
 * Adds a site to the list of sites.
 */
function addSite( sites, uuid, label, lon, lat, icon ) {
  sites.push(
    new OpenLayers.Feature.Vector(
      getPoint( lon, lat ),
      {
        uuid: uuid,
        siteType: icon,
        title: label
      }
    )
  );
}

/**
 * Fetches the list of all sites from the query service and adds them
 * to the site layer.
 *
 * @param siteLayer - The site is added to this layer.
 */
function addSites( siteLayer ) {
  $.ajax({
    url: MAP_SERVER_SERVICE + 'sites',
    type: 'GET',
    dataType: "json",
    contentType: "application/json; charset=utf-8",
    success: function( data ) {
      var sites = [];

      $.each( data, function( key, val ) {
        addSite(
          sites,
          val.uuid,
          val.label,
          val.longitude,
          val.latitude,
          val.site_classification );
      });

      siteLayer.addFeatures( sites );
    }
  });
}

About 100k of JSON data is returned from the map server service (i.e., the Ajax call).

Question

How would you load ~700 features into IE8 through OpenLayers without bogging it down?

The sites cannot be rendered on map tiles because the icon is dynamic (i.e., different search results will produce different tiles at the same location).

Ideas

Some ideas I have considered include:

Upgrade – Probably not feasible to ask everyone who will be using the application to upgrade to IE9.
Browser – Might not be feasible to ask people to download Chrome (or Firefox) just to view a web application.
Empty Map – This would make the map quite responsive initially, but would still bog down when showing a large number of sites. This is the most trivial of solutions.
Optimize – There is no guarantee that even after optimizing OpenLayers that the site will be more responsive.
Restrict Area – Define a region of “visibility” (like a rectangular section) and only display sites that fall within that region. This significantly reduces the usefulness of the tool.
Site Limits – Constrain the number of sites returned for any given search. This would produce misleading results.
Hide Overlapping Features - Hide sites that are partially hidden by other sites. This is probably produces optimal results. This is a type of cluster strategy.

What other possibilities are there to speed up IE8?

For what it's worth, the postDraw and createNode (30.25 seconds) methods of OpenLayers take up most of the application time, according to the IE8 debugger.

Update #1

To illustrate the amount of data:

There are sites that fall outside of the province, but they are literally few and far between.

Zoomed in view:

Update #2

I thought I could "hide" some of the features, but this does not help:

function hideFeatures( layer ) {
  var features = layer.features;

  for( var i = 0; i < features.length; i++ ) {
    features[i].style = { visibility: 'hidden' };
  }

  layer.redraw();
}

It looks like they'll actually have to be removed from the layer to be performant.

Related Links

https://stackoverflow.com/questions/6528243/how-to-hide-vector-features-in-openlayers
https://stackoverflow.com/questions/6578846/get-distance-in-pixels-between-2-lonlat-objects-in-openlayers
http://openflights.org/blog/2009/10/21/customized-openlayers-cluster-strategies/

Thank you!


texts: Moving max using calc and movingFun on a stack of 12 rasters r

texts:
I have downloaded the netcdf files of climate water deficit (def) from TerraClimate, years 1981-2010. Each raster is a monthly def for my region of interest, that being entire India.

I created a stack of 12 rasters, which is a stack of monthly def for year 1981 only and want to find the moving/rolling max value, considering a window of 3 months. I have a very straightforward single line of code-

trial<-calc(rasterstack, funtion (x) movingFun(x,3,max, circular=FALSE, na.rm=TRUE, type='around')

It works and everything. However, I do not understand how are there 12 rasters, when there should be only 10 rasters i.e. no rolling max values for 1st month Jan 1981 and last month dec 1981, since I speciefied circular=FALSE and 'around'.

Why are there 12 rasters instead of 10? Any other ways to calcualte moving/rolling max values?


texts: Limit to download images from Open Access Hub

texts: I am running a code https://github.com/diogocaribe/harpia/blob/dev/app/SentinelDownloadTile.py to download all files of Sentinel 2 from my area of interest. However, now a days the download doesn't work but the code is the same. Is there limitation of how many files we can get from Open Access Hub?


texts: QGIS 3.10, Python IDE - Application path not initialized

texts:
I installed QGIS 3.10 on Win10-64 and I try to use the qgis core and processing functionalities from PyCharm IDE.
I used to do this successfully with previous QGIS versions, but this time the output puzzles me.
I initiate my batch file like this:
@echo off
SET OSGEO4W_ROOT=C:\OSGeo4W64
call "%OSGEO4W_ROOT%"\bin\o4w_env.bat
call "%OSGEO4W_ROOT%\bin\qt5_env.bat"
call "%OSGEO4W_ROOT%\bin\py3_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\Python37\Scripts
path %PATH%;C:\qgis-plugins\qgis-otb-plugin
path %PATH%;%OSGEO4W_ROOT%\apps\qgis\python\plugins
path %PATH%;%OSGEO4W_ROOT%\apps\qgis\python\plugins\processing
path %PATH%;C:\Users\...\AppData\Roaming\QGIS\QGIS3\profiles\default\python\plugins
path %PATH%;%OSGEO4W_ROOT%\apps\saga
path %PATH%;%OSGEO4W_ROOT%\apps\gdal2

set QGIS_PREFIX_PATH=%OSGEO4W_ROOT%\apps\qgis

set GDAL_FILENAME_IS_UTF8=YES

set VSI_CACHE=TRUE
set VSI_CACHE_SIZE=1000000
set QT_PLUGIN_PATH=%OSGEO4W_ROOT%\apps\qgis\qtplugins;%OSGEO4W_ROOT%\apps\qt5\plugins

set PYTHONPATH=%PYTHONPATH%;%OSGEO4W_ROOT%\apps\qgis\python
set PYTHONHOME=%OSGEO4W_ROOT%\apps\Python37

set QT_QPA_PLATFORM_PLUGIN_PATH=%OSGEO4W_ROOT%\apps\Qt5\plugins\platforms
set QGIS_PREFIX_PATH=%OSGEO4W_ROOT%\apps\qgis

start "PyCharm aware of Quantum GIS" /B "C:\Program Files\JetBrains\PyCharm Community Edition 2018.3\bin\pycharm64.exe" %*

cmd.exe

Then I try with this little test script:
import os
import sys

from qgis.core import *

for alg in QgsApplication.processingRegistry().algorithms():
        print(alg.id(), "->", alg.displayName())

Which returns

Application path not initialized
Application path not initialized
Process finished with exit code 0

I have built the project interpreter (virtual env) through the base interpreter "C:\OSGeo4W64\apps\Python37\python.exe" inheriting global site packages.
Just a heads up, the cmd gives the following PATH appended variables:
C:\Windows\System32>echo %PATH:;=&echo.%
C:\OSGEO4~1\apps\Python37
C:\OSGEO4~1\apps\Python37\Scripts
C:\OSGEO4~1\apps\Python37
C:\OSGEO4~1\apps\Python37\Scripts
C:\OSGEO4~1\apps\qt5\bin
{app}
C:\OSGEO4~1\apps\Python27\Scripts
C:\OSGEO4~1\bin
C:\WINDOWS\system32
C:\WINDOWS
C:\WINDOWS\system32\WBem
C:\OSGEO4~1\apps\msys\bin
C:\Program Files\R\R-3.5.1\bin\x64
C:\Program Files\RStudio\bin
C:\OSGEO4~1\apps\qgis\bin
C:\OSGEO4~1\apps\grass\grass78\lib
C:\OSGEO4~1\apps\Qt5\bin
C:\OSGEO4~1\apps\Python37\Scripts
C:\qgis-plugins\qgis-otb-plugin
C:\OSGEO4~1\apps\qgis\python\plugins
C:\OSGEO4~1\apps\qgis\python\plugins\processing
C:\Users\...\AppData\Roaming\QGIS\QGIS3\profiles\default\python\plugins
C:\OSGEO4~1\apps\saga
C:\OSGEO4~1\apps\gdal2


texts: Merge multiple linestring features into a single multilinestring

texts:
Using QGis or ogr2ogr how do I convert a vector file containing multiple linestrings features into a single multilinestring feature?

At the moment I've been doing this manually by converting the data to GeoJson and editing it manually.


texts: vb.net arcobjects determine 'create features' on which layer

texts:
When I create new feature in editing session, I need to determine the new created feature is in which layer because after that have further more action on the layer. I wonder can this be done in vb.net arcobjects?

I am using vs2012 and arcmap 10.2.


texts: Contours exported from ArcGIS to DXF yield an extremely large file. Is there an alternative?

texts:
I am exporting elevation contours from ArcGIS to a DXF file so that I can FTP the result to another office. The resulting file is extremely large. When  export the same contours to a shapefile (from file geodatabase) it not only takes a fraction of the processing time but is 6x smaller in size then the DXF file. This is only the major contours, I still need to create a DXF file from the minor contours and the process is crashing all together.
I need to deliver the contours in DXF format.  I tried compressing the file but it is still too large. Is there another way?


texts: Returning current attribute name in QGIS Expression Dialog

texts:
In Expression Dialog, I can use $geometry to returns the geometry of the current feature.
But how could I return the current attribute of the current feature? Something like $attribute.
The purpose of this question is to save the time finding the attribute name form "fields and values" list and click the attribute name multiple times.


texts:
I want to perform a view shed analysis in QGIS using GRASS 7 r.viewshed.

When I double click on r.viewshed in GRASS GIS 7 commands in the processing toolbox of QGIS I get the following message:

  Missing dependency.This algorithm cannot be run :-(  This algorithm requires GRASS GIS 7 to be run. Unfortunately, it seems that GRASS GIS 7 is not installed in your system, or it is not correctly configured to be used from QGIS Click here to know more about how to install and configure GRASS GIS 7 to be used with QGIS

GRASS commands seem to work such as r.los but the GRASS 7 commands don't.

I am using a MAC OS X Mavericks with the latest version of QGIS 2.4.0 Chugiak.


texts: Generating wind direction raster using its U and V components in ArcGIS Desktop?

texts:
I have two raster layers representing the U and V components of a vector wind. I could generate the magnitude raster by applying (u*u+v*v) in Map Algebra in ArcGIS.

Now I need to generate the direction raster with the following function:

if(u > 0 && v > 0,( atan(u/v) + 180),if(u < 0 && v < 0,(atan(u/v) + 0),if(u > 0 && v < 0,( atan(u/v) + 360), if(u < 0 && v > 0,( atan(u/v)+270)))))

but I don't know how to apply it in ArcGIS. Any help please?


texts: Best practices for planning hardware requirements for ArcGIS Server deployment scenario

texts:
In this question I would like to ask for sharing of some experience in planning ArcGIS Server hardware part in Web mapping projects (Flex, Silverlight).

Let suppose that we are developing RIA application and deploying ArcGIS Server with services.

I have noticed that very often on the stage of early planning of GIS projects customers would like to know approximate hardware requirements for ArcGIS Server part.

Usually it sounds in these ways:

We need to serve N editor user and M viewers of spatial data. What server do we need (CPU, RAM)? How many ArcGIS servers?
We already have dedicated server which we are planning to use for ArcGIS Server. How many users (editors/viewers) will be it capable to serve?

It is common that on this stage you don't have access to spatial data so it is not possible to perform any tests. But we can estimate rough amount of Map/Feature and Dynamic/Cached services.

I have studied Esri System Design Strategies document but I have not found any numbers or formulas to make such rough estimation. As I understand the best precise way is to obtain at least part of the data and perform real tests.

So, what is the best practices in such situation? Is it possible to perform such very rough calculations to answer these questions without tests with real data?

May be you already faced such problems and have some experience?


texts:
I have a problem with proxy.cgi. I read other posts but didn't find a solution. When I make my WFS request in console, the following error comes up: "POST http:// localhost:8081/cgi-bin/proxy.cgi?url=http%3A%2F%2Flocalhost%3A8080%2Fgeoserver%2Fwfs 404 (Not Found)", in other words the application don't find my proxy.cgi file. My settings are:

OpenLayers.ProxyHost="/cgi-bin/proxy.cgi?url=" in OpenLayers.js;
uncomment tag servlet and servlet-mapping in web.xml file of server Tomcat;
add privileged="true" attribute to tag Context in context.xml of server Tomcat;
change first row of proxy.cgi in #!C:/Python27/python.exe (installation path on Windows 7).

Any idea how to resolve the problem?


texts: How to Change Attribution Parameter of WMS Layer Dynamically in Openlayers

texts:
This might be a noob question. I have a wms layer in OpenLayers 2.13.1 with attribution: "some_text" set when it's added. This is working fine, the attribution shows up fine in the map div.

My question is how to change the content of the attribution param when the layer is already added. So far I have tried myWMS_layer.mergeNewParams({attribution: "some_updated_text"}); but no luck. I've also find updateAttribution function in OpenLayers docs, but don't really know how to use it. I've tried myAttributionControl.updateAttribution(); and nothing happened.

Thanks in advance!