texts:
I'm trying to process two raster on the Raster Calculator of PyQgis and I wrote a script, but QGis 3.2 crashes everytime I run the script.

import processing
import os
import time
from qgis.analysis import QgsRasterCalculator, QgsRasterCalculatorEntry
folder="C:/++ELABORAZIONI_NUOVE/entrevistas/PUNTOS_ENTREVISTA_PROVA/mask/"
filelist=os.listdir(folder)
entries=[]
feedback = QgsProcessingFeedback()
for file in filelist:
    if file.endswith('.tif'):
        layer=QgsRasterLayer('C:/++ELABORAZIONI_NUOVE/entrevistas/PUNTOS_ENTREVISTA_PROVA/mask/','layer')
        lyr=QgsRasterCalculatorEntry()
        lyr.ref=layer.name()+'@1'
        lyr.raster=layer
        lyr.bandNumber=1
        entries.append(lyr)
        comarca=QgsRasterLayer('C:/++ELABORAZIONI_NUOVE/expert_panel/comarca_alt_emporda.tif','comarca')
        cmc=QgsRasterCalculatorEntry()
        cmc.ref=comarca.name()+'@1'
        cmc.raster=comarca
        cmc.banNumber=1
        entries.append(cmc)
        expression='(('+ entries[1].ref +') = 1)*('+ entries[0].ref +')' +'+ ((('+entries[1].ref+') = 0) - 999)'
        print (expression)
        folder2='C:/++ELABORAZIONI_NUOVE/entrevistas/PUNTOS_ENTREVISTA_PROVA/CALC/'
        calc = QgsRasterCalculator(expression, folder2 + file, 'GTiff', layer.extent(), layer.width(), layer.height(), entries )
        calc.processCalculation()


texts: Filter by bound on Google Earth Engine not covering the whole of bound feature

texts:
I tried to filter Sentinel-2 image collection by bound, then I visualized  the first in the collection after clipping to the bound feature.
The problem is that the returned images does not cover the whole bound feature, probably because the feature extends beyond one image row.
Is there a way to fix this? maybe with composites?
This is the code I used:
var start = ee.Date('2020-7-5');
var end = start.advance(1,'month')
var dateRange = ee.DateRange(start,end)

var s2 = ee.ImageCollection('COPERNICUS/S2_SR')
                          .filterBounds(table)
                          .filterDate(dateRange)
                          .sort('CLOUDY_PIXEL_PERCENTAGE');
var img = s2.first() ;
print(img);
var param = {'bands':['B4', 'B3','B2'], 'min':400 , 'max':3500 }
Map.centerObject(table,12);
Map.addLayer(img.clip(table), param, 'se img');
Map.addLayer(table,{},'table')

The map window


texts:
Esri has a base map available through a web service:

http://services.arcgisonline.com/ArcGIS/rest/services/World_Street_Map/MapServer

Do Google or Bing, or anybody else have a similar service?


texts: Creating a variable distance buffer in QGIS

texts:
I have a QGIS project with a shapefile (CRS:WGS84) overlayed on a googlemaps satellite layer.
My shapefile consists of many GPS coordinate points measured in decimal degrees, and another field for accuracy of the GPS reading, measured in meters.
I want to create a variable distance buffer around the points so that each point has a buffer the size of the respective accuracy reading.
For example if a specific point has an accuracy reading of 6.0 meters, I want a 6.0 meter buffer around that point.

I've tried creating a variable distance buffer normally but it renders totally off. I think this is probably due to the lat/long an CRS being set for decimal degrees, but the accuracy measurement in meters. However I have no easy way of transforming that accuracy into decimal degrees without transforming the entire map.

Is there a way to create a variable distance buffer in meters around points in decimal degrees?


texts: Attaching pdf files to qgis2web webmap accessible through a server

texts:
I have made a webmap using QGIS2Web.
I am able to attach pdf documents using    the line

google report

When the webmap source files (shown in the image below) are on a local drive,it is possible to access the pdf. However, when i move them to a server, i can't access the pdf files.

Is it possible to have them in a folder in similar location to that of the source files and if so how do i edit the code to read them?


texts: How to get rid of "demo" from geowebcache layer link

texts:
I wan to use standalone geowebcache. I know how to serve cached layers but one thing really annoys me and it is the "demo" in the link e.g.

http://myserver.com/geowebcache/demo/layer

I would like to have the link in form without the "demo" or to be able to change it to something else.


texts:
This had been asked previously on Stack Overflow.

I've tried different methods/formulas but results are not correct. 40.764296,-73.97302 should result into 930,3590 pixel on map image. I've Manhattan image map rotated 28.34. I could not find correct formula for converting lat/lng to pixel and vice versa. I May be missing something important?

Manhattan map image (1816x8160) having following lat/lon of corners.

TopLeft: (-73.9308,40.8883)
TopRight: (-73.8584,40.858)
BottomLeft: (-74.0665,40.7024)
BottomRight: (-73.9944,40.6718)
Map is not true north and rotated at 28.34, also its UTM Zone 18N (78W to 72W). Here are further details about this map taken from PDF Maps iOS app.

Size (pixels): 1816 x 6160
Pixel Resolution: 3.829 meters
Bounds (pixels): (-1624, -3518) x (7866, 7719)
PROJCS["WGS 84 / UTM zone 18N",
    GEOGCS["WGS 84",
        DATUM["WGS_1984",
            SPHEROID["WGS 84",6378137,298.257223563,
                AUTHORITY["EPSG","7030"]],
            TOWGS84[0,0,0,0,0,0,0],
            AUTHORITY["EPSG","6326"]],
        PRIMEM["Greenwich",0,
            AUTHORITY["EPSG","8901"]],
        UNIT["degree",0.0174532925199433,
            AUTHORITY["EPSG","9122"]],
        AUTHORITY["EPSG","4326"]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",-75],
    PARAMETER["scale_factor",0.9996],
    PARAMETER["false_easting",500000],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AXIS["Easting",EAST],
    AXIS["Northing",NORTH],
    AUTHORITY["EPSG","32618"]]

How to convert lat/lon to x y and vice versa?


texts: Use custom Transverse Mercator projection for poles

texts:
Due to some reason (internal storage system), I want to define a custom Transverse Mercator projection for some parts of the world where the origin is nearby an area of interest.

But I also need to cover the poles, and as I said I want to stick to Transverse Mercator projection - if possible.

Is there any reason why I should use another projection like UPS instead?


texts: How can I monitor the open ArcSDE database connections without sdemon? I need to see at any time how many sessions are consumed from my app and if they are related to ArcGIS. I use Oracle 10.2.


texts: Perform a group by query on geodatabase in ArcMap

texts:
I have a geodatabase called predictions.gdb and a feature class called output_gis_120_Project in that database. I imported this feature class as a layer in ArcMap. The attribute table with columns lat,lng,attrx,attry,avg_price which is in the form

>

lat1,lng1,attr1,attr2,avg_price1
lat1,lng1,atrr3,attr4,avg_price2

etc... I would like to perform a sql query which is

  SELECT lat,lng,SUM(avg_price) from output_gis_120_Project GROUP BY
  lat,lng;

I read that ArcGIS can only have group by statements in a subquery so I went to layer--->properties-->definition-->query builder and rewrote my query as

  "select * from output_gis_120_Project where lat,lng,sum_price in
  (SELECT lat,lng,SUM(avg_price) sum_price from output_gis_120_Project
  GROUP BY lat,lng)

I keep on getting an error saying that there was an error with the expression. I tried putting the fields in brackets [] and only selecting avg_price in the subquery, but this didn't help either.


texts: How to Close all Layers but keep Tables on MapInfo using MapBasic

texts: I'd like to know how to close all the layers on MapInfo but keep the tables open using MapBasic code. The idea is to be able to keep the tables open so that I can create multiple different maps using the same tables. Basically, for every map I create and export as an image , I want to close the layers for that map and then start creating my next map with the same tables as before. I don't want to Close all and then reopen the tables.


texts: Georeference raster in WGS84 and project in North azimuthal in QGIS

texts:
I have a raster (JEPG) image from a full planet surface (in this case Jupiter).
I use QGIS and its plugin Georeferencer to add control points on the image. I use the four corners. Then I use Gdal to warp the new georeferenced image into north azimuthal projection. But the final image isn't a circle at all! So I probably missed something during the process.
Here is the full process I followed:

Open the raster in Georeferencer QGIS's plugin
Set the coordinate reference system as WGS84 (EPSG:4326)
Select four points (corners) of the raster and set the destination value:

Select upper left corner and set X/Y to -180/90
Select upper right corner and set X/Y to 180/90
Select lower left corner and set X/Y to -180/-90
Select upper right corner and set X/Y to -180/-90

Configure the transformation setting to the following parameters:

Transformation Type: Polynomial 1
Resampling method: cubic
Compression: None

Georeference image.
Warp into azimuthal projectino with gdarwarp :
gdalwarp -t_srs "+proj=aeqd +lat_0=90 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs" jupiter_georef.tif jupiter_circle.tiff

Unfortunately, the final image doesn't match a full circle.

The georeferencer screen :
Output:

I also tried to use QGIS to warp the georeferenced image by changing the coordinate system. The output is better but the circle is cropped by a square and it has a thin line missing.

Any clue to solve this problem?


texts: Visualizing various layers through commmand line/automation on QGIS

texts: I am new to QGIS but I would like to visualise spatial data. I have a huge amount of layers and I want to take a snapshot of all these different layers with a script. I was hoping I could select which layer I want and maybe configure something from each layer through the command line but I am not so sure how to go about it. Actually, I am not so sure what tools QGIS offers for this task or if it offers any at all? Maybe I should look at other tools than QGIS?


texts: ArcGIS Map Export - Series of Attributes

texts:
I am working with ArcGIS 10.7 and can not find a solution to my problem.

I have got a shapefile with 1000 datasets and 20 different species of plants. Now I want to export the 20 species maps in series.
It is always the same extent just records of the species vary.

I worked with Data Driven Pages a lot but I just cant see how this has to be done.

Any ideas?


texts: How do I get the maximum value of raster cells ignoring nodata with gdal_calc?

texts:
I have two integer rasters, A and B, with the same dimensions and cell sizes. I want to compare each cell in A with each cell in B and keep the maximum value. In cases where both have nodata, I want the result to be nodata. In cases where one is nodata, I want the result to be the integer value of the one with data.

   Raster A                    Raster B
|   1  |   1  |             |   3  |   1  |
|   1  |nodata|             |nodata|   3  |
|nodata|   1  |             |nodata|nodata|
|nodata|   4  |             |nodata|   1  |

--------------------------------------------
                Desired result
               |   3  |   1  |
               |   1  |   3  |
               |nodata|   1  |
               |nodata|   4  |

've tried combining them in many different ways using gdal_calc but somehow I can't seem to get the answer I'm looking for. My first instinct was to use numpy.fmax, which states that the NaN value is ignored. But this calculation

--calc "numpy.fmax(A,B)"

only gives me the results I want in cells where both datasets have a value. A nodata in either results in nodata in the result:

              numpy.fmax result
               |   3  |   1  |
               |nodata|nodata|
               |nodata|nodata|
               |nodata|   4  |

Perhaps I'm confused about how nodata translates to NaN in numpy? What is a straightforward way to accomplish this with gdal_calc?


texts:
I am flummoxed on how to use GRASS's v.kernel.

I have a vector layer of around 2.5 million points. I want to make a heat map using v.kernel to show concentrations, since I have variable instances with overlapping points, sometimes huge overlaps.

I've already gotten this vector layer in GRASS, and it displays just fine.

I've tried using GRASS's v.kernel command based on what I've seen here and on other forums, and I can't get it to do anything besides output a raster that's just a pink square.

Here's the command I'm using:

v.kernel --verbose input=master_grass7 output=master_grass7a_heatmap stddeviation=.0001

I've varied the stddeviation to all sorts of values from 1000000 to .000001, and it had no effect.

I've read the v.kernel documentation repeatedly and don't really understand what it's getting at. At least, the instructions are on esoteric concepts, nothing practical. I've also checked the source code, and I'm not really understanding it, either. Yes, I can read C. The problem is it depends on a lot of stuff defined elsewhere in GRASS GIS.

I've also done a lot of Google searching, and I can't find a comprehensive guide. All that I'm getting are scattered copies of the v.kernel doc/man page or people who apparently got it to work without a fuss.

I've also checked up on the concept of kernel density estimation (KDE), and even then I don't see how to use the v.kernel command. That command appears to be a specific interpretation of KDE; its switches don't appear to correspond well to generic KDE concepts.

So back to the main question here: how can someone who is not intimate with GRASS product development use the v.kernel command? Is there a plain language translation available?


texts:
Hi I am beginner at QGIS. I am trying to export the attributes "surface" and "parking lane" of every street of the city Dresden.

I tried it with "overpassturbo", but it doesn't give me every street, only few. If i try to export with the OSMExport command, i have to minimize my area, it doesn't work for the whole area of the city.


texts: Generating elevation 'database' from NED IMG data

texts:
I'm starting to work with the NED 1/3 and 1/9 arc-second data. I've been using gdal in Python and QGIS to do some basic conversions of the raw .IMG data from USGS. My goal is to get elevation data of the San Francisco Bay Area that I can look up with coordinates. It seemed like the most straightforward way to do this was to use gdal to convert to .XYZ format - this resulted in a really big file from a small piece of the area I'm using as a test.

Eventually I'd like to do this for the state of California, the USA, etc. so scalability is a concern. Although I envision some sort of database, if gleaning the info directly from the IMG or some other format is better I'm open to the idea.

A similar poster reads data from the files directly - would this be a good route to go from a memory/speed point of view? What are the tradeoffs? If this is a better method than converting to raw data directly, can anyone point me to some resources on how to accomplish this in python?

What is the best way to get to where I want to be?


texts: Interpolated raster not aligned with the vector data

texts:
I have loaded some tab delimited text data into QGis to interpolate the points. When I get to the interpolation (IDW) part everything breaks down. The output raster file does not align with the original data. The difference in alignment differ depending on how many rows and columns I choose for the output file. If I choose the standard 300 x 300 the raster is south of the original data and if I choose like 400 x 400 it goes north of the original data..

I'm a bit new to QGIS, perhaps I'm not checking a box somewhere or something!?


texts: Arcmap: Unable to detect empty value in row using Field Calculator

texts:
I am using the following pre-logic Python code to try to catch Null, empty, or "0" string values in the "street" column and then return a "V" value if it  is one of those. While the code runs and works on some rows, it doesn't seem to be catching all empty cells. Am I missing a case? Or could it be that there are multiple spaces, and if so, how can I check for those?

Pre-Logic:

def check(Street):
  if Street is None or Street == "" or Street == "0" or Street == None or not Street or Street == " ":
    return "V"
  else:
    return " "

Verfied =

check(!Street1!)


texts: Traceback (most recent call last): - trying to run script to generate site boundaries from excel spreadsheet

texts:
I am trying to run a script to create a polygon layer in a gdb and keep coming up with this error -

Traceback (most recent call last):
  File "I:\GIS Data\ArcToolBox\SiteBoundaryCreation_FMG_Format.py", line 87, in  polygon.created_date = surveyStartDate < --- this is the line the error refers to

File "d:\programs\arcgis\desktop10.3\arcpy\arcpy\arcobjects_base.py",
line 35, in __setattr__return setattr(self._arc_object, attr, ao)

RuntimeError: ERROR 999999: Error executing function.

Failed to execute (SiteBoundaryFMG).
Failed at Tue Aug 28 16:12:43 2018 (Elapsed Time: 3.35 seconds)

The following section of code is where line 87 sits -

#new variable Polygon.  Set it to create a new row in variable areacur (defined on Line43) at the cursor point (set in areacur)
polygon = areacur.newRow()
#set the shape property of the polygon variable to pntarray (which was set on line 107)
polygon.shape = pntarray
#set the various field to variable defined and assigned above
polygon.Site_Name = siteID
polygon.Site_Type = primT
polygon.Date_created = surveyStartDate **<-- this is the line the error refers to**
polygon.Project = projectName
polygon.Survey_ID = surveyIDName
polygon.created_user = userName
polygon.created_date = userCreatedDate

is it something to do with the name in the line?  or is it unrelated? the line the error is referring to is an optional field when using the tool in the arcmap toolbox (see attached image - Survey Start Date (Optional)). would this have something to do with the error?


texts: What is the complexity of neigborhood search in PostGIS?

texts:
I have a point with coordinates x and y. I want to know all the points in a square around this point, with side L. For this, in postgresql/PostGIS, I use:

SELECT * FROM table_of_points WHERE (ST_Transform(the_geom,srid) && ST_MakeEnvelope(" + str(x-L) + ", " + str(y-L) + ", " + str(x+L) + ", " + str(y+L) + ", srid )

with srid being the Spatial Reference System Identifier.

What is the computational complexity of this query in PostGIS?


texts: Trying to visualize a geojson on react-leaflet

texts:
I'm trying to visualize a GeoJSON layer on Leaflet using React-Leaflet library.
My mental process is:

Click on the button which calls a method named getHospitals()
getHospitals() gets the geojson layer from the backend using axios
The geojson layer is saved in a variable
Then, create a reference which will passed to the geojson component in order to 'update' the component (I don't know if it's correct)
Update the component state with the layer and the reference (this.setState())
The component renders the new geojson

The problem is the geojson doesn't appear on the map. This is the code:
import React, { Component } from "react";
import {
  MapContainer,
  TileLayer,
  LayersControl,
  ScaleControl,
  GeoJSON,
} from "react-leaflet";
import "leaflet/dist/leaflet.css";
import axios from "axios";

import "./Mapa.css";

export class Mapa extends Component {
  state = {
    geoJSON: null,
    reference: null,
  };

  /* Definición de estilos del mapa */
  style = {
    height: 600,
    margin: 10,
  };

  getHospitals = async () => {
    let hospitals = await axios.get(
      "https://cartovis-server.herokuapp.com/hospitales"
    );

    let reference = React.createRef();

    this.setState({
      geoJSON: hospitals.data,
      reference: reference
    });
  };

  render() {
    return (

        <MapContainer center={[38, -1]} zoom="5" style={this.style}>
          <GeoJSON
            attribution="Capa de Hospitales de ESRI"
            data={this.state.geoJSON}
            ref={this.reference}
          />
          <ScaleControl />
          <LayersControl>
            <LayersControl.BaseLayer checked name="OpenStreetMap">
              <TileLayer
                attribution='&copy; OpenStreetMap contributors'
                url="https://{s}.tile.openstreetmap.org/{z}/{x}/{y}.png"
              />
            </LayersControl.BaseLayer>
          </LayersControl>
        </MapContainer>

Hospitales

          <button className="btn btn-primary" onClick={this.getHospitals}>
            Ver
          </button>

    );
  }
}

export default Mapa;

Any idea?


texts: Changing transformation of already reprojected raster in ArcGIS Desktop?

texts:
I am pretty confused by projecting data from one coordinate system to another, especially if multiple transformation options are available (ArcGIS 10.4).

I downloaded Landsat LT5 data (origCS = GSC: WGS84, EPSG:4326) and I reprojected it to final PCS: S-JTSK (desiredCS = EPSG: 5514). Multiple transformations are available in ArcGIS:

I therefore applied the NGA transformation (S_JTSK_To_WGS_1984_NGA, called _nga) and further processed and classified my rasters.

Now,  I am thinking to change the way how Landsat has been projected, i.e change transformation #1 to transformation #2, keeping the original (WGS84) and final coordinate systems (S-JTSK) the same.

If that is possible and "legal", what is correct approach to do that?

As processing of my projected rasters has taken lot of time, I don't want to start again with the raw Landsat data, and I prefer to play around with transformation of classified rasters.

My thoughts how to change transformations of final rasters are:

simply transform the raster_nga to raster_t2 (from one transformation to another, staying in S-JTSK) - is this possible?
project one landsat_wgs to lansat_t2 (by transformation 2) and use it as master scene  to resample final classified rasters?
something else??
is this needed or even not at all??

I am working on study area of size 400 km2, using single Lansat scene per time.


texts: How to create and export complex symbology?

texts:
I've got punctual geometries that I'd like to represent in a desktop and a web application, using complex symbology, like this :

My idea is to use a symbology editor to create a SLD file that defines the symbols.
I used QGis to create the symbology and then tried saving it into a SLD file. The problem is that symbols are simplified and instead of the symbol I created, I get this one:

Do you know about a SLD editor that supports complex symbology?

Do you have any other suggestion different from creating a SLD file?


texts:
I'm using geotools v14.3 to load a shapefile (Java Swing application). I want to allow the customers to use the mouse to click on the map area to draw new points, lines and polygons. First, I'm testing the points. So, a customer should be able to draw a point on the clicked location.

My current approach

After following few articles, I first capture the clicked mouse location with a mapmouselistener.

@Override
    public void onMouseClicked(MapMouseEvent ev) {
        System.out.println("mouse clicked : " + ev.getWorldPos());
        DirectPosition2D p = ev.getWorldPos();
        drawMyPoint(p.getX(), p.getY());

    }

Point drawing function is as follows.

void drawMyPoint(double x, double y) {
        SimpleFeatureTypeBuilder b = new SimpleFeatureTypeBuilder();
        b.setName("MyFeatureType");

        // set crs
        b.setCRS(DefaultGeographicCRS.WGS84);
        // add geometry
        b.add("location", Point.class);

        // build the type
        final SimpleFeatureType TYPE = b.buildFeatureType();

        // create features using the type defined
        SimpleFeatureBuilder featureBuilder = new SimpleFeatureBuilder(TYPE);
        GeometryFactory geometryFactory = JTSFactoryFinder.getGeometryFactory();
        com.vividsolutions.jts.geom.Point point = geometryFactory.createPoint(new Coordinate(x, y));
        featureBuilder.add(point);
        SimpleFeature feature = featureBuilder.buildFeature("FeaturePoint");
        DefaultFeatureCollection featureCollection = new DefaultFeatureCollection("internal", TYPE);
        featureCollection.add(feature); //Add feature 1
        // Add feature 2, 3, etc

        Style style = SLD.createPointStyle("Star", Color.BLUE, Color.BLUE, 0.3f, 15);
        Layer layer = new FeatureLayer(featureCollection, style);
        layer.setTitle("NewPoint");
        map.addLayer(layer);
        mapPane.repaint();

    }

Problem

I can see the layer name being added to the map layer table. But the point is not visible on the map.


texts: Creating world files for a batch of aerial photos

texts:
I am looking for a way to geolocate a bunch of aerial photos by creating a world file per photo. The photos are taken in the same orientation and from the same altitude.
The work flow would be to create one world file with the correct first 4 lines of code giving the correct rotation and scale. And then to calculate all the top left corner points for all the photos via a formula.

To create a world file for a aerial photo it first needs to be correctly referenced. I thought of referencing 1 photo and then making copies of the world file for all the othe photos. Except after I georeferenced the first photo and exported the first world file I realized the gps point of the world file is not on the top left corner of the photo, instead it is slightly off to the side, this is because the photo was rotated so the program used the top left corner of the geotiff where there is no photo data just a black triangle. I cannot use these values form the world file because the reference is off.

Is there a way of exporting a world file by first specifying where the gps point will be and that the top left pixel of the photo will also be on that point or a formula that calculates lines 1,2,3 and 4 of the world file for me?


texts: Set symbols for a layer in ArcMap with two classes (greater than zero and less than zero)

texts:
I need to symbolize features in a feature class with two colors thus defining what features have positive (>0) and negative (<0) values. The only two ways I found this to be possible are:

To use a static field with the values calculated beforehand

OR

Use definition queries using two identical layers in the TOC.

Since it is more time consuming to update labels for two layers, I have approached this problem in another way by using Graduated colors symbol setting with two classes (Manual classification). I have added two extra features in the feature classes with really big and really small values (to cover possible values that can come later on while the feature class will be edited) and then edited the intervals in the Range columns in the Symbology tab of the Layer Properties. Afterwards, those two features are deleted.

Now features' values will fall within these two ranges and will be symbolized properly.

Considering the need to update multiple documents with these settings and having different datasets as layers, are there any other approaches (without using Python to automate the hack described above)? All other GIS.SE posts I've found so far refer to the methods outlined above.


texts: Coordinate System (Malaysia RSO Projection) not lining up in QGIS and ArcGIS for Desktop?

texts:
I am using QGIS 1.9.0 on Linux Mint; and I am new to (Q)GIS

I have three layers (two shapefiles, one OpenLayers layer [Bing Maps]) that I am trying to line up located in peninsular Malaysia.

Shapefile_1 is a sub-district border (WGS84; EPSG:4326)
Shapefile_2 is river system within the sub-district (GDM2000 Peninsula RSO; EPSG:3375)
OpenLayers (set as the ProjectLayer CRS; WGS84 Pseudo Mercator; EPSG:3857)

When importing the layers in QGIS, Shapefile_1 is close to the borders of the sub-district shown in the OpenLayers layer -- but clearly not matched.

Shapefile_2 is incorrectly located hundreds of kilometres north east off the coast.

When the same layers, with the (apparently) same CRS information are displayed in ArcGIS, they line up correctly.

Any insights would be greatly appreciated.


texts: Topology Editing using ArcView (now called Basic) license level of ArcGIS Desktop?

texts:
I am building a geodatabase from data collected on a GPS handheld in the field.

The GPS device collects decently accurate information but I need to make sure all the polygons, lines, and points are aligning correctly.

I currently have the an Arcview license and don't believe I can create my own topology.

Is there any work around or extension that I can download, instead of upgrading my ArcGIS license to create and edit the topology in this way?


texts: Project (lat,lng) GPS point to a line segment (road segment)

texts:
Is there a way to project a point A to a B,C line segment and find a point D ? The points are latitude, longitude points, so I believe we should incorporate the haversine distance formula.

Edit: I am trying to do something line this, but with GPS coordinates.


texts: Performance issues with (clustered) point features

texts:
How do I get a point features layer containing roughly 2800 labeled features to perform smoothly when exported (and converted to GeoJSON) with qgis2web?

Specifics: I use qgis2web 3.4.0 and QGIS 3.2.3 to create web maps and am experiencing performance issues with the following map.

    http://www.rv-k.de/LK_Esslingen/Fahrradwegweisung/Abstimmung/Online.html

Things I have tried:

File size Overall file sizes are not very large. Still, I tried to further reduce file size by changing export precision settings. No performance gains.
Number of features per layerThe layer "Knoten" contains roughly 2800 point features and seems to cause the performance issues because broken up into 5 different layers (with 700 features per layer max) and turning off labels for each layer, the map performs fine. (Other layers contain significantly less features.)

    http://www.rv-k.de/LK_Esslingen/Fahrradwegweisung/TestGIS/5Layerstest/5LayersGIS.html

Clustering features
Still, I need all features to display as one layer and with labels. Therefore, I tried clustering using https://github.com/Leaflet/Leaflet.markercluster. There are still performance issues because all 2800 features are loaded at the same time once you go beyond zoom level 16. Still, I need the markers to display individually at this point and therefore set (disableClusteringAtZoom: 16).

    http://www.rv-k.de/LK_Esslingen/Fahrradwegweisung/TestGIS/Clustertest/WebGIS.html


texts: How do I edit a table in GRASS?  When importing a vector into GRASS from QGIS, all the attributes are lost and I need these to continue with my analysis.


texts: Using Spherical Standard Deviation to identify ledges on steep cliffs?

texts:
I have been implementing the approach highlighted here: Finding regions or zones in raster DEM of similar slope or aspect values in ArcGIS for Desktop?

to calculate firstly the "spherical variance" and then the "spherical standard deviation (SD)" on a 50 cm DEM of an island.

My intention is to identify the ledges that may exist on the cliffs of the island and the relatively flat areas.

However, I have found it difficult to identify and justify how to classify the resulted spherical SD values. In particular, I understand that the closer the values are to zero, the more likely the variation of slopes in these areas is smaller (i.e. gentle and flat areas). In my understanding, I could classify such areas as the ledges I am looking for. However, I am unsure up to what values of spherical SD ledges can be classified (i.e. SD=0.1 or SD=0.2 and so on).

I seek answers pointing me to the threshold that I could set in my pursuit of classifying spherical SD values as flat areas and ledges. Ideally, a book or a scientific paper that has done something similar would be perfect.


texts: Set project variable to a field value when iterating over vector features with graphical modeler in QGIS 3

texts:
I have a Data layer containing features with some fields (e.g. rows with color and name). In the Graphical Modeler I want to iterate over the features and assign the field values to some project variables. With the Vector Features Input it should be possible to iterate over features but I couldn't find out how to read the actual field values and assign it to a variable.
Eventually I want to batch-process many image exports in different variants with the print layout. I created a single layout that can be modified with variables, e.g. it sets colors via expressions (assigning project variables).
How can I access a field value when using Vector Features as iterator?


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

texts:
I am looking for a way - given an IWorkspace object - to check if it was created from a Personal-GDB or a File-GDB. I tried using IWorkspace.WorkspaceFactory to check if it's an instance of e.g. AccessWorkspaceFactory but unfortunately this doesn't work for the fgdb. According to .NET the fgdb workspace was created by an AccessWorkspaceFactory, too. Duh.
So far I've only come up with the idea that one could check if it's an pgdb by trying to create the workspace using the according factory. The same goes for the fgdb, obviously. Like so:

try {
   IWorkspaceFactory factory = new AccessWorkspaceFactoryClass();
   factory.OpenFromFile(workspace.PathName);
   // if we got that far it seems to be a personal gdb
   pgdb = true;
} catch (Exception) {
   pgdb = false; // nope, not a personal gdb (or valid path or whatever)
}

But this solution doesn't seem to be very elegant. Are there any data structures to check where the workspace came from?


texts: Out of Range error when creating thumbnails of MXDs

texts:
I am tryting to create thumbnails of all the mxd's in a directory but no mater how high I set the value of the "While" loop it gives me an "Index out of range" error...

Also can I get it to look at subdirectories as well?

best,

#Make thumbnails from all mxd's in current directory

import arcpy, glob, os, sys
from arcpy import env
from arcpy import mapping
env.overwriteOutput = True

#path = sys.argv[1] # for using in toolset, use directory parameter
path = os.getcwd() # Script in same directory as files being processed
mxdList = glob.glob(path + "\*.mxd")
env.workspace = path
x=0

print 'Reading mxd files from ' + path

mxd_List = glob.glob('*.mxd')
count_Mapdocs = len(mxd_List)
print mxd_List,count_Mapdocs

while x < 50:
 mapDocs=mxd_List[x]
 mapDoc=arcpy.mapping.MapDocument(mapDocs)
 framesList=arcpy.mapping.ListDataFrames(mapDoc)
 mapDoc.activeView=framesList[x]
 mapDoc.makeThumbnail()
 mapDoc.save()
 print "Created thumbnail for " + mxd_List[x]
 print "x= "+ str(x)
 x=x+1
 print "x updated =" + str(x)


texts: How to add MS SQL table using python console?

texts: I have MS SQL Server 2008 R2 database. How to add MS SQL Table using python console in QGIS?


texts: Computing x and y distances between two raster maps using GRASS GIS?

texts: I am working with a maximum water elevation raster in which, similarly to using r.grow.distance, I am trying to extract the distance from the each null cell to the nearest water elevation cell. However I need both the x and y distances, which are not provided by r.grow.distance. I looked into the source code for this function and found that the x and y distances are calculated within the code but r.grow.distance can only provide results in euclidean, squared, maximum, manhattan, or geodesic formats. Is there an easy way to extract these x and y distance values? If so, I am working on a real-time storm surge model so it would be ideal to be able to calculate these values in as little time as possible. I also looked into r.distance but this does not seem like an option that would work for me.


texts: Generating regular sample points in a polygon

texts:
Is there a way I can generate regularly distributed sample points within a polygon using ArcGIS 10?  I originally used Hawth's Tools to do this in 9.3, however I don't think this is possible in 10.

Thanks!


texts:
I'm trying to figure out how closely the QGIS browser matches ArcCatalog functionality, and I am confused about whether it is supposed to be possible to move copy or move data files using the browser.
This page on the Browser panel (https://docs.qgis.org/3.10/en/docs/user_manual/introduction/browser.html) says,

The QGIS Browser panel is a great tool for browsing, searching, inspecting, copying and loading QGIS resources. Only resources that QGIS knows how to handle are shown in the browser.

which seems to suggest that copying should be possible, but when I attempt to drag-and-drop shapefiles and similar, the mouse cursor turns into a circle and slash icon when I hover over a possible destination directory, and right-clicking only seems to offer options for "Add selected layer to Canvas", "Properties...", and "Delete file..."
Am I missing something, or is does the QGIS Browser simply not offer the ability to copy files?


texts: replace ü/ ö/ ä / spaces and brackets in Attribute Table

texts:
I've got a Polygon Feature Class of communities in switzerland. Unfortunately, the names in the attribute Field "NAME" include spaces, brackets and ä/ ö/ ü characters. For further processing, i need to remove the spaces and brackets and to replace ä --> ae, ö --> oe, ü --> ue.
I've tried to do this using update cursor but i don't get it to work:

with arcpy.da.UpdateCursor("AVR_GEM","NAME") as cursor:
    for row in cursor:
        row = str(row).replace("ä","ae")
        row = str(row).replace("ö","oe")
        row = str(row).replace("ü","oe")
        row = str(row).replace(" ","")
        row = str(row).replace("(","")
        cursor.updateRow(row)

Can anyone show me a quick way to get the job done?


texts: How does the parameter "offset" work with the ST_LocateAlong/ST_LocateBetween PostGIS function?

texts:
I try to understand the ST_LocateAlong/ST_LocateBetween PostGIS function. The documentation tells me this:

  geometry ST_LocateAlong(geometry ageom_with_measure, float8 a_measure,
  float8 offset);

  Return a derived geometry collection value with
  elements that match the specified measure. Polygonal elements are not
  supported.

  If an offset is provided, the resultant will be offset to the left or
  right of the input line by the specified number of units. A positive
  offset will be to the left, and a negative one to the right.

I understand that a point will be returned if feed in a LineM and a measure.
But what is the use of an offset? Will the output point be shifted?

Similar to the offset in ST_LocateBetween - Will the output line substring be shifted? The documentation of ST_LocateBetween does not touch the parameter offset.

I know that I would just have to try it but I am still missing these skills, nevertheless I do need the information.


texts: Using extra parameters on ArcGIS print templates service

texts:
I am using ArcGIS Server 10.5. And I am creating my own print template MXD files on ArcMap to print maps on web page. Then I am publishing them on ArcGIS Server. But I want to print extra information on my print template.

First page of print template will include the map view, legend and so on.
Second page will include a table information that I post using service.

If this is possible, how can I do it?


texts: How to return multiple results from the Google Maps geocoding API?

texts:
The Google Maps API Geocoding Developers Guide page shows an example of how to geocode an address.

How can I return multiple address suggestions, so that I can pick the correct response?

For example, a search for "Richmond" should return multiple Richmonds across the world. However, this request only returns a single response (when I use a valid API key):

https://maps.googleapis.com/maps/api/geocode/json?address=Richmond&key=YOUR_API_KEY

The example given at https://developers.google.com/maps/documentation/geocoding/intro#RegionCodes even suggests that multiple results should be returned for address=Toledo but I'm still only seeing one result for Toledo.

Is there a parameter to specify the number of results to return?


texts: AGOL "Update Data" from Excel import numeric columns with some NULL rows?

texts:
I need to append some new features into an ArcGIS Online hosted feature layer.  The data has been supplied to me in an Excel file (including X/Y columns).  I was intending to import it as-is using the hosted feature layer's "Update Data" button.
However, when I get to the field-matching phase of the process, there are some numeric fields in the hosted layer for which it will not let me select the corresponding column from the spreadsheet.  It appears as though this is because that column is blank for some records.
If I change the blank values to '0' (zero), then it will let me select those columns.  So I guess if not every value in the column is a number, then it assumes it must be a text field (and I can select it to import into text fields).
However, these columns need to be imported into number fields, and these blank values MUST NOT be imported as zeros, as that would simple be incorrect and result in wrong analysis.  They need to be imported as NULLs.  I've tried changing the values to null, <NULL> or #N/A but none of these worked - only 0 worked.
The hosted layer field is configured to allow null values.
Is there any way to import directly from an Excel spreadsheet into an AGOL hosted feature layer where NULLs can be imported for a numeric attribute?
(Otherwise, my work around may have to be to convert it to a Shapefile or fGDB first, massage the data there, and import that.  I'm really hoping to avoid this!)


texts:
Can anyone suggest a quick (and preferably open source) way to create an image from a shapefile.  For example, I have a shapefile containing parcel boundaries, I want to create a png file depicting these parcels.  The symbology required for this image will be bare-bones (solid outline and fill).

I know that shp2img from MapServer would do the trick.  I am just curious if there is anything else available.

Any references to pre-compiled tools or api's/sdks welcome.


texts: How to publish shapefiles with GeoServer?

texts:
I have multiple shapefiles of my country, its states and cities. I want them to display in map i.e. First layer is country, then states-layer and then cities layer.

I only know that to do this, first of all these files have to be uploaded in GeoServer, installed on Server, then somehow GeoServer is requested to give desired shapefiles.

Please help me how it can be done i.e. upload files in GeoServer, How to request it and add shapefile into web page.


texts: QGIS3.8 Hotspot Analysis Plugin not installing

texts:
I'm hoping someone on here might be able to help me. I'm fairly new to working with QGIS I'm completely lost. For my PhD project I have to use QGIS to produce georeferenced heatmaps based on informants' hand-drawn maps. An academic from another institution has provided me with step-by-step instructions of what to install and the steps to follow, but I've hit a big problem with installing the Hotspot Analysis plugin in QGIS.

I am working on a Macbook. So far I have installed Python, as well as pip, numpy, scify, and pysal (as instructed). However, when I try to install the Hotspot Analysis plugin I keep getting this same error message:

  Couldn't load plugin 'HotspotAnalysis' due to an error when calling
  its classFactory() method ModuleNotFoundError: No module named
  'pysal.esda'

Does anyone have any idea how to rectify this? I've tried uninstalling and reinstalling everything. I've also tried installing the plugin from a zip file, but I still get the same message.


texts: GDAL with GEOS support, pandas and VTK using python on Windows?

texts:
So far I have been using OSGeo4W shell (or OSGEO install) as really convenient way to use python together with additional packages, namely gdal, shapely, numpy, scipy and matplotlib.

I was happy till I came to the point where I need to start use also pandas and VTK packages.

I came across Anaconda, which allows me to easily install everything including the needed pandas and VTK packages. But when I run my scripts which were working well with OSGEO, with Anaconda it comes with:

ERROR 6: GEOS support not enabled.

I have installed GDAL from https://binstar.org/osgeo/gdal where it is written: "Built with all command line tools and with GEOS and HDF5 support". But obviously something is wrong. I also tried install additionaly first geos and after that gdal from https://conda.binstar.org/osgeo but ended up with the same error.

I found a post talking about this issue linking to special package build gdal_geos on binstar. But there is verion only for linux-64.

So ideally both ways or at least one of these: Can I easily add pandas and VTK support to OSGeo4W shell? And/Or how to make work GDAL with GEOS support on Anaconda install?

Update: I have installed Fiona package to Anaconda and replaced OGR part (load (and buffer with Shapely) objects from shapefile) by it and the code works. If I switch back to OGR, error comes up again.


texts: Changing Legend Symbology to include a 'Sum' field using QGIS 3.16.3 Python Console

texts:
I currently have a vector layer that contains intensity polygons and are colour co-ordinated depending on the level of intensity. Within this 'Intensity' file, there are two columns Intensity and Affected_P. The file currents shows the intensity in the legend for the symbology but I would like to include the sum statistic created - Affected_P too the legend as well. So, for example, the legend with the intensity then the Affected_P field in brackets - 2.4 (60,000). Ideally, done within the Python Console.
The fields/columns look like this. Currently, there are multiple Affected_P values for one intensity. I am not too sure if it would add all the values automatically when adding to the specific intensity legend you want or that you would have to add in a piece of code to this.

The legend looks like this currently but I would like the Affected_P column shown next to the intensity number for each colour

Invalid syntax error when I try to add the code to python, under where I have changed the symbology.
#Change symbology
default_style = QgsStyle().defaultStyle()
color_ramp = default_style.colorRamp('Spectral') #Spectral color ramp
color_ramp.invert()

layer = QgsProject().instance().mapLayersByName('Intensity')[0]
field_index = layer.fields().lookupField('Intensity')
unique_values = list(layer.uniqueValues(field_index))
categories = []
for value in sorted(unique_values):
    symbol = QgsSymbol.defaultSymbol(layer.geometryType())
    category = QgsRendererCategory(value, symbol, str(value))
    categories.append(category)
renderer = QgsCategorizedSymbolRenderer( "intensity" || ' - ' || sum("affected_p","intensity"), categories)
renderer.updateColorRamp(color_ramp)


texts: Custom marker/icon based on attribute data in leaflet/geojson

texts:
I would like to add different icons to my leaflet map based on attribute data. My layer is a geojson layer, it is coming from geoserver with jsnop. I can see my data with the default marker, but when I try to change the code to add custom style I get this error:

  Cannot read property 'labelAnchor' of undefined

It would be great if I can get some help.

var WFSLayer = null;
var ajax = $.ajax({
    url : URL,
    dataType : 'jsonp',
    jsonpCallback : 'getJson',
    success : function (response) {
       WFSLayer = L.geoJson(response, {
                pointToLayer: function(feature, latlng) {
                    var smallIcon = L.Icon.extend({
                        options: {
                            iconSize: [27, 27],
                             iconAnchor: [13, 27],
                             popupAnchor:  [1, -24],
                            iconUrl: 'leaflet/icons/' + feature.properties.pcp + '.png'
                        }
                    });
                    var myIcon = new smallIcon();
                    return L.marker(latlng, {icon: smallIcon});
                },
            onEachFeature: function (feature, layer) {
                popupOptions = {maxWidth: 600};
                layer.bindLabel('
'+feature.properties.musno+'
');
                //sidebar.setContent('
'+feature.properties.musno+'

'+'
'+feature.properties.exchange_name+'

'+feature.properties.pcp, popupOptions);
                layer.bindPopup('
'+feature.properties.musno+'

'+'
'+feature.properties.exchange_name+'

'+feature.properties.pcp
                    ,popupOptions);
            }
        }).addTo(markers.addTo(map));
    }
});


texts:
I want to solve this problem :

given a departure point
given an arrival point
given an amount of time

In which areas can I go while coming from departure to destination in the given amoint of time?

For instance, the faster route take me to the arrival in 2 hours, but I have 3 hours. Which areas can I go to with this extra 1 hour?

All I can find is isochrone or isodistance around a point.


texts: Rasterio reprojection doesn't seem to work

texts:
I would like to reproject a raster using the coordinate system of a shapefile as my destination coordinate system. I am using one of the examples in the rasterio docs as inspiration:.

import numpy as np
import rasterio
from rasterio.warp import calculate_default_transform, reproject, Resampling
import geopandas

src_file = 'output/out.tif'
dst_file = 'output/out_reproj.tif'
dst_crs = geopandas.read_file('confidence.shp').crs['init']

with rasterio.open(src_file) as src:
    transform, width, height = calculate_default_transform(
        src.crs, dst_crs, src.width, src.height, *src.bounds)
    kwargs = src.meta.copy()
    kwargs.update({
        'crs': dst_crs,
        'transform': transform,
        'width': width,
        'height': height
    })

    with rasterio.open(dst_file, 'w', **kwargs) as dst:
        for i in range(1, src.count + 1):
            reproject(
                source=rasterio.band(src, i),
                destination=rasterio.band(dst, i),
                src_transform=src.transform,
                src_crs=src.crs,
                dst_transform=transform,
                dst_crs=dst_crs,
                resampling=Resampling.nearest)

The shapefile is in EPSG::3035 and the input raster is in EPSG::4326. The output raster is still in EPSG::4326.

Am I doing anything wrong in the script?


texts: R: Writing individual variables from a RasterStack in loop

texts:
I have a RasterStack with 3 layers each one of them containing 19 variables. Using writeRaster I need to save specific variables into my directory while naming the saved files with the actual variable name and it's respective layer id.

If one were to take the variables one by one in each layer...

writeRaster(myRasterStack$var1.1,  "var1.1.grd", format = "raster")

var1 indicating the actual name of the variable in the raster object and .1 it's respective layer.

I'm working on the following loop.

variables <- as.factor(c("var1.1", "var8.1", "var1.2", "var8.2","var1.3", "var8.3"))
for (i in 1:length(variables))
{
  writeRaster(myRasterStack$[i], filename = paste0("variables", i,"/data/.grd"), format = "raster")
}

What am I missing in the loop?


texts: Import Vector data to QGIS from pdf and apply georeference

texts:
I am Using QGIS 2.14. I have a PDF created by Mapinfo, it contains multiple layers of vector data that I can open as separate vector layers in QGIS.

However no map projection data was saved to the PDF. One of the vector layers is a UTM grid. Is it possible to manipulate the vector layers using the information from the grid layer and save as correctly located shp files?


texts:
I've been having issues with displaying basemaps using the QuickMaps Services plugin and I've tried the connecting to Tiles servers through the XYZ Tiles channel to display basemaps. Using both methods I just get a blank work space with no basemap.
I have reason to believe that there is something wrong with the server connection as I get this error when I try to Check Version Number:

Also, I can't seem to access the Plug-in library and the repository is indicated as "unavailable":
Has anyone encountered what I believe is a server connection issue and has anyone been able to fix it?
SOLVED: It turns out I had Use Proxy for Web Access checked before which resulted in problems accessing the internet


texts: ESRI Geocoder placeholder text not displaying

texts:
Following the ESRI API, I added the placeholder parameter to the geocoder widget to add some informative text would instruct the user what the search box was for.  For some reason it is not displaying anything on my page - any clues?

Geocoder Widget Code:

geocoder = new Geocoder({
                map: map,
                placeholder: "Locate by address"
            }, "search");
            geocoder.startup();

and then the simple "search div" is loaded in a content page (not main map window):


texts: Getting popup with layer attributes to display after using address search/geocoder widget?

texts:
Using Web AppBuilder and ArcGIS Online, I want to let my user(s) enter an address on the "Geocoder Widget" and have an automatic popup that displays. I know how to configure the content of the popup. I know I can change the Geocoder to search for specific features/attributes.

Example: I have three school boundaries (High School districts, Middle School Districts, and Elementary School Districts) that overlap each other. I want to create an app that lets someone Enter their address and it spits out their neighborhood HS, MS, and ES (I will be adding driving directions and descriptions of the schools in the pop up).

ArcGIS for Javascript : Locate Widget Activate Popup for Layers X > Z? asks the same thing I believe but was never answered or commented on.


texts:
The "Vectorial" menu by qgis is empty. I need the data management tools. How can I make it appear?
I tried reinstalling Qgis 3.4.10, and nothing changed.


texts: Generating regularly spaced sample points in qgis based on ratio of area

texts:
I am generating sample locations in irregular shaped polygons in qgis. In my situation, there are often many polygons in one shapefile all with unique shapes that need a certain number of sample units on regular spacing.

Currently, I am only able to work with one polygon at a time using the following process:

First, I use the qgis generate a regular grid tool. As I have found, this will generate a grid over the extent of the polygon. The points are not constrained to be within the shape. I have to guess a number and hope its the correct amount.

Second, I clip the new sample unit layer and the polygon layer so that only points that fall in the polygon are in the new layer.

Third, I count the number of sample units and see if it is an appropriate number. If it is not, I start over at step one with a different number of points for the regular grid.

My idea is that it would be great to automate this process using Python. I am a novice Python programmer.

The logic that I have developed in my head is as follows: (I am going to use the terms Polygon, Sample, and SampleClipped)

I would like to use some sort of Do Until procedure for each polygon. It would first create a grid starting at 1 sample unit in the grid. On each iteration, the loop would include clipping the Sample and check SampleClipped to see if it meets a specified ratio (say 1 sample unit per 4 acres, assuming that acres is in the Polygon attributes). If 1 sample unit in the grid does not meet this ratio, it would try 2 and so on until it meets the ratio. Since often, there are more than 1 polygons in a shapefile, it would need to loop through all polygons in the shapefile.

Another option would be the same procedure above except there would be a field in the attribute table that is Minimum Number Sample Units. The code would loop until the first instance in which it is at least that number of sample units. It can be over that number but not under.


texts:
I'm having a polygon-file of all municipalities in Bavaria, with data of population density. Now I want to get a value for each polygon (= each municipality) that shows the average of population density of all neighboring polygons plus the source polygon.

I am working in QGIS. My knowledge is quite basic. Can anyone help?
As I want to get the value for all municipalities of Bavaria, I can't do this operation by hand.


texts: Nearest point on a line ( Spherical / Mercator projection )

texts:
I have a line (Ax,Ay - Bx,By) over a mercator projection (google maps) and a random point (Cx,Cy) nearest to that line, i would to know the closest point (transparent blue on the image) over that line to point (blue in the image)

EDIT: to clarify that this is in a Mercator projection (spherical projection)


texts: So I geo-referenced a few scanned maps successfully, however when it orients the scanned map correctly, a black rectangle background appears behind the map.  This doesnt allow for me to "puzzle" together different maps because the black background lays on top of other scanned map portions.  I was wondering if there's a way to crop out that background while obviously keeping the map part intact and correctly oriented/geo-referenced.


texts:
ArcGIS 10.1 - attempting to add a derived raster layer to the ToC.  This works in 10.2.  In 10.1 I get a screen flash and nothing added.

rasterized = r'in_memory\rasterized'
arcpy.PolygonToRaster_conversion(poly_shp, field, rasterized, None, cellsize)
res = arcpy.MakeRasterLayer_management(rasterized, 'tmplyr')
rasterlyr = res.getOutput(0)

mxd = arcpy.mapping.MapDocument("Current")
df = arcpy.mapping.ListDataFrames(mxd, '*')[0]
arcpy.mapping.AddLayer(df, rasterlyr, "TOP")

Where:
poly_shp = unicode vector name
field = unicode field name from the poly_shp
rasterized = str
cellsize = float

Adding:

arcpy.RefreshTOC()
arcpy.RefreshActiveView()

has no effect.

Edit:

I know that the in_memory workspace is functioning properly as the following try is hit on subsequent runs.

try:
        arcpy.Delete_management(rasterized)
        arcpy.AddMessage('Overwriting in_memory raster')
except:
        pass


texts: PyQgis open vector layer attribute table

texts:
Is there a Qgs dialog or widget that I can use in a standalone pyqgis application to open the attribute table of a QgsVectorLayer? this seems to be such a simple thing, there has to be an existing widget or dialog...

Update:
After some extensive search, there seem to be two ways to do this in pyqgis, neither worked well, I am just going to summarize them here:

Approach 1. use the QgsDualView class as below:

self.editor = QgsDualView(self)
self.editor.init(self.layer, self.map_control.canvas)
self.editor.setView(QgsDualView.AttributeEditor)

This will display attributes as below:

Approach 2. use the QgsAttributeTableView

lyr_cache = QgsVectorLayerCache(self.layer, 120, self)
tbl_model = QgsAttributeTableModel(lyr_cache, self)
tbl_model.loadAttributes()
tbl_model.loadLayer()
tbl_filter_model = QgsAttributeTableFilterModel(self.map_control.canvas, tbl_model, self)
self.editor = QgsAttributeTableView(self)
self.editor.setModel(tbl_filter_model)

This will display as below:

So, as you can see, both look very much half-baked. I did accept @Matthias Kuhn's reply as answer as it did answer my original question, but in light of the fact that it didn't actually work, I continue to hope that someone can give suggestions that could make either approach work better. I will certainly up-vote your answers and give confirmation, this will benefit the larger community.


texts: How can I avoid changes to field types in geopandas?

texts:
I am using geopandas in the following routine to create the union of 2 (polygon) shapefiles:

import geopandas as gpd
def union(fn_A, fn_B, output):
    shpA =gpd.read_file(fn_A)
    shpB =gpd.read_file(fn_B)
    union = gpd.overlay(shpA, shpB, how='union')
    union.to_file(driver = 'ESRI Shapefile', filename= output)

Shapefile A has an integer attribute "a" and shapefile B has an integer attribute "b". The resulting shapefile has attributes a and b but of them are doubles. I have read that geopandas infers a schema when saving to files, but do not know how to even get the schema of the inputs in the first place using geopandas. It appears that this only happens if the result of the union contains empty fields (i.e. areas where A and B don't overlap). If A is completely covering B, field a remains integer.

How can I avoid this conversion?


texts: Problem with access to ArcGIS Server Manager

texts:
ArcGIS For Server 10.0

Firefox can not establish a connection to the server. (Error in all browsers) ArcGIS Server Manager problem

In ArcCatalog i can connect but can't start

Configuration MapServer can not be started.

Also in computer management missing Local Users and Groups

What is my fault? Thanks


texts: QGIS equivalent of ArcMap spatial join - share a line segment with

texts:
I have two polylines shapefiles that I wish to join together spatially.
One is a road network and the other consists of bus routes which I have calculated the frequency over a certain period of time. I then wish to change the size of road segments based on this frequency.

I have found a successful way to do it in ArcMap, using the "SHARE_A_LINE_SEGMENT_WITH" option in the Spatial Join tool. Afterwards, I dissolve the joined shapefile based on the road_id column and sum the frequency column, also in ArcMap.

As shown in this great visual example, share a line segment with is quite different from intersect, contains, within or other methods offered in QGIS.

Is there some way to perform this type of spatial join in QGIS, whether with a tool, plugin, or spatial SQL query?

My python knowledge is really poor.

How to find share segments between two shapes (line) using QGIS? seems to be the same question as mine, however I can't comment to ask if there is any update.

Also, I found Understanding Join Attributes by Location in QGIS? for QGIS, where no operation seems to be what I want exactly.

This is the road network:

With bus routes on top:

When I do a simple intersect, i'm left with a bunch of routes (in red) that shouldn't be there (circled in blue). I' guessing because of the segment ends.


texts: Conditional rasters value replacement in

texts:
In Google Earth Engine I would like to replace all the negative pixel values of each image of an ImageCollection with 0.

My approach is the following:

var collection = ee.ImageCollection('NOAA/VIIRS/DNB/MONTHLY_V1/VCMSLCFG').filterDate('2014-01-01', '2015-12-01').select('avg_rad');

var replacement = ee.Image(0);

var conditional = function(image) {
  image.where(image.lt(0), replacement);
};

var output = collection.map(conditional);

Howeve, I get the following error ressage

User-defined methods must return a value.

How can I solve?


texts: How draw a route in Google Maps as i goes on clicking on Google map?

texts:
I wanted to get route in Google map whenever I click on the road map of Google. What I mean by click is that when I can click one place then other and so on, but route to each click must be generated from previous click. Route must be shortest route from one point to other. I must be able to drag that click point when required and then new route must be generated for last two click.

Can any one suggest how i should do this?

Note: No marker required.

Hope could have express my problem effectively.


texts: clear map zoom and geometries of feature layer

texts:
i have selected feature(feature layer) from combobox and it zoom to the feature . now i want to clear combobox slection and map . and zoom map to its default zoom.

//combobox selection clear
    dijit.byId("A1").reset();
    dijit.byId("A2").reset();
    dijit.byId("A3").reset();
    dijit.byId("A4").reset();
    dijit.byId("A5").reset();

    //layer selection clear
     document.getElementById('A1_layer').clearSelection();
      document.getElementById('A2_layerC').clearSelection();
      document.getElementById('A3_layerC').clearSelection();
      document.getElementById('A4_layerC').clearSelection();
      document.getElementById('A5_layerC').clearSelection();

app = {
    zoomRow: function(id, which){

      var query = new Query();
      //var thePoly, theExtent;
      if(which == "Land"){
        query.where = "Name='" + (id).toString() + "'";
        console.info(query.where);
        query.returnGeometry = true;
        A1_layer.selectFeatures(query, FeatureLayer.SELECTION_NEW, function (features) {
          thePoly = features[0].geometry;
          theExtent = thePoly.getExtent().expand(2); //Zoom out slightly from the polygon's extent
          map.setExtent(theExtent);
        });
        esriRequest({
          url: "http://localhost:6080/arcgis/rest/services/........",
          content:{
            f:'json'
          },
          handleAs:'json',
          callbackParamName:'callback',
          timeout:15000
        }).then(lang.hitch(this,function(response){
          var store2 = new Memory({data:[]});
          dijit.byId("A2").set('store',store2);
          var data = array.map(response.features,lang.hitch(this,function(feat, index){
            var name = feat.attributes.nam;
            var dataItem = {
              id:index,
              name:name
            };
            return dataItem;
          }));
          store2 = new Memory({data:data});
          dijit.byId("A2").set('store',store2);
          document.getElementById('A2').value = "Select Room";
        }));
      }
<input id="A1" data-dojo-type="dijit/form/ComboBox" value="Select landing" onchange="app.zoomRow(document.getElementById('A1').value, 'Land');" data-dojo-props="maxHeight: 200" style="overflow:auto; width:200px; background-color: #E7FCCA "/ >

  <input id="A2" data-dojo-type="dijit/form/ComboBox" value="Select room onchange="app.zoomRow(document.getElementById('A2').value, 'Room');"style="overflow:auto; width:200px ;background-color: #E7FCCA" />

but i got error "Cannot read property 'geometry' of undefined" at thePoly = features[0].geometry;


texts:
I'm trying to get proj strings from cartopy crs', such as cartopy.crs.PlateCarree().  I can manually infer the proj string from the init() function in the source code but I need to query it, so I can apply it to a tiff file with gdal.
I'm trying to avoid adding modules outside from scipy, numpy, gdal, cartopy, OWSLib, MatPlotLib, and Pillow.
Is there a solution?


texts: What is the correct projection to use with tm_wgs84_sans_antarctica.shp?

texts:
I am using tm_wgs84_sans_antarctica.shp with mapnik to generate tiles. I am using the tiles with Google Maps Javascript API V3. But the projection is way off. Here is roughly what my code looks like:

import mapnik
m = mapnik.Map(
      tile.size()[0],
      tile.size()[0],
      )
s = mapnik.Style()
r = mapnik.Rule()

polygon_symbolizer = mapnik.PolygonSymbolizer(mapnik.Color('#f2eff9'))
r.symbols.append(polygon_symbolizer)
s.rules.append(r)

m.append_style('My Style',s)

ds = mapnik.Shapefile(file='tm_wgs84_sans_antarctica.shp')
layer = mapnik.Layer('world')
layer.datasource = ds
layer.styles.append('My Style')
m.layers.append(layer)

bbox = tile.bounds()
bbox = mapnik.Box2d(bbox[0], bbox[1], bbox[2], bbox[3])
m.zoom_to_box(bbox)
im = mapnik.Image( *tile.size() )
mapnik.render(m, im)

As you can see I am not setting any srs either in map or in layer. I tried different combination of projection strings without really understanding what they mean, but either got no map or wrongly projected map. Please let me know if any other information is required.

Thanks.


texts: Using one massive map service vs multiple smaller services in ArcGIS Server?

texts:
I have a Silverlight application that consumes around 30 map service layers at different types.  The service layers are categorized based on type and only applicable to certain screens of the application.

For organization's sake, I split each of the layer categories into separate map services, but I'm considering whether or not that makes sense.

Should I just lump all of the map layers into a massive map service or keep them separate?

I'm curious about server load vs application performance.


texts: Getting Negative values in Regression based Biomass model!

texts:
I am using Multi Linear Regression (stepwise method) to estimate above ground biomass. Different spectral indices and other variables mostly derived from Landsat 8 OLI are used to check potential relationship with the biomass samples. (The biomass samples range between 7 to 37 tons). The adjusted R Square hovers around 0.6 which is respectable.

The resultant equation contains coefficients which are statistically significant. However, the input equation produces raster with values ranged between -24 to 57.

I am not sure what might be the source of these negative values in the raster and how should I deal with these negative values?


texts:
I know it is curious, but I want to run GeoServer 2.3.0 on Tomcat 5.5.20. When I deploy the war file it deploys successful but I can't start the GeoServer.

Message:

FAIL - Application at context path /geoserver_2.3.0 could not be started

Unfortunately I must use Tomcat 5.5 and until GeoServer 2.1.4 I can deploy and start fine.

You know how can I solve this or you have any experience.


texts: Unable to create Polygon using Geodjango admin interface

texts:
I am new to GeoDjango. I have created Admin interface for my Model. It loads correctly. I am also able to draw the Polygon on map(At client side) but when I save It tells me Invalid input and throws some error in console.

Model

from django.contrib.gis import models

Poly(models.Model):
    poly_gon = models.Polygone()
    objects = models.GeoManager()

Admin interface

from django.contrib.gis import admin
from myapp.models import Poly

admin.site.register(Poly)

Error at Client side(On the polygon field in red colour)

  "Invalid geometry value."

Error at server side(in terminal, after starting server with runserver command)

Error in server side: GEOS_ERROR: ParseException: Expected 'Z', 'M', 'ZM', 'EMPTY' or '(' but encountered : ')'

GEOS_ERROR: ParseException: Expected 'Z', 'M', 'ZM', 'EMPTY' or '(' but encountered : ')'

Error creating geometry from value 'SRID=4326;POLYGON()' (Error encountered checking Geometry returned from GEOS C function "GEOSWKTReader_read_r".)

I was able to enter data to my Poly from shell by following code.

from django.contrib.gis.geos import GEOSGeometry
from myapp.models import Poly
p = Poly(poly_gon=GEOSGeometry("""POLYGON((-87.6705551147461 41.89135645852043, -87.64171600341797 41.89288988217029, -87.63690948486328 41.880110226947934, -87.66368865966797  41.87806524488436, -87.68016815185547 41.887267148816726, -87.6705551147461 41.89135645852043))"""))
p.save()


texts: How to create a function for finding the largest of 2 polygons by area in PostGIS?

texts:
I tried my level in creating the function but it shows error while running.

I have enclosed the code here :

create function largest_poly (geom1, geom2)
returns geometry as $largest_poly$
begin
from plpygis import Geometry
ST_Area(geom1 ST_Polygon),
ST_Area(geom2 ST_Polygon)
 poly largest_poly=Null
for
 (poly ST_Area: geom){
 {If (largest_poly==Null)|| ST_Area.geom() > largest_poly.ST_Area())
    largest_poly=ST_Area
 }
 return largest_poly;
}
end;
$largest_poly$ language plpygis;


texts:
I've been uploading a huge amount of contour information to a PostGIS / PostGres table. Uploaded using QGIS' db manager to a separate table a chunk at a time, and then using a quick bit of SQL to append those chunks onto the main contour table.

pgAdmin is estimating it to have 492,408 rows containing five columns, which are:

Altitude (double) e.g. 129.5
Class (char) e.g. Basic
geometry (geometry, a multiline string)
OID (integer)
pid (integer and primary key)

Now, if I try and add it into QGIS or Manifold (using Add PostGIS layers) it will basically take ages if it can add at all.

I've added a spatial index to it, as well as an index on the pid. Indexes that I've added are:

CREATE INDEX "Contours_50cm_Idx"
ON "Contours_50cm"
USING gist
(geometry);

and

CREATE UNIQUE INDEX "PID_INDEX"
ON "Contours_50cm"
USING btree
(pid);

Does anyone have any idea of anything else that could help optimise this?

Google has revealed people sometimes use partitioning on tables. Users won't have to see all of this data at a time - it covers quite a large area geographically.


texts: Why does one polygon disappear when merging in edit mode?

texts:
I want to merge some polygons into one using the merge tool. But when doing so, one of the polygons just disappears (see screenshot 2). Why does this happen?
When I use the simplify tool before, it does not happen.

My system: QGIS 2.6.0; Linux Mint 17


texts: QGIS topology checker - check if point features are covered by line features

texts:
I want to use the topology checker to check if point features are covered by line features. For example are the bus stops covered by the bus lines.Therefore I use the topology checker with the following rule: "points" must be covered by "lines".

"Haltestellen" = bus stops and
"Buslinien" = bus lines

But I get the error for all points that "point is not covered by segment".

I even get the error for points that I edited on the line with snapping.

Here are my questions:

Is it possible to check if point features are covered by line features? With topology checker? In general?
What does the rule: "points" must be covered by "lines" check? Under what cicumstances is a point feature covered by a line feature?

I tried this in QGIS Versions 2.18; 3.4 and 3.10

My OS is Windows 10.

I tried with shapefiles and geopackage layers

The CBS of point and line layers are the same. Here I tried the coordinate systems with epgs: 31468 (DHDN / 3-degree Gauss-Kruger zone 4) and 25832 (ETRS89 / UTM zone 32N).

Can anyone explain or point me to the right direction?


texts: Exporting Georeferenced PDFs using ArcGIS Desktop 9.3?

texts:
I am using ArcGIS Desktop 9.3 and I want to export maps in georeferenced PDFs.

Earlier I knew that there is an advanced option in file>export map>options>advanced where I was able able export map georeferencing information.

However, when I try to export following the above instructions I cannot find (export map georeference information) in advanced option.

What is the problem?


texts: Help on creating a presentation for Highschool students

texts:
I was commissioned to do a presentation on GIS (well it'll mostly focus on transportation GIS) for about 50 high school students.

I'm looking for funny graphs / images that I can include in this presentation.
What are some activities (split them up in groups, or interactive questions/discussions) I can do to engage them and get them interested?
Any suggestions? comments, or concerns?
Mostly these kids will be forced to be there with a lot more interest in the person beside them, or what's for supper, or an xbox than GIS. How can I get their attention and keep it?

The presentation will run about 45 to 60 minutes + 30 minutes for activity / fun time. We'll be in a hall with 5 to 6 round tables

I know this should be CW, but the feature has been disabled in this site


texts:
My question is how to test each element in one list with each element in another list? I have lot of lines segments and using lengthCalculator to find all segments less than some specific length. Further with help of TopologyBuilder I get 'to_node' attribute related to LINE output port and 'node_number' related to NODE output port. I have to detect all dangles. So I am trying to test if there is a "node_number =! to_node" If so, then I know at this segment is dangle. Agree? for this I used transform Tester again, but this is not working. My next idea is to create list of all nodes and another list of all segments < specific_length, and then check (test) each node element in one list with to_node element in teh second list. Will this work? Not good idea?
This question is linked to: Find dangling lines .
Tnx


texts: QGIS: Expression WHEN property1 = 'FOO' THEN property2 = 'BAR'

texts:
Our Foundation is adding additional properties to points.

If property1 = FOO , then add property2 = BAR

    WHEN
        property1 = 'FOO' THEN property2 = 'BAR'
    END

When I run this the property1 column is deleted. Can you help me with the syntax, here, please?

Assumptions/Knowns:
This is a spatialite layer. property1 and property2 columns both exist. I am using QGIS 2.4. I opened the attributes tables to use the Expression dialog. The property1 column is selected when the expression is executed. property1 and property2 are columns, and I am aware that " is sometimes used to enclose columns. I do need to store these values for use in other programs where they cannot be so easily derived from the data.

I saw several examples that changed the current value of property1: WHEN property1 = A THEN 1 but this only changes the value in the query column, and I am looking to insert another value property2 based on the value of property1.


texts: ModelBuilder unable to overwrite existing features with same names

texts:
I have created a model that is meant to be run regularly to update an existing data set with new data. Ideally, the model will be able to simply overwrite the old data with the new data which would have the same file name. However, when I try to run the model, I get an error that it cannot run because files with those names already exist. I already have "Overwrite the outputs of geoprocessing operations" checked on, so I'm not sure what else can be preventing it from running.

Are there any other parameters that need to be set to enable automatic overwriting?

Everyone on the internet seems to be having the opposite problem that I am!

Start Time: Tue Feb 12 16:27:35 2019
ERROR 000210: Cannot create output O:\GIS\ArcGIS_Pro\Data_Analysis_and_Planning\gswiBasin
The table already exists.
The table was not found. [gswiBasin]
Failed to execute (Copy Features: Copy files into working database).
Failed at Tue Feb 12 16:27:36 2019 (Elapsed Time: 0.86 seconds)


texts: I have 15 years of monthly air pollution values in multiple monthly tables I interpolated and then summarized by table for each census division in ArcGIS.  I want to summarize the mean monthly values into yearly then map out the yearly values in ArcGIS.  (i.e. take 12 tables and calculate mean yearly values).  Each table has a pollution value and unique census id.  Any ideas the best way to do this?


texts: Why gdalwarp cannot make transparent black nodata?

texts:
I have a MODIS ariel image geotiff and I want to remove black nodata pixel (this is Aerosol Of Thickness) 1 band image.

I've done like this C:>convert -transparent black MOD04.A2003226.0325.006.2013074065606.tif output1.tif

and the output1.tif is large x20 times input MODIS tif file. But it is blacker than input, no transparent black color.

Some one can help me? Please!

This is input sample (i need remove all black nodata point)
Input image

And this is output converted (it's worse)

output image


texts: Using field calculator to find fields in row that have same value as reference field then return string gives ERROR 000539

texts:
I  am trying to create a field "conflict", and the condition is when other fields in that row have the same maximum number the function returns "Yes". I am having trouble making this work in the field calculator in ArcGIS Pro.This is what my attribute table looks like and the error I get (see pictures below).
I eventually want to see which fields(Koala, GHFF, Education, Cropping, Conservation, Cattle, and Carbon) have the same value as the "Max value" field, then use the conflict field to indicate if there is more than 1 field with the same value as the maximum value.
I'm fairly new to coding.
def conflict (Ed,Con,max):
    if max in (Con and Ed):
        return ("Yes")
    else:
        return ("No")
conflict(!Education!,!Conservation!,!MaxValue!)


texts: How to convert Geostatistical analyst (GA) file to Raster FAST?

texts:
I have been trying to convert a number of Geostatistical analyst files that i were created using Geostatistical wizard in ArcGIS 10. I want to save these files as Raster and use the
File

Am i missing something, is there a way to do it faster ?


texts: Intersection of convex polygons on a spherical surface

texts:
There are a number of algorithms available for the intersection of two convex polygons, but I am interested in an algorithm to find the intersection of two convex polygons on the surface of a sphere, where the edges of the polygons are great circles. The polygons might be arbitrarily large (i.e. no approximation for small areas).

Is there a description of such an algorithm available? Better, is there any Python or C code out there that would do this, and has a liberal open source license or is public domain?


texts: Masking dynamic leader lines for pie charts in QGIS

texts:
I am trying to duplicate
the example in this post Moving Pie Chart in QGIS and creating leader line with some modifications. I'm trying to combine this with transparency and I'd like to mask off the leader lines based on the pie chart size (so the line appears to connect to the circumference and not the center of each pie).

I think I need a masking solution because the pie chart size varies dynamically and is not always predictable (ie, there is a minimum size bound below which pies are presented at a particular size, but not a maximum size). Any ideas? QGIS 3.16


texts: Polygon st_union in sf drops some polygons from output

texts:
I have written the following code that downloads state boundaries and population data from the US Census and then recursively merges the least populated states with their least populated neighbors. The idea is to come up with 20 districts of more equal population than the 50 states. As implemented below, some polygons disappear. Anyone know how to fix that problem?
library(sf)
library(tidyverse)
library(tidycensus)

key <- "your key goes here (see tidycensus docs)"

census_api_key(key)

states <- get_estimates("state", "population", geometry = TRUE) %>%
  st_transform(2163) %>%
  spread(variable, value) %>%
  st_simplify(TRUE, 1e3) %>%
  st_buffer(1e3) %>%
  rename(id = GEOID, name = NAME, population = POP, density = DENSITY)

merge_states <- function(x) {
  i <- which.min(x[["population"]])
  j <- st_intersects(x[i, ], x[-i, ])[[1]]
  k <- j[which.min(x[["population"]][j])]
  if (length(k) < 1) {
    l <- seq.int(1, nrow(x))[-i]
    m <- which.min(st_distance(x[i, ], x[l, ]))
    k <- l[m]
  }
  x$geometry[k] <- st_union(x$geometry[k], x$geometry[i])
  x$population[k] <- x$population[k] + x$population[i]
  x <- x[-i, ]
  return(x)
}

zones <- states

while (nrow(zones) > 20) zones <- merge_states(zones)

#plot(zones[, "id"], key.pos = NULL)

ggplot(zones) + geom_sf(aes(fill = id)) + guides(fill = FALSE)


texts: How to hide labels in intersecting layers using label expressions in QGIS?

texts:
We are using views to build dynamic tables of utilization and enrollment. When mapping, however, the schools that meet both thresholds both display, and I can't figure out how to remove the labels from one layer if the same feature is in the 2nd layer:

I have seen other posts that discuss using the get_feature() expression builder, but I can't see how to get that to work.

Here's an attempt trying to hide the labels from the enrollment layer if the schnum is present in the Current Utilization layer using a label expression:

case
when schnum in (get_feature( 'Current Utilization', 'schnum', '*'))
then ' '
else abbreviation
end

But that doesn't seem to work, and I can't find much documentation to explain if this is even possible.


texts:
Whilst reading the documentation about using the PostGIS function ST_PixelAsCentroids() I was confused somewhat by this (ST_PixelAsPoints(rast, 1)).* syntax and as a newbie to SQL I was wondering if someone could explain this convention and why it's used?

It's obviously used as a method to unbundle complex results but I don't even know what it is called to google it properly to get a formal understanding of it.

I realise this post may be better suited to a SQL forum on Stack but I'll see what GIS Stack conjures up.


texts:
I downloaded some rasters of climatic variables included in TerraClimate (Abatzoglou et al. 2018), from http://www.climatologylab.org/terraclimate.html

I know that each raster is huge, but my computer has been trying to open just one of them and is lasting too long: since 10 hours ago or more. I wish to know if anybody out there was able to open them in R.

I did the following:

library(raster)
library(ncdf4)
library(foreach)
library(doSNOW)

cl <- makeCluster(8, type="SOCK") # for 4 cores machine
registerDoSNOW (cl)

setwd("MY/WORKING/DIRECTORY/")
temp <- nc_open("ppt_1958.nc")

I also tryied, with the same result:

temp <- raster("ppt_1958.nc")
temp <- stack("ppt_1958.nc")

My computer is running under Windows 10, 64 bits, has i7 processor and 32 GB RAM.

  Abatzoglou, J.T., S.Z. Dobrowski, S.A. Parks, K.C. Hegewisch, 2018,
  Terraclimate, a high-resolution global dataset of monthly climate and
  climatic water balance from 1958-2015, Scientific Data,


texts: Field Calculator is not calculating selection from Make Feature Layer in arcpy

texts:
I am trying to use field calculator to calculate 3 different fields in python from a Make Feature Layer variable.  It calculates the first field (MATERIAL_T) but then it doesn't calculate the second two fields (COMP_YR_N and Hierarchy_Material).

cal = arcpy.MakeFeatureLayer_management(servLoc, "ServiceLocation", "MATERIAL_T IS Null", outLoc)

arcpy.CalculateField_management (cal, "MATERIAL_T", "\"Unknown\"", "PYTHON")
arcpy.CalculateField_management (cal, "COMP_YR_N", 1800, "PYTHON")
arcpy.CalculateField_management (cal, "Hierarchy_Material", "\"Unknown\"", "PYTHON")


texts: How many ArcGIS for Server service instances is optimal?

texts:
What is the best number of ArcGIS for Server instances?

When I increase the number of instances will the services run faster?


texts: How to read oracle not spatial data using pyqgis

texts:
I'm working on a QGIS plugin. I'm trying to join a vector layer and a non-spatial table both in Oracle spatial.

I add the vector layer to the map, but I can not add the table without geometry so i can join the layers. I have to pull it out and add in CSV format then I add it as CSV file, do I?


texts: I have table with column name,Latitude and Longitude.I want to update this with polygon name and its centroid latitude longitude.I know there is a function in postgis st_centroid(geometry). But I want to enter x value and Y value in a separate column.How can I do this?


texts: Creating grid from 4 sided polygons with equally divided facing sides using FME

texts:
I have a set of polygons, each one comprises 4 sides. Each facing sides are equally divided into segments (right figure).
Is there anyhow to automatically create a grid as shown in the figure on on the left ?
I am using FME


texts: Reducing File Size without losing quality

texts:
I am using QGIS and currently have Tiff files that when uncompressed are upwards of 1 GB.  When compressed, they are around 250 MB.  I am trying to reduce this file size dramatically whether it be by converting to a different file type (jpeg, png) or doing something else from scratch.  The Tiff files are georeferenced aerials so if I convert to jpeg or png then I will also need a world file to reference my images.  The only way I have been testing this is through Raster > Conversion > Translate where I select my input layer and my output file type.  I am trying to reduce these files to around 100 MB if possible.

If there is a different or better way to reduce these files let me know.  The image is an example of the type of image I am using.


texts: QgsVectorLayer vs QgsVectorDataProvider for addFeatures() function in QGIS

texts:
Has the addFeatures function moved between QgsVectorLayer and QgsVectorDataProvider?

Do both work? Is one deprecated. The apis do not provide that info. Nor does a google search.

Which is preferred and why?

What am I misunderstanding here? What is the situation that is changing which means I may need to change this?


texts: Is there a way to view the list of assets in a shared Earth Engine asset folder?

texts: If my co-worker shares an asset folder with me (e.g. user/coworkername/somefolder) and provides writer permissions, then is there a way that I can view the contents of this folder within the 'Asset Manager' on the Google Earth Engine interface. Otherwise, I have to know the exact name of each asset in order to import it, rather than being able to browse the assets within that folder to identify the one that I need.


texts: Does anyone know about a web client for GIS from the Vaadin framework?

texts:
Reviewing the universe of GIS tools on the Web, I see that several customers have already developed GIS on the Web, so the is no need to start from scratch. Now I have come across the Vaadin framework that appears appropriate for someone like me who does not have much experience with web technologies (html, javascript, css, etc.).

Could you tell me if there is any web client for GIS from the Vaadin framework?


texts: openlayers 3 how to get map scale using zoom level

texts:
is there a way to get the map scale on openlayers 3 the map change zoom level ?

i found a site that converts from map distance to scale

zoomlevel : 20
to
map scale : 1:50000.

https://www.maptools.com/scale_calculator

is there a way in openlayers to get the map scale using the zoom level as map distance


texts: Spatial Joining multiple features (from one layer) to one feature os other layer on QGIS 3.4

texts:
My objective is similar to the one especified on
Joining Multiple Features to One Feature
The thing is that I'm on QGIS, and seems that tools are slightly different, so that confuses me on how to use the answers given to that question.
I'm working with a shapefile of urban lots (this is the shapefile I'm trying to copy the attributes to) and a shapefile with several coverage areas related to cultural patrimony elements in town. These coverage areas are the features overlapping inside the same layer.
Also tried the one-to-many spatial join mentioned by our friend, but the data is already very large (the lots attributes counts almost 2 million lines and 20 columns), so that is not much interesting to my needs, just like isn't much helpful the solutions with many processing steps.


texts: Disable DELETE for PostGIS Role/USER yet allowing them To INSERT new Features?

texts:
I intent to block my USERS from deleting features on PostGIS Database. I disabled DELETE in pgAdmin as below but the USERS are now unable to save feature on the database using QGIS3 as client. On trying to save any changes in QGIS it give  Commit Errors: Could not commit changes to Layer...

 How then do I stop users from DELETING features yet giving them privileges of ADDING NEW features?


texts: I try to update QGIS to 3.14 version with osgeo4w-setup-x86_64 but 3.14 does not appear in list. It only shows 3.12


texts: Performing buffer analysis using Python in QGIS?

texts:
Can we create a buffer of distance 0.004 using python in qgis.

I tried Buffering in pyQGIS?, but it isn't working for python in qgis 2.10.1-Pisa.

The above is the diagram of buffer I am asking. It's located in the vector->geoprocessing tools->Buffer(s)
Any of you could help me out with the equivalent python code for the above Buffer in QGIS.
code:

from qgis.utils import iface
from qgis.analysis import QgsGeometryAnalyzer
QgsApplication.initQgis()
uri = QgsDataSourceURI()
uri.setConnection("localhost","5432","gis_data","postgres","password")
uri.setDataSource("schema_name","table_name","wkt","")
vlayer = QgsVectorLayer(uri.uri(),"LayerName1","postgres")
QgsMapLayerRegistry.instance().addMapLayer(vlayer)
mc = iface.mapCanvas()
layer = mc.currentLayer()
QgsGeometryAnalyzer().buffer(layer, "file:///C:/Users/***/Downloads/aa.shp", 0.0004, False, False, -1)

Everything works fine till the LayerName1 loads and getting errot at Buffer statement the last line of code error:
Feature geometry not imported (OGR error: Pointer 'hGeom' is NULL in 'OGR_G_ImportFromWkb'.
            )


texts:
I'm working with Retrofit to connect my application with my API. Retrofit is using Jackson 2+, and ESRI is using Jackson 1.95+
So, I can't figure out how to solve this.

I have an application using Google Maps. I'm using Retrofit to connect my Android application with my REST service, so I can manage HTTP requests, send and receive information from my server hosted in AWS. I want to add a module, so the user can decide if he wants to use GoogleMaps or his own map imported from ArcGIS Online.

First, I spent hours trying to import these to IDEA 14, I was getting DEX problem (Library Repeated), looking deeply I found that Retrofit was using a newer version of Jackson, and ESRI is using an older one, so, I erase the ESRI Jackson Library and it seems to work, but when I run my app I'm getting a black screen on my MapView.

Any ideas?


texts: Pandas to GeoPandas : Polygon Geometry, WKTReadingError

texts:
I read a PostGIS table with pandas : df = pd.read_sql_table(con=engine, table_name)

The column geom contains the geometry of my entity, it is a polygon geometry.
I want to convert the DataFrame to a GeoDataFrame with :

from shapely import wkt
df.loc[:, 'geom'] = df.loc[:, 'geom'].apply(wkt.loads)

Here is the answer:

ParseException: Unknown type: '01060000206A080000010000000103000000010000001700000026FDBD944CB828410C5A48CD852B5A41F629C7744AB82841BA15C252842B5A412C82FF9548B82841A73E9011842B5A417B9FAAFA22B82841B1355BE08C2B5A411A321E0DEFB728416631B13B952B5A41CA349AE401B828419F58A7909F2B5A411C06F37746B82841F4A5B77AA12B5A41172994DD58B82841BFD71000A72B5A414AD40B1694B8284127DA5543B32B5A417E1B62DCB8B8284177D84486BB2B5A41D540F321DDB82841C7F0D87FBB2B5A413C122F0F07B928419A7B4826B62B5A413A3B191C39B92841516B9A59B42B5A415918225F3EB92841028239FDA62B5A41D2FBC68FBFB9284168D0D0E9AD2B5A415F402F8CD3B92841761728E0AB2B5A418D7DC98EFBB92841096CCE8CA72B5A4103D2FE27FDB92841234BE6CCA32B5A41A67C08029EB92841B5183C119C2B5A41710514BA57B92841DCB8C5AD9C2B5A4125E659B9F8B82841959A3DF2942B5A416D74CE5FF0B828418D2958F3902B5A4126FDBD944CB828410C5A48CD852B5A41'
WKTReadingError: Could not create geometry because of errors while reading input.

I also tried :

df.loc[:, 'geom'] = df.loc[:, 'geom'].map(shapely.wkt.loads)

And I got the same anwser as above.
Any ideas to proceed the conversion ?


texts: Google Earth Engine: Extract rule based values from List and Array

texts:
I'm wondering how to extract values from a ee.List and ee.Array that meet a certain criteria in the Google Earth engine environment.

var array1 = ee.Array([1,1,2,2,3,3,4,4])
var list1 = ee.List([1,1,2,2,3,3,4,4])

I was thinking in the lines of filtering but I don't know the property to filter on, should I .set() that beforehand somehow?

var output = list1.filter(ee.Filter.gt('Item',2))

I assume I should use .slice for the array, but can't find an example that uses slice on a filtered array or conditional statement.

Wondering if somebody can guide me in the correct direction!

Regards


texts: Is there a way to use wms like a tiled map service in sharpmap?

texts:
In sharpmap , is it possible to load a default wms like a tiled map service?

I was able to use a WMS but currently it loads the map as a single image.

I desire a behavior  like OpenLayers 3.

Edit:

My attempt with BruTile:

      const String url = @"http://ows.terrestris.de/osm/service?dummy=dummy"; //dummy is needed because BruTile appends every parameter with ampersand

        BruTile.TileSchema schema = new BruTile.TileSchema();
        schema.Srs = "EPSG:25832";
        schema.Format = "image/png";
        schema.Axis = BruTile.AxisDirection.Normal;
        schema.Name = "UTM32";
        schema.Height = 256;
        schema.Width = 256;
        schema.OriginX = 271930.43;
        schema.OriginY = 5209532.85;
        schema.Extent = new BruTile.Extent(271930.43,5209532.85, 956095.08,6111268.07);

     Dictionary<String, String> parameters = new Dictionary<String, String>();
        parameters.Add("STYLES", String.Empty);

        var request = new BruTile.Web.WmscRequest(new Uri(url), schema, new List<String> { "OSM-WMS" }, new List<String>(), parameters, "1.1.1");
        var provider = new BruTile.Web.WebTileProvider(request);
        var tileSrc = new BruTile.TileSource(provider , schema);

var layer1 = new SharpMap.Layers.TileLayer(tileSrc, "OSM");
            layer1.SRID = 25832;

That is what I have tried so far, but it does not work. Sharpmap shows nothing, not even an exception.

var ti = new BruTile.TileInfo();
ti.Extent = new BruTile.Extent(271930.43, 5209532.85, 956095.08, 6111268.07);
var uri = request.GetUri(ti).ToString();

This gives me a valid wms request, but I get a 407 Proxy Exception when I try the following.

 var ti = new BruTile.TileInfo();
 ti.Extent = new BruTile.Extent(271930.43, 5209532.85, 956095.08, 6111268.07);
 var tile = provider.GetTile(ti);

In Shapmap it is possible to configure a Proxy for SharpMap.Layers.WmsLayer but not for SharpMap.Layers.TileLayer.
Is there another way to set up a proxy?