texts: Issues converting image to numpy array in Earth Engine

texts:
I am trying to build a numpy array fromm the values of an image band in GEE to export for analysis. The roi is quite large (over the South Pacific), so several images make up each day. I have adapted this code for exporting the numpy array, and used this code to produce a mosaic of all the immages for each day in the date range.
My code:
import ee
import numpy as np

poly = ee.Geometry.Polygon(
  [[[124, -60], [124, -10], [310, -10], [310, -60]]], None,
  False)

start = ee.Date('2019-12-30')
finish = ee.Date('2020-01-17')

collection = ee.ImageCollection("COPERNICUS/S5P/OFFL/L3_AER_AI").filterDate(start, finish).select('absorbing_aerosol_index').map(lambda image: image.clip(poly))

# Difference in days between start and finish
diff = finish.difference(start, 'day')

# Make a list of all dates

range = ee.List.sequence(0, diff.subtract(1)).map(lambda day: start.advance(day,'day'))

# Funtion for iteraton over the range of dates
def day_mosaics(date, newlist):

    # Cast
    date = ee.Date(date)
    newlist = ee.List(newlist)

    # Filter collection between date and the next day
    filtered = collection.filterDate(date, date.advance(1,'day'))

    timeBand = ee.Image(date.millis()) \
        .divide(1000 * 3600 * 24) \
        .int() \
        .rename('t')

    # Make the mosaic
    image = ee.Image(filtered.mosaic()).addBands(timeBand)

    # Add the mosaic to a list only if the collection has images
    return ee.List(ee.Algorithms.If(filtered.size(), newlist.add(image), newlist))

# Iterate over the range to make a new list, and then cast the list to an imagecollection
newcol = ee.ImageCollection(ee.List(range.iterate(day_mosaics, ee.List([]))))
# print(newcol)

listOfImages =newcol.toList(newcol.size())
firstImage = ee.Image(listOfImages.get(0))
secondImage = ee.Image(listOfImages.get(1))

# convert to numpy array

band = secondImage.select('absorbing_aerosol_index').unmask(-9999).sampleRectangle(region=poly)

band_arr = band_arrs.get('absorbing_aerosol_index')
np_arr = np.array(band_arr_b4.getInfo())

# convert to nan

np_arr[np_arr == -9999] = np.nan

import matplotlib.pyplot as plt
plt.imshow(np_arr)
plt.colorbar()

The output is a 50 x 186 array-I don't know the original pixel count in the image (or how to find it), but this seems too small. Furthermore, when quicly I try to visualise this array, the result is the following:

It appears to be upisde-down and half of the AOI is missing.
What am I doing wrong?


texts: Why create datastore in GeoServer REST API can not create database on PostgreSQL also?

texts:
I am going to create datastore (PostGIS type) in a workspace by GeoServer REST API. Moreover I want it create database with exact datastore name on PostgreSQL simultaneously. Therefore I add  <entry key="create database">true</entry> to my datastore.xml file (based on GeoServer REST API). It can create datastore in workspace without any error. However it cannot create database same datastore name on PostgreSQL! It seems it did incomplete! This is because if I login to GeoServer web administrator interface when I click store name. I can see form is completed! And everything is checked already but it didn't save. If I press save button it will finish its task completely. And it can create database on PostgreSQL.
<?xml version="1.0" encoding="UTF-8"?>
    <dataStore>
    <name>G0BUXtNQiQ</name>
    <type>PostGIS</type>
    <enabled>true</enabled>
    <workspace>
      <name>G0BUXtNQiQ</name>
      <atom:link xmlns:atom="http://www.w3.org/2005/Atom" rel="alternate" href="http://localhost:8080/geoserver/rest/workspaces/G0BUXtNQiQ.xml" type="application/xml"/>
    </workspace>
    <connectionParameters>
      <entry key="port">5432</entry>
      <entry key="user">postgres</entry>
      <entry key="passwd">postgres</entry>
      <entry key="dbtype">postgis</entry>
      <entry key="host">localhost</entry>
      <entry key="database">G0BUXtNQiQ</entry>
      <entry key="schema">public</entry>
      <entry key="create database">true</entry>
    </connectionParameters>
    <__default>false</__default>
  </dataStore>

This is log file:

2020-09-07 17:39:46,990 WARN [geofence.cache] -
org.geoserver.geofence.cache.CachedRuleReader$NoAuthException: Can't
auth user [admin] 2020-09-07 17:39:47,045 INFO
[auth.GeoFenceAuthenticationProvider] - User admin NOT authenticated
2020-09-07 17:39:49,127 WARN [geofence.cache] -
org.geoserver.geofence.cache.CachedRuleReader$NoAuthException: Can't
auth user [admin] 2020-09-07 17:39:49,287 WARN [geoserver.security] -
Failed login, user admin from 0:0:0:0:0:0:0:1 2020-09-07 17:39:49,287
INFO [geoserver.security] - Brute force attack prevention, delaying
login for 1842ms 2020-09-07 17:39:51,695 INFO [geoserver.security] -
Strong cryptography is available

SYSINFO: GeoServer 2.17.2,
PostgreSQL 12.4 (Debian 12.4-1.pgdg100+1) on x86_64-pc-linux-gnu, compiled by gcc (Debian 8.3.0-6) 8.3.0, 64-bit


texts: Accessing GeoServer with encrypted password?

texts:
In Accessing secure Geoserver layer with username password I have found how to access the GetCapabilities of a protected service adding the user/password in the URL: The answer works perfectly:

http://user:password@localhost:8080/geoserver/wms?request=GetCapabilities

The problem is that I can't put the user/password information in a public metadata file for obvious reasons.

In the same question, the user who answered it said that the password can be encrypted. That would be great.

How can I get my encrypted password of GeoServer? I know that in "Settings" you can check the "Encrypted" box. What I want to know is the string of my encrypted password to add it to the URL.

The GeoServer documentation (http://docs.geoserver.org/latest/en/user/security/passwd.html) says it will look like this:

crypt1:KWhO7jrTz/Gi0oTQRKsVeCmWIZY5VZaD. crypt1

It doesn't say which tool I should use to get the string (only the encrpytion type, and I don't know what to do with that).


texts:
I have a shapefile with Mediterranean countries, and an attribute table such as the one below:

In an excel spreadsheet, I have population data for all of these countries, which I'd like to assign.

Is there a way I can export my excel data so that population of Italy is assigned directly to Italy in the attribute table for example?


texts: One way Replica deleted feature not synchronizing when I synchronize

texts:
I have created one-way replicas from 50 geodatabases to one geodatabase (Arcsde. With PostgSQL 9.2) After that, I update my feature class in parent geodatabase, but when I synchronize the replica the deleted feature has not been synchronized to the child geodatabase.  What can be done to make the synchronize work?

  here is replica log :


texts: Python Add-in using arcpy.mapping.ListLayers() does not update when run multiple times

texts:
I am developing a Python Add-in for ArcMap 10.4 that, upon a click on a map, opens a wxPython GUI that contains a combo box listing the layers in the map document.

Here's a current snippit from my script:

class ToolClass2(object):
"""Implementation for DesignTool.tool (Tool)"""
dlg= None
def __init__(self):
    self.enabled = True
    self.shape = "NONE" # Can set to "Line", "Circle" or "Rectangle" for interactive shape drawing and to activate the onLine/Polygon/Circle event sinks.
    self.cursor=3 #set cursor to crosshair
def onMouseDownMap(self, x, y, button, shift):
   global xy
   xy = (x,y)
   global layerlist
   layerlist=[]
   mxd=arcpy.mapping.MapDocument("CURRENT")
   for lyr in arcpy.mapping.ListLayers(mxd):
       layerlist.append(lyr.name)
   del mxd

I then reference layerlist within my wxpython class to have the layers display as choices within the combo box. When I run it for the first time after opening ArcMap, it runs correctly and displays all of the layers within the map document. However, when I attempt to run it again after adding or deleting layers, it does not update the list - even if I start a new blank document.

How can I write my script to update upon each run of the tool?

I am just assuming that this is an issue with my map document, but perhaps it is an issue with my WxPython dialog?

UPDATE: I now believe that this issue is coming from the fact that the self.Destroy() command in wxPython does not work in ArcMap. As a workaround for this, I was using self.Show(False) to close the UI (similar to an example shown here: https://geonet.esri.com/thread/38336), but when run again, self.Show(True) displays the same UI, as it was never actually destroyed. Are there any known alternatives to self.Destroy() that will actually get rid of the UI and not just prevent it from being displayed?

import arcpy
import pythonaddins
import wx

class ToolClass2(object):
    """Implementation for DesignTool.tool (Tool)"""
    dlg= None
    def __init__(self):
        self.enabled = True
        self.shape = "NONE" # Can set to "Line", "Circle" or "Rectangle" for interactive shape drawing and to activate the onLine/Polygon/Circle event sinks.
        self.cursor=3 #set cursor to crosshair
    def onMouseDownMap(self, x, y, button, shift):
        global laylist
        laylist=[]
        mxd=arcpy.mapping.MapDocument("CURRENT")
        for lyr in arcpy.mapping.ListLayers(mxd):
           laylist.append(lyr.name)
        del mxd
        if self.dlg is None:
            self.dlg = TestDialog()
        else:
            self.dlg.Show(True)
        return

class TestDialog(wx.Frame):
    def __init__(self):
        wxStyle = wx.CAPTION | wx.RESIZE_BORDER | wx.MINIMIZE_BOX |wx.CLOSE_BOX | wx.SYSTEM_MENU | wx.CB_DROPDOWN
        wx.Frame.__init__(self, None, -1, "Design Form", style=wxStyle, size=(330, 370))
        self.SetMaxSize((330, 370))
        self.SetMinSize((330, 370))
        self.Bind(wx.EVT_CLOSE, self.OnClose)
        panel = wx.Panel(self, -1)
        wx.StaticText(panel, -1, "Choose Layer:", pos=(8,64))
        self.LayerCombo = wx.ComboBox(panel, -1, value=laylist[1], pos=(180, 64), size=(120,21), choices=laylist)
        self.Bind(wx.EVT_BUTTON, self.OnSet, id=self.btnSet.GetId())
        self.Show(True)

    def OnClose(self, event):
        self.Show(False) # self.Destroy() doesn't work

    def OnSet(self, event):
        Layerpath= str(self.LayerCombo.GetValue())
        self.Show(False)

app = wx.PySimpleApp()
app.MainLoop()


texts: Changing the geolocation of a point in the map

texts: I uploaded a vector layer with points = cities. After zooming in I found out that some cities are on wrong places. How can I change the geolocation for a specific point, city?


texts:
I recently unearthed a GPS unit I purchased a few years. It's a Garmin GPSmap 60 (picture) and has a usb port on it. Is there a way that I can link it to my mac? I do have windows installed also, and could use that if necessary.

I tried plugging it in, and downloading data via google earth pro without success. I also hunted around online to see what I should install but did know what I was searching for exactly (a driver? a software program? ).


texts: Compacting file geodatabases for editing performance?

texts:
Is there a way to quantify the effect of compacting your ESRI file geodatabase(s) on editing performance? I currently use ArcInfo(10.0).

I have a 769 MB fgdb that holds my most used/edited datalayers either in feature datasets or stand alone feature class(es) and am considering a redesign of how I have things set up for my most common workflows and editing tasks. Other than this link to the pertinent ESRI's help page, I really haven't been able to find anything pertaining to this question.


texts:
I want to calculate FAPAR in different land cover types using max and min value of NDVI.

I have a raster image with 9 Land cover types and max, min value of NDVI in each land cover type but I don't know the way how to calculate it in different LC types from NDVI image.

What kind of function can be used to do this?
What software is more useful for this kind of calculation?


texts: PyQGIS scripts outside of QGIS GUI running processing algorithms from GRASS provider

texts:
I'm writing a script that I would like to work outside of QGIS GUI. The script uses the processing algorithm standalone just fine for the native algorithms, however I am experiencing errors when trying to run grass algorithms. The code I am currently using to initialise the environment come from the following post:
Using QGIS3 Processing algorithms from standalone PyQGIS scripts (outside of GUI)

import sys
from qgis.core import QgsProcessingRegistry
from qgis.analysis import QgsNativeAlgorithms
from qgis.core import (
     QgsApplication
)

QgsApplication.setPrefixPath('/usr', True)
qgs = QgsApplication([], False)
qgs.initQgis()

# Append the path where processing plugin can be found
sys.path.append('C:/OSGeo4W64/apps/qgis/python/plugins')

from processing.core.Processing import Processing
Processing.initialize()
QgsApplication.processingRegistry().addProvider(QgsNativeAlgorithms())
import processing

This code works to initialise the environment and allows me to run QGIS native processing operations in standalone scripts. The issue arises when i try to run the following code:

# Grass Process
parameters = {
    '-f': False,
    'GRASS_RASTER_FORMAT_META': '',
    'GRASS_RASTER_FORMAT_OPT': '',
    'GRASS_REGION_CELLSIZE_PARAMETER': 0,
    'GRASS_REGION_PARAMETER': None,
    'areas': 'C:/TEMP/areas.tif',
    'direction': 'C:/TEMP/direction.tif',
    'format': 0,
    'input': 'C:/TEMP/DEM.tif',
    'output': 'C:/TEMP/DEM_Filled.tif'
}
processing.run('grass7:r.fill.dir',parameters)

The process runs successfully from within QGIS running from the python console, however when I run the script externally I get the following error from processing:

File "C:/OSGeo4W64/apps/qgis/python/plugins\processing\tools\general.py", line 96, in run
    return Processing.runAlgorithm(algOrName, parameters, onFinish, feedback, context)
File "C:/OSGeo4W64/apps/qgis/python/plugins\processing\core\Processing.py", line 181, in runAlgorithm
    raise QgsProcessingException(msg)
_core.QgsProcessingException: There were errors executing the algorithm.

I am running the scripts in the python environment initialised by: C:\OSGeo4W64\bin\python-qgis.bat


texts:
After installation of Java SDK for ArcGIS version 10.0, authorization of the SDK I still can't run the simple code in agent analyst:

def IWorkspace():

self.IAoInitialize = AoInitialize()

self.IAoInitialize.initialize
(self.esriLicenseProductCode.esriLicenseProductCodeArcInfo)
self.IAoInitialize.shutdown();

self.IWorkspaceFactory = FileGDBWorkspaceFactory();

self.IFeatureWorkspace = IFeatureWorkspaceProxy(self.IWorkspaceFactory.openFromFile(self.dataPath, 0));

when running the agent analyst it fails executing last line:

Exception in thread "Thread-14" AutomationException: 0x80040228 - The Product License has not been initialized. in 'ESRI GeoDatabase'
    at com.esri.arcgis.datasourcesGDB.FileGDBWorkspaceFactory.openFromFile(Unknown Source)
    at default_package.GISModel.IWorkspace(Unknown Source)
    at default_package.GISModel._schedule2(Unknown Source)
    at BAB_SYNTH_6.execute(Unknown Source)
    at uchicago.src.sim.engine.ScheduleGroup.execute(ScheduleGroup.java:81)
    at uchicago.src.sim.engine.RandomScheduleGroup.execute(RandomScheduleGroup.java:52)
    at uchicago.src.sim.engine.Schedule.execute(Schedule.java:349)
    at uchicago.src.sim.engine.BaseController$1.run(BaseController.java:397)
    at java.lang.Thread.run(Thread.java:619)

What should I do?


texts: Classification error "Layer error: Property 'class' of feature '00000000000000000000_0' is missing"

texts:
I am trying perform a Random Forest classification using an uploaded shapefile of polygons with only one class = 0.  Get error Cover crop classification: Layer error: Property 'class' of feature '00000000000000000000_0' is missing.

Suspect not reading the "class" field (attribute) in the shapefile correctly, but the field prints.

https://code.earthengine.google.com/?scriptPath=users%2Ftracyd%2FDelawareChange%3ACoverCropClassification


texts: Postgis spatial join performance (point on area)

texts:
I need to update a single attribute of a large table of point data with a result of a spatial join (point on area operation) on a polygon table. I'm using the below query but wondering if this is the best performance approach I could take.
I have geometry indexes on both tables and this is a one time operation.

update point_table pt set name = po.name
from polygon_table po where ST_contains(po.geom, pt.geom);


texts: Why is ellipsoidal distance greater than cartesian distance?

texts:
I created some linestring features in EPSG:32632 (WGS 84 / UTM zone 32N) exactly 10 000 m long.

QGIS calculates all line lengths (horizontal, vertical, diagonal) as 10 001.9 m. The measure tool explains "Project CRS transformation is turned on and ellipsoidal calculation is selected. The coordinates are transformed to the chosen ellipsoid (WGS84), and the distance is calculated in meters."

When I set the ellipsoid to "None / Planimetric" in project properties, I get the expected results (10 000 m).

I don't understand the phrase "coordinates are transformed to the chosen ellipsoid (WGS84)", is it a transformation to unprojected degrees (4326?)? Since 32632 already is based on WGS84, is there anything to transform? If there is a kind of "great circle" calculation involved, I would have expected lengths can only get smaller.

Is QGIS's calculation correct/meaningful, is it plain wrong or do I see rounding errors?

This are the geometries I tested:

LineString (370000 5615000, 370000 5625000)
LineString (366464.46609406732022762 5616464.46609406732022762, 373535.53390593267977238 5623535.53390593267977238)
LineString (365000 5620000, 375000 5620000)
LineString (373535.53390593267977238 5616464.46609406638890505, 366464.46609406638890505 5623535.5339059317484498)

Edit as addition to MappaGnosis's answer: I forgot that in UTM the scale is less than 1 inside the standard lines, this seems to be the simple answer for the title question.


texts: Sen2Like GDALRasterizeOptions Failed to process SRS definition

texts:
While running sen2like tool in single-tile-mode an error related to osgeo gdal occurred (Python):
[INFO ] 2021-04-07 13:01:06 - reader - Sentinel2MTL Class
[INFO ] 2021-04-07 13:01:06 - reader - Product: sen2like-master/data/Sentinel2/33TWF/S2A_MSIL1C_20200313T095031_N0209_R079_T33TWF_20200313T102505.SAFE
[INFO ] 2021-04-07 13:01:10 - sentinel2 - SAT_AZ , SAT_ZENITH, SUN_AZ, SUN_ZENITH
[INFO ] 2021-04-07 13:01:10 - sentinel2 - UNIT = DEGREES (scale: x100) :
[INFO ] 2021-04-07 13:01:10 - sentinel2 - /home/ec2-user/wd/14572/S2A_MSIL1C_20200313T095031_N0209_R079_T33TWF_20200313T102505.SAFE/tie_points.tif
[INFO ] 2021-04-07 13:01:10 - sentinel2 - Generating nodata mask from band B01
[INFO ] 2021-04-07 13:01:11 - image_file - Written: /home/ec2-user/wd/14572/S2A_MSIL1C_20200313T095031_N0209_R079_T33TWF_20200313T102505.SAFE/nodata_pixel_mask_B01.tif
[INFO ] 2021-04-07 13:01:11 - sentinel2 - Generating validity mask from cloud mask

Traceback (most recent call last):
File “sen2like-master/sen2like/sen2like/sen2like.py”, line 435, in
sys.exit(main(with_multiprocess_support=True))
File “sen2like-master/sen2like/sen2like/sen2like.py”, line 430, in main
start_process(tile, _products, args, start_date, end_date)
File “sen2like-master/sen2like/sen2like/sen2like.py”, line 253, in start_process
process(_product, args.bands)
File “sen2like-master/sen2like/sen2like/sen2like.py”, line 284, in process
product.mtl.get_valid_pixel_mask(os.path.join(config.get(“wd”), product.name, ‘valid_pixel_mask.tif’))
File “/home/ec2-user/SageMaker/sen2like-master/sen2like/sen2like/core/readers/sentinel2.py”, line 404, in get_valid_pixel_mask
xRes=res, yRes=res)
File “/home/ec2-user/SageMaker/custom-miniconda/miniconda/envs/sen2like_ker/lib/python3.6/site-packages/osgeo/gdal.py”, line 1163, in Rasterize
(opts, callback, callback_data) = RasterizeOptions(**kwargs)
File “/home/ec2-user/SageMaker/custom-miniconda/miniconda/envs/sen2like_ker/lib/python3.6/site-packages/osgeo/gdal.py”, line 1150, in RasterizeOptions
return (GDALRasterizeOptions(new_options), callback, callback_data)
File “/home/ec2-user/SageMaker/custom-miniconda/miniconda/envs/sen2like_ker/lib/python3.6/site-packages/osgeo/gdal.py”, line 3530, in init
this = _gdal.new_GDALRasterizeOptions(*args)
RuntimeError: Failed to process SRS definition: EPSG:32633

I don’t understand the error, the installed version of gdal is the 3.0.1 as required.
gdal_rasterize documentation says “The <srs_def> may be any of the usual GDAL/OGR forms, complete WKT, PROJ.4, EPSG:n”. So the parameter should accept “EPSG:32633”
Any idea?


texts: Simulating urban land use change using a Cellular Automata - Markov Model in ArcGIS

texts:
I am trying to simulate urban land use change using a Cellular Automata - Markov Model.

Is it possible to do this in ArcGIS (or any ArcGIS extension)?

I do not have access to IDRISI. I am fine with using any ArcGIS extension or other free software. My urban data has 5-7 land use classes that I need to simulate for. The output should be in map form.


texts:
I have downloaded a series of Landsat-8 level 2 images. My main purpose is to obtain reflectance out of them; however, the scale of the different bands does not look like reflectance (e.g. Blue Band goes from -2000 to 9393).

I am assuming those numbers correspond to digital numbers (DN) (I am not 100% sure), so I am wondering if there is a procedure to convert this DN to reflectance for this level 2 data.


texts: Labeling line Nodes with m and z values in QGIS

texts: Is there a way to label the nodes of line features with the z and m values using QGIS?


texts:
ArcGIS 10.2, ESRI Maps for MS Office
Why a heatmap that was created in Excel and looks fine and is still fine on a preview of a map if I publish it it ArcGIS Online  but turns into black points if I open it in ArcGIS Online map viewer ?


texts: use proj4 to specify Robinson projection with R ggmap and ggplot2 packages?

texts:
I want to project this map in robinson projection:

library(ggmap)
world <- map_data("world")
ggplot() + geom_path(data = world,
                              aes(long, lat, group = group))

And I would like to change the projection to "Robinson" (following advice from answer to my previous question: What projection does the global climate region map from Wikipedia use?

I had a hard time finding a default implementation of this projection, I worked out the following for using the proj4 library:

library(proj4)
robinson <- project(cbind(world$long, world$lat),
                    proj = "+proj=robin +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs")

I have tried a number of approaches, including:

# using ggmap::get.map()
get_map("world", projection = mapprojection(robinson))
# using ggplot2::coord_map
coord_map(projection = robinson)
# and sp::coordinates:
library(sp)
coordinates(world) <- ~ lat + long
gridded(world) <- TRUE # returns error
proj4string(world) <- CRS(robinson)

but none of these work. Is it a typo, or I missing something fundamental about this method?


texts: Low-cost centimeter-level accurate satellite positioning (GNSS/GPS)

texts:
I am looking for a low-cost (< US$500) GNSS/GPS setup with a centimeter precision without much HW hacking. I am not able to produce my PCB or do any soldering (though I would do that if there is no other way) so a kind of a easy-to-assemble setup would be welcome. I know about the $900 Piksi thing but that is still too expensive for me. It seems like cm precision should be possible for much less - like employing a 50 USD raw GPS sensor with an antenna and ordinary PC with RTKLIB software.

I am not sure if it is better to use two GPS sensor setup for RTK (one base station and one for rover) or whether I can get the corrective DGPS data elsewhere (my region is Czech Republic - there seems to be national grid here allowing to stream correction data for reasonable cost).

My application will be in a passenger car so I will not be limited with power source - no low power needed although that would be nice. I will be using the position readings within OpenCV - so I need to get the data into C/C++ code. The application is data collection so I can use raw GPS post-processing.

EDIT: I was recommended a $150 U-BLOX LEA-6T USB board that would make a $300 sum for dual unit setup without the host HW (probably 2 x PC or Raspberry Pi?) and software. As software RTKLIB could probably be used although I have zero experience with either (HW & SW) of these.

Also, there seems to be the NAVIO board from EMLID in development and hopefully to be released soon. This will probably offer an "out of the box" RTK solution (both HW & SW) at an estimated price of $165 for board (see comments on that link). They made some RTK demonstration (theory here) and also some Real-time OS tests.


texts: PyProj - Transform coordinate in "HK80 grid" to "WGS 84 / UTM zone 49N or 49S"

texts:
Am trying to use python pyproj to do the above conversion (convert/transform from Hong Kong GRID to Hong Kong UTM), but am getting different values from the expected result.

My expected conversion result is: "North 836148, East 822745" to "North 828672.7449, East 12488562.2742". See my codes below for both zones (N and S):-

# https://epsg.io/102141

inProj = Proj(init='epsg:2326') # PCS: Hong_Kong_1980_Grid or EPSG:102140
outProj = Proj("+proj=utm +zone=49 +ellps=intl +units=m +no_defs") # PCS: Hong_Kong_1980_UTM_Zone_49N or EPSG:102141
# outPro = Proj(init='epsg:102141') # This returns RuntimeError

# Expected result: 828672.7449, 12488562.2742

x1, y1 = 822745, 836148
# x1, y1 = 836148, 822745

x2, y2 = transform(inProj, outProj, x1, y1)
print (x2, y2)

Or

# https://epsg.io/32749

inProj = Proj(init='epsg:2326') # PCS: Hong_Kong_1980_Grid or epsg:102140
# outProj = Proj("+proj=utm +zone=49 +south +datum=WGS84 +units=m +no_defs") # PCS: Hong_Kong_1980_UTM_Zone_49S or EPSG:32749
outProj = Proj(init='epsg:32749')

# Expected result: 828672.7449, 12488562.2742

x1, y1 = 822745, 836148
# x1, y1 = 836148, 822745

x2, y2 = transform(inProj, outProj, x1, y1)
print (x2, y2)

Could it be that my expected result is incorrect or I am missing something?
I want to certisfy the expected output by going from: "Hong_Kong_1980_Grid_System" to "Hong_Kong_1980_UTM_Zone_49N / 49S"?

North   East    North1      East1
836148  822745  828672.7449 12488562.27
831905  822495  827480.7854 12473200.37
838208  816918  810095.1039 12481647.01
838127  816314  805922.5972 12482194.51
837452  820146  815300.3898 12484301.71
833764  816330  808333.7182 12478896.42
841476  830964  834944.8892 12495195.37

The result needs to match above table from Hong Kong 1980 grid to WGS84/UTM 49S


texts: Delete several tables in QGIS DB manager

texts: I have a lot of tables in a PostGIS-Database and I want to delete most of them. In the QGIS DB manager, I am able to delete one by one (right-click -  delete). However, that's quite annoying. How to delete several tables in the same time.


texts: open-source code for the map matching algorithm

texts:
For map matching of the GPS data to the network data, there is a algorithm from Schussler, N. & Axhausen, K. (2009a).
Is this algorithm released under open-source license ? Where can I find this ?

"I know SF City implemented one of these algorithms to process gps tracks from a smart phone app they are using to track cyclists (cycletracks); and they released all of the code under an open-source license."- found in the below question.

What algorithm should I use to remove outliers in trace data?


texts: Method to determine of coordinates in GeoJSON feature are Lat/Lng or UTM?

texts: Is there a @turf or any other helper method to quickly determine if coordinates in GeoJSON feature are Lat/Lng  or UTM?


texts: Why does a coordinate system transformation warp distances?


texts: INASAFE using US Census data? or similar QGIS hazard analysis - thresholding?

texts:
INASAFE seems like an interesting QGIS plugin but as I wade through the docs, I really need some sample data to play with in order for me to understand.  There is a tutorial which uses Indonesia data, but I can't find that dataset either.

If anyone is using this and has any example datasets, it would be very appreciated.

OK, after working thru the INASAFE tutorial it was unclear how to use Census data, so I waded thru the python scripts.  It became clear that INASAFE couldn't deal with my data.  I'm not familiar enough with either QGIS or python to add my capability to INASAFE so I decided to try my luck doing a similar analysis using QGIS (newbie here, but a bold one).

INASAFE appears to work well with OSM building data and Indonesia has done a terrific job mapping what appears to be thousands of their buildings in OSM, but my county has very few and census data has much more detailed data.

So I took my 2010 census block data, which had fields of POP10 for # of people in block and HOUSING10 for # of housing units in block and added a calculated AREA field.  Since I had everything in Google Mercator, this area is in some nonsense units, but I don't care about that.  As long as cutting the polygon in half gives an area half as big, then it will work for my purposes.

Next I did an intersect with some circa 2002 FEMA 100-yr floodplanes (remember when FEMA floodplanes were free!).  On the intersect layer I added a calculated AREA field called FloodAREA.

I calculate the flooded Population as POP10*FloodAREA/AREA (ratio of area of flooded polygon to original polygon).   As expected, some times FloodAREA=AREA (completely flooded block), and sometimes was less (partially flooded block), but never more (sanity check).

I ran basic statistics on the intersect layer and got the surprising result that 45,000 people (10% of pop) in my county live in 100 year flood plane.

I was wondering about thresholding - if polygon has like 0.1 people flooded, probably shouldn't count it - people and houses are not uniformly distributed in block.  Especially in flood plane, I suspect lowest spots are less likely to have houses built on them.  So, if a block has 250 people, for example, and it shows 25 people flooded - this really is saying that lowest 10% of block is flooded. Housing density rarely seems to approach 50%, so it's very possible (likely?) that when 10% of a block is flooded, there are not 10% of the houses in the flooded area.

So 2 questions:

Are there any ideas or standards out there for discarding partial people from the sum or requiring minimum % of block flooded before counting?
How do I get a sum of a field in QGIS with a thresholding formula applied (like sum(if Flooded_Pop>1 and FloodAREA/AREA>0.05)?


texts:
I tried to merge two DEM-rasters without any overlapping parts using "Mosaic To New Raster" tool. However the value range of output raster doesn't look correct. Assuming that (A1,B1) and (A2,B2) are value ranges of the first and the second raster respectively, I was expecting value range of output be (MIN(A1,A2),MAX(B1,B2)) but wasn't. I tried merging with different tool settings but I couldn't get expected result.
Is my expectation wrong? If not where do I make mistake?


texts: Web app - how to replace ESRI branding with a custom graphic?

texts:
Alough there apears to be an easy way to remove the "Powered by ESRI" graphic in JSON I would like, rather than removing it, to substitute it with a custom PNG.

When I set "Logo": false, in "mapOptions" the Esri logo turns off but if I set the logo to an existing valid PNG, rather than displaying it, the custom graphic is somehow over-ridden and the "Powered by ESRI" graphic is displayed.

Any ideas how to replace the default with a custom graphic?

Default graphic:

Setting the logo to a valid PNG has no effect:

Setting the logo to false works!

Result:


texts: Access ArcGIS Server Rest Service with jQuery

texts:
I am trying to use jQuery to access an ArcGIS Server rest service. I keep getting:

Request header field X-Requested-With is not allowed by Access-Control-Allow-Headers in preflight response

ESRI support swears that it is a jQuery problem because they can access those web services via their API, which is Dojo.

From my study on the web it appears to be an issue of headers. In my JavaScript I have tried numerous combinations of headers, for example:

$.ajax({
     type: "Get",
     //url: "get_esri_feature_service.php",
     url: "https://geohub.vancouver.ca/arcgis/rest/services/GDH/School/FeatureServer/0/query?where=1%3D1&geometryType=esriGeometryEnvelope&spatialRel=esriSpatialRelIntersects&units=esriSRUnit_Foot&outFields=objectid%2Cschoolname%2Caddress&returnGeometry=true&f=json",
     dataType: "json",
     crossDomain:true,
     headers: {'Access-Control-Allow-Origin': '*',
               'Access-Control-Allow-Methods': 'POST,GET,OPTIONS',
               'Access-Control-Allow-Headers': 'Origin, X-Requested-With, Content-Type, Accept'
              },
     success: getFeaturesSuccess,
     error: ajaxError
 });

Still it doesn't work. Any suggestions?


texts: Georeferencing PNG in ArcMap based off of data frame corners?

texts:
First, in Arcmap I set the data frame exactly equal to the page layout (2'hx3'w) and exported my hillshade, DEM, etc one after the other from layout view so that they were all the exact same size and extent.

Then, I did my editing in Photoshop to create a bump map/hillshade backdrop.

Now I am trying to bring the image back in as a backdrop for a trailmap. I don't want to georeference based off common features if I can just use the exact corners of the data frame from my original export arcmap (I am already familiar with the process of georeferencing using common features). It seems like using the exact corners would prove more accurate and this process would be easier to streamline for future users.

If I am understanding correctly I can find the exact corners of the data frame in the Data Frame Properties>Data Frame>Extent>Fixed Extent. Is there a way to apply these coordinates to the corners of my image?


texts: ArcGIS Interpolation IDW and Spline, extreme values

texts:
I am trying to interpolate air pollution data that is displayed in 1x1 Km points in ArcGIS. I am having two issues:

When I try to interpolate using IDW, the software freezes at 94% "building search tree," and I have had not been able to produce an IDW surface
When I try to interpolate using the spline method, a raster is executed but with extreme negative (-91) and positive (82) values that are not within my data range (0-17)


texts: Create a 3D topographic map with locations marked on the map using R

texts:
I am currently trying to create an interactive 3D surface plot using plotly in R. My dataset includes a matrix of Latitude and Longitude and I have a DEM for the region I would like to make this surface plot for. However, I am unable to do the same.

I have a dataframe Points.sp that contains a list of Lat/Long along with elevation data for areas where I have sampled.

> head(Points.sp)
                         Points     Long      Lat Elevation
1                                    NA       NA        NA
2                      Nair212       NA       NA        NA
4                      CESF543 75.72130 12.88470   876.000
5                     CESF1587 75.93559 11.93745  1494.523
6          tysoni_WILD15AMP650 75.82200 12.38000  1051.000
7           tysoni_SDBDU201274 75.65570 12.22010  1176.000

Here's what I have tried:

Points_new <- data.frame(Points.sp$Long, Points.sp$Lat, Points.sp$Elevation)
Points_new <- as.matrix(Points_new)
p <- plot_ly(z = ~Points_new) %>% add_surface()
p

In the above plot, I manually loaded the elevation values for each of these points, but the problem being that I need these points on a DEM for the region, and everytime I try to load the DEM, I get an error saying - 'z' must be a numeric matrix.

The DEM here is a raster for the region and I obtained it from the SRTM dataset and cropped it to contain my region only.

 elev <- raster("E:\\MasterThesis\\Absence Approach\\Data\\Elevation_Data\\alt\\alt")
 cr1 <- crop(elev, extent(WG.new)) #WG.new is a mask for a smaller region
 newelev<- mask(cr1, WG.new)

 > newelev
  class       : RasterLayer
  dimensions  : 1577, 624, 984048  (nrow, ncol, ncell)
  resolution  : 0.008333334, 0.008333334  (x, y)
  extent      : 72.87501, 78.07501, 8.11667, 21.25834  (xmin, xmax, ymin, ymax)
  coord. ref. : +proj=longlat +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +no_defs
  data source : in memory
  names       : alt
  values      : -14, 2625  (min, max)
  attributes  :
         ID COUNT
  from: -431     2
  to  : 8233     1

If I try loading just the DEM alone by converting the raster to a matrix, I am able to view the DEM, but now I need to add the Lat and Long for locations from the Points.sp dataframe.

newelev <- as.matrix(newelev)
p <- plot_ly(z = ~newelev) %>% add_surface()
p

Any suggestions? I would ideally want points as small circles on the DEM for the region along with the DEM smoothed and stretched out.

*****UPDATE*****

Based on @wittich's elegant solution, I tried what he posted, but I am getting an image like the screenshot posted below.


texts: Add legend in a SplitPanel in Google-Earth Engine

texts:
I have made a code with a Split Panel (1) , now i what to add a legend like this example:

https://code.earthengine.google.com/b5d35a19689ec5698e59f36a9b3ae6f8

this is my original code that I want to use:

(1) https://code.earthengine.google.com/e1438ab2c27221380f9b7c68c8c8b0a1


texts: Information of particular locations (points) in a map

texts:
For a marathon race I need to create a map with specific information about material needed in particular positions (eg crossroads). It's about 150 positions. You see one part in the picture below.

Untill now my company have used Apple Pages and we had to adapt every changes by hand. I would now like to use GIS.

I'm used to ArcGIS which is too expensive to afford for the company I work for. Thus, I came to QGIS.

Is it possible to create such a map in QGIS or how would you recommend to proceed?

Honestly, I am not even sure how to do it in ArcGIS.


texts: Removing roads from Google satellite hybrid map in OpenLayers?

texts:
I am trying to display Google satellite map in Openlayers. So, I added Google Satellite Hybrid map. My problem is, I need to display labels(like POI, place name etc) but don't need any roads & its lables.

So, I just tried this link & I am able to get exact output like below image. In this link they are using Google Map API. So, its working fine.. But I am using Openlayers.

For Street map, I did using below code and working fine..

var map;
var stylez = [{
        "featureType" : "road",
        "stylers" : [{
                "visibility" : "off"
            }
        ]
    }
];

var gmap = new OpenLayers.Layer.Google("Google Streets", {
        type : 'styled'
    });

var styledMapOptions = {
    name : "Styled Map"
};

var styledMapType = new google.maps.StyledMapType(stylez, styledMapOptions);

function init() {
    map = new OpenLayers.Map('map');
    map.addLayers([gmap]);

    map.setCenter(new OpenLayers.LonLat(51.523913, 25.256328).transform(
            new OpenLayers.Projection("EPSG:4326"),
            map.getProjectionObject()), 14);

    gmap.mapObject.mapTypes.set('styled', styledMapType); //Here setting the cusom style
    gmap.mapObject.setMapTypeId('styled');

}


texts: To Iterate through all layers and all attributes for each layer

texts:
I am working in the Python console in QGIS 3.10. I want to iterate through all layers and for each layer iterate through all attributes and set them all to 'Hidden'.

I want to do this in preparation for creating a webmap using a plugin. The data sets I am using have come from various places - GSI, BGS, OSI, EMODNet - and there are a lot of attribute fields in the 14 layers.
I am new to Python and have managed to list the layer names in the console, but am failing to loop through the attributes - I keep getting attribute values - thousands of them.

Can anyone help explain how I can get to list the attribute names and thereby amend their widget to 'Hidden'?

project = QgsProject.instance()
for id , layer in project.mapLayers().items():
    for field in layer.fields():
        print(layer.name() , field.name())

OK, not to worry, I found the answer.

The following code worked,

project = QgsProject.instance()
for id , layer in project.mapLayers().items():
    for field in layer.fields():
        print(layer.name() , field.name())
        idx = layer.fields().indexFromName(field.name())
        setup = QgsEditorWidgetSetup('Hidden', {})
        layer.setEditorWidgetSetup(idx, setup)

This question was cancelled because it was deemed off topic, which I don't get, but never mind.


texts: Converting Entity to GeoJSON with Cesium

texts:
I know there is a way to load entities from a GeoJSON using GeoJsonDataSource but is there any way to do the opposite: get any entity as a GeoJSON representation?

My goal is to use Turf.JS with Cesium.JS geometries.


texts: How can I select points within a part of a raster?

texts:
I have a map with huge amount of points (9562 villages) and two rasters: population density and land elevation.

I am trying, without success for the moment, to select just those villages with high population density and low land elevation.

Does anyone have the solution?

I tried to do it using a site suitability analysis, but for that I need the v.kernel GRASS command which is not working on my pc (another problem described here: v.kernel is missing under Processing Toolbox GRASS GIS 7 commands).


texts: Why are there differences in output position between the ArcGIS versus NCAT transformation [NAD83 to NAD83 HARN]?

texts:
My question is specifically referencing the file-based grid transformation methods that are now available in ArcGIS. I thought this grid-file method was exactly the same method that NGS's NCAT used when moving from NAD 83 to NAD 83 HARN, however I am not getting the same outputs. They are very close and seem to always been in the correct direction, but ArcMap's outputs fall short when compared to my outputs from NCAT. Either way, these state-specific methods (at least in my state) appear to be far more accurate than any others that are listed in the transformation drop-down menu.

Since this doesn't require concatenating transformations at all, why would Arc's outputs be different than NCAT's?


texts: Detecting "lack of movement" on GPS device

texts:
We're using a GPS tracking device mounted in vehicles and later display the movements in our application.
We're having a small issue however with vehicles when they are NOT moving.

Despite standing still, the GPS device will send information with varying positions, sometimes it'll even report the vehicle is moving (every once in a while it'll report a vehicle is moving quite fast). The end result is an ugly representation of the movement in our application, and worse, incorrect calculations with regards to how long the vehicle was moving and was standing still (we've got some statistical analysis going on as well).

I know the problem isn't new: Google Navigation on my Android phone also has trouble when I stop at a crossroad every now and again, detecting that I'm now driving in the opposite direction (when in fact, I'm standing still).

But we really need to have some method of telling the vehicle isn't moving, especially for prolonged periods of time.

The problem is compounded in some vehicles which are kept under roofs for the night, which causes the GPS to go wild due to a weaker signal (which is still strong enough to get a position however).

We can detect if the vehicle's engine is turned on or off, but we cannot assume it's not moving with the engine off (there were incidents where a damaged vehicle was towed... and of course there's a case of potential theft which also needs to be reported).

What's the best way to approach this problem?

EDIT:

To add some more information:

The GPS device does have an accelerometer, but that only returns a binary information (moving / not moving) and it's either TOO sensitive or just plain isn't working.
We do have access to things like number of satellites or quality, but we're unsure how to utilize that information. Thus my question. :)
The problem isn't about detecting if a vehicle is moving or not in real-time. We collect the data, and later on do some statistical analysis and display it. We do show the current position of the vehicle, but that is of little importance. So basically we need to be able to tell a vehicle was or wasn't moving by looking at historical data.


texts: Using ArcMap add x y data, points are in wrong location

texts:
According to ArcGIS esri when I add a CSV with the "add XY" tool, * "The tool automatically adopts the data frame’s coordinate system"*.  Since I'm importing a csv with no defined coordinate system, I would expect that importing a csv with only lat/long values and a name, that it would produce a data set of points in the origional dataframe's coordinate system at the lat and long given in the table.  That is not what's happening.
 I'm trying to add a CSV table or Excel file (tried both) of a list a radio sites around LA.  The project is NAD83, as I believe that's an FCC standard.  I'm trying to add the table using the ADD XY Data tool, placing the Longitude column in the X field and Latitude in the Y field ( Fig 1).  On the map the points end up half way to Hawaii however.  If you look at the identities window (Fig 2),  you'll see that my lat/longs in the identities table don't match the locations.  For example "Avalon" has a lat of 33.34101 in the table, but 27.196035, and the long is about 20 degrees off also (See green and orange circled values in Fig 2).  In my table below, there are only those values, where you can see where the lat and long values are coming from.  I've also tried a table where i removed the hyper link column, but the result is the same.  Also (fig 4), shows my layer properties for my point data feature.
Does anyone have idea what's going on?
Im thinking the CSV file I'm bringing in has no projection defined.  I would think the lat and long should just line up to the project coordinates.
.

Fig 1

Fig 2

Fig 3

Fig 4


texts: Raster Layer in QGIS - nightlights data not loading

texts:
Using QGIS, I want to load the nightlights data available on https://eogdata.mines.edu/download_dnb_composites.html. My output is an all black screen with numbers shown on the bottom left in the Layers section.
I have tried changing the symbology to all Reds palette but the output is one colour throughout only.
I am not sure what is the problem. Can someone suggest what to do?


texts: Intersect raster or laz (pointcloud) with vector polygons

texts:
I am a vector trained GIS person and wish to add the value of a pointcloud to my vector data. A friend of mine suggested cloudcompare but said it was also possible using my beloved QGIS. I work with postgres in combination with qgis often.

Because of my vector way of thinking i am looking for tools to intersect or add a column to the points with the vector polygon id.

I have the feeling that i am not thinking in the right direction. I am doing this in my free time and it will take a while for me to figure this one out. So i am asking for tips and advise to extract the data of a pointcloud in vector polygons.


texts:
I'm trying to get the difference between two polygon layers using the QGIS difference function. These are the two files:
Dropbox

When I view the attribute tables of both I can see there are 3 polygons that are different between the two. When I run the difference tool I get an output that has fifteen polygons the same as the original file with 15 polygons. Is there something I'm missing? I want to get a separate shapefile with only the 3 polygons that are the difference between the two.


texts: In terms of SQL, what exactly is Query Layer's "Unique Identifier" doing?

texts:
When you make a link with a database using Query Layer in ArcMap, you must select a Unique Identifier, which ensures that the resulting Layer only has unique records in the Attribute table.

In some cases, I have seen used Query Layers in such a way that the Unique Identifier asked by ArcMap is used to make sure that records are distinct (reducing from 4000 to 500 records for instance)

But in terms of SQL what exactly is it doing? Is it a "where" clause, selecting by distinct, selecting the min or the max of what it encounters? I know it's created an ObjectID based on the unique identifiers but I don't know how it's assigning that ObjectID to rows and deciding which rows are dropped from the attribute table.


texts: How to represent large data using Google Map API in JavaScript

texts:
I am trying to represent data using bins as the data was too large to display each and every point simultaneously.

The problem with this also existing if the no of bins are getting larger than 15K, map hangs while zooming in and zooming out.

So I thought of using WMS, I am new to this, I have set up GeoServer, but confused how to add layer to this and how to access that layer.

Can I add layer dynamically to GeoServer as the scenario is: User Uploads a geopoints information file and then views the data in that file.

I am creating a GeoJSON data format from the uploaded file of user. So is it possible to upload that GeoJSON data as a layer at runtime with GeoServer and then ready to view.

Is the huge data representation possible using OpenLayers?


texts:
I am working with some land ownership data. I would like to identify parcels the have small adjacency between other parcels.

In a nutshell I want to identify where there is a small overlap between sections. In the example below you can see where one section overlaps another by a small amount of about 15 meters. In most situations this value of overlap is zero.

Any suggestions how I can detect these situations? I am on arc 10.2 with Spatial Analyst.


texts: numpy array coordinates to polyline in shapefile Python

texts:
I am calculating the shortest path in a numpy array which was converted from a raster (.tif).

I have a list representing a path like [(1,3),(2,3)(3,3)].

How can I transform this into a polyline in a shapefile?

Do I have to use rasters for latitude and longitude and then look up the values at the positions of the path elements or is there an easier solution?


texts: Is it possible to use a WHERE clause with Grass v.in.ogr for DGN sources?

texts: I am trying to import a .DGN with Grass v.in.ogr in QGIS.  I am able to import the DGN perfectly, but I would like to filter the elements using the WHERE option.  From searching, reading and testing I think it should be running correct, but it is erroring out.  Any help would be appreciated.  The attribute column holding the information I would like to Filter on is called Level and it is an Integer. I want to filter it on a value of 1 (see attached image).


texts: My question is whether QGIS can be downloaded and installed on a server for multi-user access. Rather than installing the software locally, can QGIS be stored on a server and accessed by many?


texts: PyQGIS raster image shifted from vector shapefile

texts:
I previously asked a question about my raster image being larger and shifted over from my vector image here, but now I realize that the raster image is only shifted to the left and up from the vector map (both image and map are defined in decimal degree coordinates). The vector layer is a world shapefile and the raster image is a .png image that is overlain on the shapefile.

Below I have two examples where the raster image is shifted when uploaded on the shapefile in PyQGIS:

The blue and green are part of the vector shapefile (blue being the ocean and green being the land). The purple, white grid, and green contours are part of the raster image. The larger red 'X' represents the set of coordinates that were defined on the vector map. The smaller red '+' represents the same set of coordinates, but on the raster image. Here is where I can tell that the raster image is shifted compared to the vector map because the red 'X' and the '+' represent the same coordinates (i.e. 75.23694W, 39.86N for the second image).

I am wondering how to fix this? Is there anyway to easily correct the raster location to the correct area?

My PyQGIS code to upload the raster image is below:

def addRasterImage1(self):
    x = self.lon0
    y = self.lat0
    self.layers = []
    fileName = 'Desktop'

    dir = os.path.dirname(fileName) if fileName is not None else "." #path of the saved raster image
    formats = ["*.%s" % "PNG"]

    file = myfile_list[k_zt]  #raster image to load onto vector image
    print file
    if file:
        layer = QgsRasterLayer(file, fileName)
        print "layer",layer

        if not layer.isValid():
            return

        # add layer to the registry
        QgsMapLayerRegistry.instance().addMapLayer(layer)

        # # set extent to the extent of layer
        self.canvas.setExtent(layer.extent())
        self.layers.insert(0, QgsMapCanvasLayer(layer))

        print "layers number:",len(self.layers)
        self.canvas.setLayerSet(self.layers)
        self.canvas.refresh()
        self.MyWorldLayer(ext = False)

def MyWorldLayer(self, ext=bool):
    file = mypath + "world.shp"
    fileName="world"
    layer = QgsVectorLayer(file, fileName, "ogr")
    if not layer.isValid():
        return

    symbols = layer.rendererV2().symbols()
    symbol = symbols[0]
    symbol.setColor(Color.fromRgb(113,170,0))
    QgsMapLayerRegistry.instance().addMapLayer(layer);

    if ext == True:
        self.canvas.setExtent(layer.extent())
    if ext == False:
        rect = QgsRectangle(xmin,ymin,xmax,ymax)
        self.canvas.setExtent(rect)

    self.layers.insert(0, QgsMapCanvasLayer(layer))
    self.canvas.setLayerSet(self.layers)

Where self.MyWorldLayer(ext = False) is the vector shapefile map.
Please let me know if you have any suggestions on how to fix this!

Update:

I am trying to use this solution here, but this is for a raster of points. Can this be used for a .png raster?

** Here is a link that you can access to see the .png raster and the world shapefile used as the vector layer layers link**


texts: Where can I get the grid representation for a projection on meters?

texts:
In PostGIS I have the following projection:

aut_srid = 21781

text = "PROJCS["CH1903 / LV03",GEOGCS["CH1903",DATUM["CH1903",
SPHEROID["Bessel 1841",6377397.155,299.1528128,AUTHORITY["EPSG","7004"]],
TOWGS84[674.374,15.056,405.346,0,0,0,0],AUTHORITY["EPSG","6149"]],
PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],
UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],
AUTHORITY["EPSG","4149"]],PROJECTION["Hotine_Oblique_Mercator_Azimuth_Center"],
PARAMETER["latitude_of_center",46.95240555555556],
PARAMETER["longitude_of_center",7.439583333333333],
PARAMETER["azimuth",90],PARAMETER["rectified_grid_angle",90],
PARAMETER["scale_factor",1],PARAMETER["false_easting",600000],
PARAMETER["false_northing",200000],UNIT["metre",1,AUTHORITY["EPSG","9001"]],
AXIS["Y",EAST],AXIS["X",NORTH],AUTHORITY["EPSG","21781"]]"

I found this picture in the PDF at the end of the link but I need a map file with the grid representation, something in shapefile, tab or similar.


texts: How do you remove the ArcSDE connection from an MXD entirely?

texts:
I want to remove all connections to the SDE from a number of MXDs that I am using. I have re-set the data sources in each layer in the MXD to a local gdb but still get prompted with the Database connection window asking for user name and password whenever I open the MXD:

It is still possible to open the MXD but I have to go through several instances of selecting Cancel before it opens.
What is still causing this window to appear?


texts:
How can I setup QSWAT in newer versions of QGIS ? The link here says QSWAT requires QGIS 2.6. I tried installing QSWAT in new versions and it couldnot be done.
Is there a way to install QSWAT in new QGIS may be version 3 or so.


texts: Renaming multiple layers in the legend from an attribute in each layer in QGIS

texts:
I have 72 layers that have been generated via affine transform of a single layer that has 72 categories. But the 72 layers have been given names that are sequential numeric. The original categories are named by a field UNIT_NAME. I want the 72 generated layers to have this attribute as a layer name for each of them.
Can I get affine transform function to name the generated layers using the value from the UNIT_NAME attribute of each layer; or can I rename all the layers to have this value, in a script?
I am using QGIS 3.16.1 Hannover on Windows 7.


texts: What does "Feature doesn't have a geometry" mean in QGIS?

texts: When I open the .shp file and click on specfic item in the attribute table, in the canvas, it said "Feature doesn't have a geometry". Is it the same as "geometry is not valid", or does it simply mean there is no related data?


texts: How do I create pyramid tiles from .ecw file?

texts:
I have an .ecw file. I want to pyramid this file and use in my app. How Can I do it with GDAL?

Note: My file is greater then 4 GB.


texts:
I wanted to know if there is a way to split a polygon using JavaScript.
So I have a web map which uses Leaflet API, and the user needs to cut polygons which are shown in the map using a straight line. I was hoping there would be a JavaScript library something JSTS or Turf but if that isn't possible an algorithm would do just fine.


texts: Get the exact extent when exporting a PDF with arcpy

texts:
I'm trying to export to PDF the extent of my selected feature in ArcGIS pro using arcpy.

My code so far:

# Open the aprx project
aprx = arcpy.mp.ArcGISProject(APRX_PATH)

# Get the map
m = aprx.listMaps("*")[0]

# Get the layer
lyr = m.listLayers("MY_LAYER")[0]

# Get the layout
lyt = aprx.listLayouts("*")[0]

# Select some feature in my layer
lyr = arcpy.SelectLayerByAttribute_management(lyr,"NEW_SELECTION","MY_CONDITION")[0]

# Get the mapframe_element
mf = lyt.listElements("mapframe_element")[0]

# Choose the output Height and Width
mf.elementHeight = 11.7
mf.elementWidth = 8.3

# Set the extent
mf.camera.setExtent(mf.getLayerExtent(lyr, True, False))
mf.panToExtent(mf.getLayerExtent(lyr, True, False))

# Clear the selection
arcpy.SelectLayerByAttribute_management(lyr, "CLEAR_SELECTION")

# Export the PDF
mf.exportToPDF(os.path.join(PDF_PATH,'MAP.PDF'),150)

I get something like that

But i need the exact extent

How can I avoid the margin in the extent?


texts: Are there any open source network analysis tools?

texts:
A researcher in my field* recently created an "Urban Network Analysis Toolbox", an ArcGIS extension built upon Network Analyst. It calculates a number of measures on urban networks, such as betweenness and reach.

He was criticized for relying on ArcGIS, a commercial tool not available to many researchers, especially outside the US. But ArcGIS was a convenient platform, since it contains the Network Analyst extension, simplifying the work required for the analysis.

Are there open source network analysis tools that do the work of the Network Analyst extension? Are there other methods for creating the network datasets required as input for the toolbox? My understanding of the toolbox is that it uses the network dataset to calculate a distance matrix for a given set of points on the network. The toolbox does not use any of the higher level problem solving tools in Network Analyst.

*Quantitative Urban Morphology, perhaps? No one seems settled on a name for the broader field. The particular brand of research I do is with Space Syntax.


texts: How to add Raster legends in QGIS map composer?

texts: Is there a plugin or method to add legends for raster layers (with a custom colourmap) as an item in the print composer?


texts: Cut part of raster with another raster (mask)

texts:
I have a raster which I want to cut. A second raster smaller, should be exactly the part where I want to cut the first layer (so I cut the same size as raster 2).

The idea is that raster 2 become the mask for the cut.
How can I do this?


texts: Can I set a geometry using QGIS editing tools for a feature without geometry?

texts: I have a PostGIS point layer with mixed features (with and without geometry). When a non geometry feature is selected from attribute table I don´t know how to set its geometry with QGIS editing tools. The add spatial object is only for new features.


texts: Selecting and Identifying Single Feature using QGIS Plugin?

texts:
given is a layer with some nodes on it and every node has some attributes(name, id, etc...). For showing up the informations of a single node I can use the "Identify Feature Tool". Choose the tool, click the node and see the informations.
I am trying to develop a Plugin with Python, which is doing almost the same: Start the Plugin, click on a node and get the informations about this node(output in the python console is enough for the first steps).
Okay, so this is the status:
I was following this tutorial: http://www.qgisworkshop.org/html/workshop/plugins_tutorial.html which works not for QGIS 2.0.
The 'MouseClickEvent' is working, I get the X,Y-coordinates. But I am not able to select and to get the informations of a node. That is my selectFeature function, which is called by result = QObject.connect(self.clickTool, SIGNAL("canvasClicked(const QgsPoint &, Qt::MouseButton)"), self.selectFeature):

def selectFeature(self, point, button):
    #pntGeom = QgsGeometry.fromPoint(point)
    #pntBuff = pntGeom.buffer( (self.canvas.mapUnitsPerPixel() * 2),0)
    #rect = pntBuff.boundingBox()

    layer = self.canvas.currentLayer()
    selected_features = layer.selectedFeatures()
    for i in selected_features:
        attrs = i.attributeMap()
        for (k,attr) in attrs.iteritems():
            print "%d: %s" % (k, attr.toString())

This let me see nothing. I was reading and trying out a lot of things and at the moment in feel a bit lost. I lost the overview.

Do I need these geometry functions to capture the node by clicking the mouse button?

What functions do I need to get the attributes of a single selected feature?


texts: i'm trying to run "Filter" operation using AcGIS on a stack layer of 331 bands. However, the result is always a raster of only one band. How can i get operation to be done on the whole stack?


texts: Google Maps / StreetView coverage layer in QGIS

texts:
I've come across a useful method for adding Google layers as a raster in QGIS: High resolution, printable alternative to OpenLayers plugin for QGIS?

Does anyone know what the serverUrl is for the google streetview coverage map? I've found this url but I can't get it to work in the XML:

https://mts1.googleapis.com/vt?hl=en-US&lyrs=svv|cb_client:apiv3&style=40,18&x=9&y=15&z=5


texts: Installing module within code with pip and ArcPy?

texts:
I am creating an arcmap add-in which relies on a third party module, usaddress. I am using pip to install the usaddress module on the user's computer. When I run the add in however, it says that usaddress was successfully installed but then fails to import the usaddress module with "No module named usaddress". Once I close arcmap and all python scripts, it then imports the module for some computers. Other computers I have tested, it still doesn't install the module at all. But when I run the code independent of the add-in, it successfully installs the module. This is not ideal as I can't do it on every user's machine.

Any ideas as to why this is happening?

See the relevant code below:

import pip
import arcpy
required_pkgs=["usaddress"]
installed_pkgs=[pkg.key for pkg in pip.get_installed_distributions()]
for package in required_pkgs:
    if package not in installed_pkgs:
        arcpy.AddMessage("Installing {0}..................".format(package))
        pip.main(['install',package])
        arcpy.AddMessage("Succcessfully installed {0}..................".format(package))
    else:
        arcpy.AddMessage("Upgrading {0}..................".format(package))
        pip.main(['install',"--upgrade", package])
import usaddress


texts: Buffer with fixed distance create oval, not circular buffer

texts:
Using QGIS 2.16.3, I wish to create a buffer zone around a single point in my multiple point file with a fixed radius of 5 km. Therefore, I made a selection and created a new layer only containing this point.
Next, I used Vecor -> Geoprocessing Tools to create a buffer zone with a fixed radius. Interestingly, the result is an oval buffer zone instead of a circular one, thus the radius is clearly not fixed, but flexible. What do I have to alter to make the bufer zone circular?

I have attached pictures of my input in the buffer zone tool, as well as a picture of my result.

Tool entry


texts:
What is a good R equivalent of the GRASS function r.grow?

The closest I've been able to find is the script suggested by whuber in response to Creating randomly shaped clumps of cells in a raster from seeds of 1 cell/pixel?, which is written to address a different problem, so I'd like to know whether there are more efficient tools to simply grow outward uniformly.


texts: Converting world informations (PGW file - pixels) to geographical coordinates?

texts:
I'm trying to load an image on a Leaflet map and to display it at the right place using the world information exported by QGIS.

For image overlays, Leaflet needs the coordinates of the upper left point and of the lower right point (boundaries):

var imageUrl = 'http://www.lib.utexas.edu/maps/historical/newark_nj_1922.jpg',
imageBounds = [[40.712216, -74.22655], [40.773941, -74.12544]]; L.imageOverlay(imageUrl, imageBounds).addTo(map);

In QGIS, when I export my image, I always use the 1:25000 scale and 72 DPI, which is a constant in my workflow.

It outputs a world file (PGW) that contains the information:

8.86711729729699094 //xscale
0 //yskew
0 //xskew
-8.86711729729699094 //yscale
491355.98836270271567628 //xpos
6600748.93414135184139013 //ypos

Knowing all that information, is there a way to convert the world file information into geographical coordinates I could use in Leaflet?


texts:
I'm looking for a downloadable listing of place names for Iran, in Farsi script.

I downloaded GeoNames for Iran, but many of the place names in there seem to be only in Latin script (of 253,000 names, 100,000 are in Farsi script).
OSM has been suggested by Brad, to be queried using Overpass (Turbo). But I'm not sure whether it's possible to extract places in a particular country. Some tags like is_in or is_in:country are helpful, but are not consistently used. Maybe there's a relation I need to find?
Global Gazetteer from Falling Rain Genomics seems to have lots of data from several public sources, including 200,000 place names in Roman script and 200,000 in Farsi and/or Arabic. But I don't see any way to extract the data other than web-crawling/scraping.


texts:
I am using QGIS3
I have GeoSON layer that is not showing up on my canvas.
When I zoom to the layer it zooms to a point that is very far from the rest of my map features.
I would like the 'browfeild' layer to be in line with the other features which are where the red dot is on the satellite map !area where browfeild layer should be]2
'Browfeild' layer properties are as follows:
Name
browfeild copy
Path
/Users/Joan-dark/Downloads/browfeild.geojson
Storage
GeoJSON
Comment
Encoding
UTF-8
Geometry
Polygon (MultiPolygon)
CRS
EPSG:4326 - WGS 84 - Geographic
Extent
504833.7039953188505024,158488.5467496835626662 : 559097.5517462189309299,199856.3195061662117951
Unit
degrees
Feature count
3,017
Extent
CRS
Spatial Extent
CRS: - Projected
X Minimum: 0
Y Minimum: 0
X Maximum: 0
Y Maximum: 0
Temporal Extent
Start:
End:
appreciate any help. thanks


texts:
Problem: I have a set of LAS files that doesn't display in Arcmap.

Details: I try to use LAS to tin in Arcmap to convert my LAS points to a TIN. However, when browsing to the folder with the files, it doesn't recognize them as LAS files: they do not show up even if I choose *All.

I have tried a couple of conversion tools available in LAStools, but without success. I have also tried to use FME to convert them to shapefiles, but the only thing I get is the bounding rectangle.

I seem to remember hearing something about this being a problem in Arcmap 10.1 with data that has been extracted using FME (which I'm not sure these are, but still), but now I can't find  anything about it online.

Has anyone else encountered this? Any possible solutions, workarounds, or at least a bug report?

To answer @Andre Silva's comment: I'm trying to convert the point cloud to a raster to do hydrological analyses with it, e.g. calculate catchments etc. I have used other las files without issues before, so my interpretation is that this problem is specific to these files. As to the exact details of the files, they're <50 MB .las files that's been filtered for ground points only.


texts:
I have point geometry in PostgreSQL/PostGIS which I query as follows:

SELECT ST_AsText(geom) as pointgeom
FROM pointable.

I'm using Python as server side script.

I want to pass the result as GeoJSON response to overlay on Google Maps but I'm new to GeoJSON.

If this make sense, can any one shade light on the code that encode the pointgeom to GeoJSON format and handover to JavaScript?


texts:
When dealing with MODIS swath products (or any swath products from VIIRS, SeaWifs, etc) such as those available here, is there any way of a-priory knowing the location a given file covers?

Typical filenames are for example T2016009004000.L1A_LAC.bz2 (Acquired on day 9 of 2016 at 0.40 AM)


texts:
I am trying to insert a new polyline dataset within my Oracle server database. The dataset in question is located within a local GDB, and is currently c420mb in size. It is taking hours and hours to run the Feature Class to Feature Class tool to upload the dataset to the SDE

Whilst composing this question, the tool completed the upload, but took took 3.5 hours using the normal settings.

Previously, when this type of situation has occurred before, I have been told to set the Configuration Keyword as "SDELOB" before running the tool. This usually works and uploads the dataset within seconds rather than hours.

So why is there a difference in upload performance? What do the keywords actually mean, and when am I suppose to use them, if at all?


texts: overpass query took long time and no response

texts:
I have a list with 900 (lat,lon) sets.
I used a loop to query each (lat,lon).
It stuck on 3rd set (48.55888888888889,12.404166666666667) and had no response.

My language is Java and reference this example(example).

The stuck line is HttpResponse response = getHttpClient().execute(post);

below is my query

[out:json][timeout:100];(way(around:5000,lat,lon)["waterway"="river"];way(around:5000,lat,lon)["waterway"="riverbank"];way(around:5000,lat,lon)["waterway"="waterfall"];way(around:5000,lat,lon)["tunnel"="*"];way(around:5000,lat,lon)["landuse"="forest"];way(around:5000,lat,lon)["natural"="water"];way(around:5000,lat,lon)["natural"="coastline"];way(around:5000,lat,lon)["natural"="geyser"];way(around:5000,lat,lon)["natural"="bay"];way(around:5000,lat,lon)["natural"="beach"];);
(._;>;);out body;

It can work on overpass turbo and take about 5 secs.
overpass_query

could any one give me some suggestions?


texts:
How can I change the format of the scaleline in OpenLayers? I have a black background and I want to change the lines in the scalebar to white. I use:

<link rel="stylesheet" href="http://www.openlayers.org/dev/theme/default/style.css" type="text/css" />
 <style type="text/css">
     div.olControlScaleLine{
        color: white;
     }

But this only change the font. Anybody kmow how?


texts: Copy and paste a feature in QGIS at the current mouse position (paste in place): not on top of the feature copied

texts:
Is there a way to copy a feature and paste it to a new position using the mouse?

At the moment, I have to do this in two steps:
 1. Paste on top of the existing, copied feature and then
 2. Move it to a new location

This is double handling (and pasting directly on top of an existing feature primes the layer to have duplicate geometries if not especially careful to move the pasted feature straight away).


texts: Where vrt files store the location of the tiles

texts: I have thousand of tif tiles and a vrt file, generated by using gdalbuildvrt. At the moment they are at the same folder level in the system, but if I open the vrt file with Notepad I cannot find any url, address or location. Does this mean that the vrt file and the tiles cannot be separated from each other and have to be together in the same folder? how can I edit that in the vrt file?


texts:
I am getting "ResultObject: Error in getting output" message when I run the following code. I know it has something to do with "os.path.join" but can't figure out what.

import arcpy
import os

    arcpy.env.overwriteOutput = True
    CADInput = arcpy.GetParameterAsText(0)
    OutputFolder = arcpy.GetParameterAsText(1)
    OutputGDB = arcpy.CreateFileGDB_management(OutputFolder, "AllPipes.gdb", "CURRENT")
    SpatialReference = r"C:\Program Files (x86)\ArcGIS\Desktop10.0\Coordinate Systems\Projected Coordinate Systems\UTM\NAD 1983\NAD 1983 UTM Zone 17N.prj"
    OutputData = arcpy.CreateFeatureDataset_management(OutputGDB, "Data", SpatialReference)
    OutputCAD = arcpy.CADToGeodatabase_conversion(CADInput,OutputGDB,"CAD_Converted","1000",SpatialReference)
    OutputCADLocation = OutputCAD.getOutput(0)

    layerfld = arcpy.AddFieldDelimiters(os.path.join(OutputCADLocation,"Point"),"Layer")
    query = layerfld + "='WM-VALVE-ID'"
    arcpy.MakeFeatureLayer_management(os.path.join(OutputCADLocation,"Point"),"Valve_lyr8",query)
    outfile = os.path.join(OutputData,"Valve")
    ValveFC = arcpy.CopyFeatures_management("Valve_lyr8",outfile)
    arcpy.Delete_management("Valve_lyr8")


texts:
I'm working on the same theme as User 32882, and I'm new here.
I have the same problem, and update his question with more details.
Old question: "I'm running into problems running the r.neighbors GRASS command in QGIS. Although it will run for small neighborhoods as soon as I increase my neighborhood above a certain number it becomes unbearably slow. I don't understand why Focal Statistics on ArcMap will do this in the blink of an eye while R.Neighbors takes ages?"
Link of Old question:
R.Neighbors with large neighborhoods?
Details:
We are exploring a circular neighborhood of ~200 pixels.
And trying to compute Maximum or Minimum statistics.
Moreover, Maximum computation seems OK, but Minimum often presents 0 results.
Why is so slow in QGIS?
What can I do to make it run faster?


texts: Unable to delete .lock files within geodatabase using Python

texts:
I am having problems with lock files while deleting and copying geodatabases for backup purposes.

print "Now backing up your Geodatabase"
distutils.dir_util.remove_tree (r"\\Msukserver\gis\geodatabase BACKUP\EVERYONE\Mining Features (MATTHEW).gdb", 0, 0) # deletes 'everyone' backup copy
distutils.dir_util.copy_tree(r"\\Msukserver\gis\geodatabase BACKUP\MATT\Mining Features (MATTHEW).gdb", r"\\Msukserver\gis\geodatabase BACKUP\EVERYONE\Mining Features (MATTHEW).gdb") # backup the last backup to 'everyone'
print "old backup copied"
distutils.dir_util.remove_tree(r"\\Msukserver\gis\geodatabase BACKUP\MATT\Mining Features (MATTHEW).gdb", 0, 0) # deletes last backup
print "old backup deleted"
distutils.dir_util.copy_tree(r"C:\GIS Home\Mining Features (MATTHEW).gdb", r"\\Msukserver\gis\geodatabase BACKUP\MATT\Mining Features (MATTHEW).gdb") # copy mining features to create new backup
print "::::::::::::::::::::::Backup Complete::::::::::::::::::::::"
raw_input("Push Return key to close...")

I am getting a 'permission denied' error on the '.lock' files. Is there a way of copying/deleting all the file except those with a .lock extension? or a way to continue/force with errors? I am aware that I am currently copying/deleting the entire directory and may need to change my code to copy/delete files only.


texts: Raster calculation in for loop returning exit code -1

texts:
I've put together a python script that does some calculations on rasters by converting them to NumPy arrays, doing the math, then converting them back to Raster and applying all the same spatial reference information, etc.

I just need to run it iteratively now on about 400 rasters.

I tried to use glob.glob to iterate through the rasters with a *tif wildcard, but couldn't get the output name to work right.

Now I'm trying to use os.path.basename with a conditional statement to filter through the files in the folder and only run my raster calculation on tif files.

When I try to run the script now, I get a "Process finished with exit code -1" message, and my print(outname) statement doesn't seem to register.

Apologies in advance if this is more of a python than an ArcPy question.

Also, thanks for your help.

import os
import arcpy
import numpy as np

#Enable overwriting output data arcpy.env.overwriteOutput = True

def cumulativecalculation(filename):

#Set geoprocessing variables
inraster = filename
des = arcpy.Describe(inraster)
sr = des.SpatialReference
ext = des.Extent
ll = arcpy.Point(ext.XMin, ext.YMin)

#Convert GeoTIFF to numpy array
a = arcpy.RasterToNumPyArray(inraster)

#Flatten for calculations
a.flatten()

#Find unique values, and record their indices to a separate object
a_unq, a_inv = np.unique(a, return_inverse=True)

#Count occurences of array indices
a_cnt = np.bincount(a_inv)

#Cumulatively sum the unique values multiplied by the number of
#occurences, arrange sums as initial array
b = np.cumsum(a_unq * a_cnt)[a_inv]

#Divide all values by 10 (reverses earlier multiplication done to
#facilitate accurate translation of ASCII scientific notation
#values < 1 to array)
b /= 10

#Rescale values between 1 and 100
maxval = np.amax(b)
b /= maxval
b *= 100

#Restore flattened array to shape of initial array
c = b.reshape(a.shape)

#Convert the array back to raster format
outraster = arcpy.NumPyArrayToRaster(c, ll, des.meanCellWidth, des.meanCellHeight)

#Set output projection to match input
arcpy.DefineProjection_management(outraster, sr)

#Setting the OutName
outname = filename + "_cumulative" + ".tif"

#Save the raster as a TIFF
outraster.save("E:/NSF Project/Salamander_Data/New_Cumulative_Output/" + outname)
print (outname)

#Loop through files in directory

for filename in os.path.basename("E:/NSF
Project/Salamander_Data/NoDataToZero/HadleyGCM/"):
   if filename[-4:] == ".tif":
       cumulativecalculation(filename)

Updated Code:

    #Set the output filename
    intername1 = os.path.splitext(os.path.basename(i))[0]
    filename = intername1[5:]
    finalname = filename + "_cumulative" + ".tif"

    #Save the raster as a TIFF
    outraster.save("E:\\NSF Project\\Salamander_Data\\New_Cumulative_Rasters\\" +    str(finalname))

    #Announce file creation
    print("Created: " + str(finalname))

#Loop through files in directory
geotiffs = glob.glob("E:\\NSF Project\\Salamander_Data\\NoDataToZero\\HadleyGCM\\*tif")

for i in geotiffs:
    cumulativecalculation(i)

Everything seems to be working beautifully with the code now, except it "breaks" after processing 10 rasters.  I think this must be a problem with one of the input files, but I'm not sure.

Error: "File "C:\Program Files (x86)\ArcGIS\Desktop10.2\arcpy\arcpy__init__.py", line 1860, in NumPyArrayToRaster
    return _NumPyArrayToRaster(*args, **kwargs)
TypeError: Cannot create raster for numpy array."

I have about 400 to get through, so I'd at the very least like to keep errors like this from derailing the process, even if I end up with some missing data.  Not sure how to accomplish that though.


texts:
I am using Network Analyst to run some analyses on global networks.  In particular, I am trying to find shortest distances between points.  However, I noticed that routes are not crossing the Pacific.  There must be some disjoint at each edge of the network, rather than connecting from one side to the other.  I assume the problem stems from representing a 3D problem in 2D.

As an example, a route from Vancouver to Japan will go south, through the Panama Canal, across the Atlantic, around the southern tip of Africa and then diagonal up to Japan.  Of course, the shortest path is directly across the Pacific.

Any idea of how to remedy this issue?


texts: How to choose cell size for grid points to fill polygons in R?

texts:
I have a shapefile consisting of hundreds of polygons depicting building footprints in a city. I want to fill these polygons with grid points. I am doing this by making a grid of points and then clipping the points grid to the polygons so that I can get a new clipped layer that is the grid points but with the invisible boundaries of the polygons.
Here is my code:
library(sp)
library(rgdal)
library(raster)
library(rgeos)

polygons <- readOGR("C:/PATH_TO_MY_FILE/POLYGONS.shp")
polygons <- gBuffer(polygons, width=0, byid=TRUE)
grid <- makegrid(polygons)
grid <- SpatialPoints(grid, proj4string = CRS(proj4string(polygons)))
grid <- spTransform(grid, CRSobj = crs(polygons))
grid.clip <- grid.t[polygons, ]
grid.clip

This results in the following:
class       : SpatialPoints
features    : 0
crs         : +proj=lcc +lat_0=40.1666666666667 +lon_0=-74 +lat_1=41.0333333333333 +lat_2=40.6666666666667 +x_0=300000 +y_0=0 +datum=NAD83 +units=us-ft +no_defs

Running plot(grid.clip) yields nothing but blankness.
To identify the CRS of the polygons:
crs(polygons)

Which yields:
CRS arguments:
 +proj=lcc +lat_0=40.1666666666667 +lon_0=-74 +lat_1=41.0333333333333 +lat_2=40.6666666666667 +x_0=300000 +y_0=0 +datum=NAD83
+units=us-ft +no_defs

I am not sure why I am getting no features from the clipped grid points, but my guess is because there were no points identified within any of the polygons. According to the following picture of my data, it looks like this is because the dimensions of the grid are significantly far off from the scale of the polygons.
I run the following to see what I was working with before the clip:
library(mapview)
mapview(polygons) + mapview(grid)

And I see this:

As you can see, the points are spread out so far beyond the scale of the polygons. How can I shrink down the spacing between the grid points so that I can actually see the points within the polygons? I recognize that there is no direct answer for what the scale should be, but I am looking for something reasonable. Since these are building footprints, and since this polygon layer has units in feet, I was thinking that maybe the grid cell size should be within feet. How do I set this? I think this involves the cell.size argument within makegrid(), but I am not sure how to use this. Also, what might be a good way to decide between 1-foot spacing vs 5-foot spacing vs 10-foot spacing for the grid points? I want the best resolution/highest reasonable point density without overloading RStudio and crashing.


texts: Error post installation ArcSDE for SQL server 2008

texts:
I had a problem during post installation of ArcSDE for SQL server 2008,
and this is the message :

  The ArcSDE Repository was unsecessfully completed. would you like to
  view this error?

then when I click on OK, what appears is:

  Connection to sde failed(0)

please can somebody help me ?

I am using ArcGIS 10,ArcSDE 10,SQL Server 2008 and for service pack I use SP2


texts:
Background:

At 3.10.5, I have a Geopackage point layer that I wish to export as a KML, and then view in Google Earth.  The point layer has a single integer field, id_unique, that needs to display as a label in GE.

I use QGIS "Export > Save Feature As..." to create the KML.

Problem:

I am unable to export a .kml that displays the label from the id_unique field.  I have tried all combinations of options in the "Save Features As..." window.  When opened in the GE table of contents, the resulting points are always named [no_name], and the map push-pins are not labeled, as shown below.

How do I export a KML that will display a field attribute as a GE label?


texts: PostGis: Calculating distance in different methods, What's the difference?

texts:
What is the most efficient (or just "best practice") way to calculate distance between 2 points, assuming we are on 2D world without big circle etc. phenomenon (distances up to 300 meters).

Method 1:

SELECT ST_Distance(ST_GeomFromText('POINT(lon1 lat1)',4326)::GEOGRAPHY, ST_GeomFromText('POINT(lon2 lat2)',4326)::GEOGRAPHY);

Method 2:

SELECT ST_Distance(ST_Transform(ST_GeomFromText('POINT(lon1 lat1)',4326), 26986), ST_Transform(ST_GeomFromText('POINT(lon2 lat2)',4326), 26986));

assuming 26986 is the correct metric SRID for the location.

Method 3:

Using ST_Distance_Spheroid


texts: Using CQL_FILTER BBOX on WMS layer causes ORA exception window SRID does not match layer SRID

texts:
I'm trying to limit WMS layer to a specific BBOX.

mainLayer.mergeNewParams({
    'CQL_FILTER' : "BBOX(geom, 6541844, 5577147, 6584851, 5620154)"
});

Which gives following exception in geoserver's log:

Internal error while evaluating [window SRID does not match layer SRID] operator ORA-06512: at "MDSYS.SDO_INDEX_METHOD_10I", line 286.

Obviously geoserver hasn't pass SRID information of my BBOX to an oracle sdo filter.
How can I deal with that? Should I edit Geoserver source or is there any Openlayers way to alter "window" SRID in cql_filter query?

Openlayers version: 2.11, Geoserver version: 2.1.2

This question is related with another question: Restrict CQL filter to an specific BBOX


texts: How to add footprint parameter into ODATA API query (Sentinel-hub search)?

texts:
As mentioned in the FAQ about OData API usage I can perform such queries to search Sentinel-2 scenes:

https://scihub.copernicus.eu/dhus/odata/v1/Products?$format=json&$filter=startswith(Name,'S2') and ContentDate/Start ge datetime'2018-08-25T00:00:00' and ContentDate/End le datetime'2018-08-26T00:00:00'

But I still don’t understand how to add footprint search parameter into this query, because there is no clean examples in the FAQ, only mentioned /Products properties:

  ContentGeometry   footprint polygon coordinates in GML format

I am asking for few examples how to query Sentinel-2 scenes using point/polygon coordinates with ODATA API.


texts: Python: Running a function on all points in a list

texts:
I'm trying to get a list of cartesian points along contour lines and convert them to polar coordinates. I'm trying to run a function to convert them to polar, one for the length and one for the angle. However when I run the function it only results in one output, so it only runs on one point in the list. Here is the code:

def Polar(x1,y1):
    return math.sqrt(x1**2 + y1**2)
# ---------------------------------------pnt only pulls one point (x,y)
for pnt in part_list:
    print Polar(pnt[0] ,pnt[1])

polarcoord= Polar(pnt[0], pnt[1])

def pheta(x2,y2):
    return math.atan (y2/x2)

for pnt in part_list:
    print pheta(pnt[1], pnt[0])

thetacoord = pheta(pnt[1], pnt[0])

thetafinal = math.degrees(thetacoord)

The part_list is the list of Cartesian points from the contour line. The list is like 20,000 points long or so. I'm wondering how to get the polar coordinate (length) and the thetafinal (angle) for all the points in the list. It should be noted I am very new to Python.


texts: How to submit user scripts for Processing Toolbox?

texts: QGIS 2.6 change-log says that now there is some kind of repository where users can submit their scripts to make them available for everyone. Unfortunately I was unable to find any details about how to do this.


texts: Downloading sidewalks as shapefile from OpenStreetMap?

texts:
How can I download the sidewalks as a layer from OpenStreetMap (OSM)?

I exported an OSM map and opened it in QGIS. I see "map lines", "map multistrings". But I don't see sidewalks as a layer.


texts: Create a raster image with integer sequence in Google Earth Engine

texts:
I'd like to create an image with pixel values that are integer sequence 1,2,3,4.... . With pixel size equal to a MODIS pixel 500m and pixel align with MODIS as well. So that I can use a grouped reducer to do zonal statistics to count how many Landsat pixels (30m) are included in a MODIS pixel(500m). What I've tried to addressing this issue are:

to generate a grid with size equal to MODIS image resolution. WORKED but facing pixel alignment issues or projection issues which I will post for another question link. I assume by generating an image will be more straightforward to do zonal statistics than create a grid polygon in Arcmap then imported to GEE.
convert MODIS image to vector to create a pixel-to-polygon layer and use polygon ID to zonal. Failed. The converted pixel polygon is not corresponding to pixel in an image, somehow the algorithm turn two neighbor pixels as one single polygon even I set the eightConnected argument in reduceTovectors as false. This is not desired.


texts:
I have a table with lat/long values in MySQL.
The sequence column gives me the order in which the path must be traced.
.
I need to create a line from these values and display the flight path via a WMS service using GeoServer.
I am making use of a SQL view with the following query to generate the paths:

SELECT  GeomFromText(CONCAT('LINESTRING (',
GROUP_CONCAT(CONCAT(longitude, ' ', latitude) ORDER BY sequence SEPARATOR ','),
')'
)) as flight_path
FROM
flight_table where id = 2345

This works well and I am able to display the paths.
However when I have a flight path that crosses the international date line then I get the following output

What would be the ideal way to resolve this in GeoServer?


texts: RuntimeError occurs randomly on pyqgis loop

texts:
I wrotee a loop in pyqgis 3 to apply processing to each polygon of a layer.
Apparently randomly, a RuntimeError occurs and crashes my script.

Here is the loop part of my script, with the first processing :

for feature in features :
    feature_id = feature.id()
    nonbati.selectByIds([feature_id])
    temp_layer = QgsVectorLayer('MultiPolygon?crs=EPSG:2154', "temporary_polygons", "memory")
    crs = temp_layer.crs()
    crs.createFromId(2154)
    temp_layer.setCrs(crs)
    selection = nonbati.selectedFeatures()
    temp_data = temp_layer.dataProvider()
    attr = nonbati.dataProvider().fields().toList()
    temp_data.addAttributes(attr)
    temp_layer.updateFields()
    temp_data.addFeatures(selection)

    # Tampon
    alg_params = {
        'DISSOLVE': False,
        'DISTANCE': 0.5,
        'END_CAP_STYLE': 0,
        'INPUT': temp_layer,
        'JOIN_STYLE': 0,
        'MITER_LIMIT': 2,
        'SEGMENTS': 5,
        'OUTPUT': QgsProcessing.TEMPORARY_OUTPUT
    }
    outputs['Tampon'] = processing.run('native:buffer', alg_params, context=context, feedback=feedback, is_child_algorithm=True)

Error always occur on the 10th line :

attr = nonbati.dataProvider().fields().toList()

Here is the error message :

Traceback (most recent call last):
File "C:\Users\moulin_z\AppData\Roaming\QGIS\QGIS3\profiles\default\processing\scripts\iteration_while_1706_fonctionne.py", line 63, in processAlgorithm
attr = nonbati.dataProvider().fields().toList()
RuntimeError: wrapped C/C++ object of type QgsVectorDataProvider has been deleted

I don't understand why this error sometimes occurs, or how fix it.


texts: How to make a minimal distance matrix for polygons using QGIS?

texts:
I have a set of polygons and I want to find minimal distances between the edge of each pair of them. It's similar task as in this question:

Finding minimum edge to edge distance of polygons using ArcGIS Desktop?

But I'd prefer to make it in QGIS or other free available program. Also I'm very beginner in using geographic information systems software so I'd be really grateful if you make your answer as simple as possible.


texts: Is it possible that QGIS is truncating my distance matrix?

texts:
I need to calculate the pairwise distances among 330 points that are in one shapefile and also pairwise distances among 948 points that are also in one shapefile. These points are in shapefiles, each has latitude/longitude coordinates that are in the projection I need to obtain distances in meters.  I want the output in list form (n*kx3).

The Distance Matrix operation seemed to work perfectly until I noticed that the two output files based on the different numbers of points both contained 107,910 lines.

The 330 point set should have 54,285 distances and the 948 point set should have 448,878 if only one unique distances are output.  However, QGIS is giving both "from A to B" and "from B to A".  Given that, the smaller data set should have 108,570 lines and the larger 897,756.  Thus, it appears that the output is being truncated.  It happens even if I do not open the files in Excel but rather bring them right into R using read.csv. Everything I have read indicates that csv files have no line limits and Excel now has a limit over 1 million lines.

The problem was occurring in Version 2.6.1 and when I realized how far behind I was, I updated to Version 2.10.1 thinking that might help my problem.  It made it worse.  Now my 330 point data set runs and produces the same truncated 107,910 pairs of distances but the 948 point set crashes QGIS.

Is there a way around this or is there something that I am missing/doing that is causing the truncation and crashing?


texts: $area is not returning the correct value in QGIS

texts: I'm working in a black QGIS project set to EPSG:27700 and a layer in EPSG:2770. When I measure the area of a polygon using the measure tool I get a value of 7594315ha. When I subsequently do this in the field calculator using $area/10000 to return the value in hectares I get a value of 3. I'm not sure why this is so small as I've never had an issue with this before and both project and layer are in the same projection.


texts: Cannot intercept Alt button when user drags a box (OpenLayers)

texts:
I am using OpenLayers.
When a user holds Shift, clicks on map and moves the mouse, the box is shown.
When user releases the mouse button, I intercept the event and use the coordinates of the box.

        var boxControl = new OpenLayers.Control();
        OpenLayers.Util.extend(boxControl, {
            draw: function () {
                // this Handler.Box will intercept the shift-mousedown
                // before Control.MouseDefault gets to see it
                this.box = new OpenLayers.Handler.Box(boxControl,
                        { "done": this.notice },
                    { keyMask: OpenLayers.Handler.MOD_SHIFT });
                this.box.activate();
            },

            notice: function (bounds) {

                var ll = map.getLonLatFromPixel(new OpenLayers.Pixel(bounds.left, bounds.bottom));
                var ur = map.getLonLatFromPixel(new OpenLayers.Pixel(bounds.right, bounds.top));
                doSomething(ll, ur);
            }
        });
        map.addControl(boxControl);

It works when with Shift button. But when I change MOD_SHIFT to MOD_ALT, my code is ignored, the bounding box is never shown on screen and "notice" function is never called.

Could you help to fix it (for either Alt or Control buttons), please?

PS. I cannot use OpenLayers.Control.Select


texts: How to work with the STRtree Index when object in the index change it's position?

texts:
I'm work with the NTS(NetTopologySuite), I get large amounts of Point geometries, for fast querying the point object, I make the STRtree Index like this,

STRtree tree = new STRtree();
foreach (var p in allPoint)
{
    tree.Insert(p.EnvelopeInternal, p);
}

While some point objects in the "allPoint" change it's position.
And I find I can't query the point which change position anymore..

How to work with it?

Thanks.


texts: Last edited information missing when used QGIS - Dxf to Shp

texts:
I am doing as constructed drawings which require the information to be in a GIS system and more specifically drainage as constructed drawings.

I must have a separate file for Pipes and Pits. First I use Autocad, save as a DXF then use Quantum GIS to convert the DXF to SHP. So the problem is that the last edited piece of information from the autocad file is missing when i convert it. I can go back and change the last edited line, polygon or polyline and the last thing touched will still be missing. Help?

I have looked on the QGIS website but haven't been able to find anything helpful.


texts: Suggested ways to do a Google-Earth-like fly-to-the-location effect with OpenLayers?

texts: What are the suggested ways to do a fly-to-the-location effect, similar to the one when opening Google Earth, with OpenLayers?  While it might be possible to draw this with OpenLayers API, we would like to know if there are better (smoother and slicker) ways to do this.